void plot_gauss(void)
{
// 
// plot various probability density distributions
//

// #include <TMath.h>
#include <TRandom.h>

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kFALSE);
// gStyle->SetOptFit(kTRUE);
// gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
gStyle->SetFillColor(0);
//

   char string[256];
   char filename[80];
   Int_t j,jj;
   
   #define npts 100000;
   # define nbins 100;

   // define parameters

   Double_t mean=6.1;
   Double_t sigma=0.5;
   Double_t norm=1;

   // define limits

   Double_t xmin=mean-5*sigma;
   Double_t xmax=mean+5*sigma;

   sprintf (string,"Gaussian: mean=%f, sigma=%f, norm=%f\n",mean,sigma,norm);
   printf ("plot_gauss: %s",string);

   TF1 *func = new TF1("func",func,xmin,xmax,3);
   func->SetParameters(mean,sigma,norm);


   TF1 *intfunc = new TF1("intfunc",intfunc,xmin,xmax,3);
   intfunc->SetParameters(mean,sigma,norm);


   TCanvas *c1 = new TCanvas("c1","c1 Plot PDFs",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->SetGridx();
   c1->SetGridy();
   // c1->SetLogy();

   Double_t ymin = 0;
   Double_t ymax = 1;
   intfunc->GetXaxis()->SetRangeUser(xmin,xmax);
   intfunc->GetYaxis()->SetRangeUser(ymin,ymax);
   intfunc->GetXaxis()->SetTitleSize(0.05);
   intfunc->GetYaxis()->SetTitleSize(0.05);
   intfunc->GetYaxis()->SetTitleOffset(1.5);
   intfunc->GetXaxis()->SetTitle("variable x");
   intfunc->GetYaxis()->SetTitle("probability");
   intfunc->GetXaxis()->SetNdivisions(5);
   // intfunc->SetMarkerColor(1);
   // intfunc->SetMarkerStyle(21);
   intfunc->SetTitle("");
   // intfunc->Draw("");
   // intfunc->SetLineColor(2);
   func->GetXaxis()->SetRangeUser(xmin,xmax);
   func->GetYaxis()->SetRangeUser(ymin,ymax);
   func->GetXaxis()->SetTitleSize(0.05);
   func->GetYaxis()->SetTitleSize(0.05);
   func->GetYaxis()->SetTitleOffset(1.5);
   func->GetXaxis()->SetTitle("variable y");
   func->GetYaxis()->SetTitle("probability density");
   func->GetXaxis()->SetNdivisions(5);
   func->SetTitle("");
   func->SetTitle("");
   func->SetLineColor(4);
   // func->DrawCopy("same");
   func->Draw();

   sprintf(string,"#mu = %.1f\n",mean);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"#sigma = %.1f\n",sigma);
   t1 = new TLatex(0.2,0.80,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"y = 5.6 #pm 0.5\n");
   t1 = new TLatex(0.2,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   // TLine *lmean = new TLine(mean,ymin,mean,ymax);
   TLine *lmean = new TLine(5.6,ymin,5.6,ymax);
   // lmean->SetNDC();
   lmean->Draw();
   TLine *lsig1 = new TLine(mean-sigma,ymin,mean-sigma,ymax);
   // lsig1->Draw();
   TLine *lsig2 = new TLine(mean+sigma,ymin,mean+sigma,ymax);
   // lsig2->Draw();
    
   

   //
   TCanvas *c3 = new TCanvas("c3","c3 Plot PDFs",200,10,700,700);
   c3->SetBorderMode(0);
   c3->SetFillColor(0);
   //c3->SetGridx();
   //c3->SetGridy();
   //c3->SetLogy();

   TH1F *pdf = new TH1F("pdf","random generation",nbins,xmin,xmax);

   for (j=0;j<npts;j++) {
      Double_t rand = func->GetRandom();
      pdf->Fill(rand);
      }

   // pdf->Fit("func","","",xmin,xmax);

   pdf->SetTitle("");
   // pdf->GetXaxis()->SetRangeUser(1.0,2.5);
   // pdf->GetYaxis()->SetRangeUser(ymin,ymax);
   pdf->GetXaxis()->SetTitleSize(0.05);
   pdf->GetYaxis()->SetTitleSize(0.05);
   pdf->GetYaxis()->SetTitleOffset(1.5);
   pdf->GetYaxis()->SetTitle("events per bin");
   pdf->GetXaxis()->SetTitle("random variable x");
   pdf->GetXaxis()->SetNdivisions(505);
   pdf->SetMarkerColor(2);
   pdf->SetMarkerStyle(21);
   // pdf->Draw();

   sprintf(string,"Entries = %.0f  nbins = %.0f  #Delta x = %.1f\n",pdf->GetEntries(),(Float_t )nbins,(xmax-xmin)/nbins);
   t1 = new TLatex(0.15,0.92,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c1->SaveAs("plot_gauss_c1.png");
   c3->SaveAs("plot_gauss_c3.png");
     
}
Double_t func(Double_t *x, Double_t *par) {

   // function to be plotted

   Double_t p1=par[0];
   Double_t p2=par[1];
   Double_t p3=par[2];
   Double_t x1=x[0];

   char string[256];
   Double_t pi=3.14159;

   Double_t distr = exp(-(x1-p1)*(x1-p1)/(2*p2*p2))/(sqrt(2*pi)*p2);

   sprintf(string,"p1=%g, p2=%g, p3=%g, x=%g, distr=%g\n",p1,p2,p3,x1,distr);
   // printf ("func: %s",string);

   return p3*distr;
}
Double_t intfunc(Double_t *x, Double_t *par) {

   // function to be plotted

   Double_t p1=par[0];
   Double_t p2=par[1];
   Double_t p3=par[2];
   Double_t x1=x[0];

   char string[256];
   Double_t pi=3.14159;

   Double_t arg = (x1-p1)/p2;

   Double_t distr = 0.5 + 0.5*TMath::Erf(arg/sqrt(2));

   sprintf(string,"p1=%g, p2=%g, p3=%g, x=%g, distr=%g\n",p1,p2,p3,x1,distr);
   //printf ("func: %s",string);

   return p3*distr;
}