void DeltaT_sys(void)
{
// plot measurements of shielding in a magnetic field.

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

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

   char string[256];
   //
   // define canvas 

   TCanvas *c0 = new TCanvas("c0","Shower Profile",200,10,700,500);
   c0->SetBorderMode(0);
   c0->SetFillColor(0);
   c0->SetGridx();
   c0->SetGridy();
    
   // gStyle->SetPadRightMargin(0.35);
   // gStyle->SetPadLeftMargin(0.35);

   Double_t xmin=0;
   Double_t xmax = 20;
   Double_t showermin=0.0;
   Double_t showermax=20;
   Double_t ymin = 0;
   Double_t ymax=0.2;
   Double_t Eg=1;
   Double_t Ec=0.01;
   Double_t Cg=0.5;
   Double_t bb=0.5;
   Double_t aa=1+bb*(log(Eg/Ec)+Cg);
   printf ("Energy=%f, Shower parameter a=%f, b=%f\n",Eg,aa,bb);

   TF1 *shower_profile = new TF1("shower_profile",shower_profile,xmin,xmax,2);

   shower_profile->SetLineColor(2);
   shower_profile->SetParameters(aa,bb);
   // gPad->SetLogy(1);

   shower_profile->SetTitle("");
   // shower_profile->GetXaxis()->SetRangeUser(xmin,xmax);
   shower_profile->GetYaxis()->SetRangeUser(ymin,ymax);
   shower_profile->GetXaxis()->SetTitleSize(0.07);
   shower_profile->GetXaxis()->SetLabelSize(0.05);
   shower_profile->GetYaxis()->SetTitleSize(.05);
   //shower_profile->GetYaxis()->SetTitleOffset(1.5);
   shower_profile->GetXaxis()->SetTitle("t");
   shower_profile->GetYaxis()->SetTitle("Profile");
   shower_profile->GetXaxis()->SetNdivisions(10);
   shower_profile->GetYaxis()->SetNdivisions(10);
   shower_profile->Draw();
   Double_t integ = shower_profile->Integral(0,20);
   sprintf (string," integral=%f\n",integ);
   printf ("Shower integral=%s \n",string);
   sprintf (string,"P(t) = #frac{b (bt)^{a-1} exp(-bt)}{#Gamma(a)}\n");
   t1 = new TLatex(0.2,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"a=%.2f, b=%.2f\n",aa,bb);
   t1 = new TLatex(0.2,0.70,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();


   #define nlayers 10;
   #define npts nlayers+1;

   // Create a binned version of the shower profile

   Double_t shower_bin[npts];
   Double_t sum = 0;

   // define, fill and plot histograms

   Float_t xhistmin = xmin;
   Float_t xhistmax = xmax;
   Int_t ndim = 6;
   sprintf (string,"ndim=%.0d\n",ndim);
   t1 = new TLatex(0.2,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();



   TH1F *profile_bin = new TH1F ("profile_bin","Profile Binned",ndim,xhistmin,xhistmax);


   for (int jj=0;jj<nlayers;jj++) {
     Double_t min = jj*xmax/nlayers;
     Double_t max = (jj+1)*xmax/nlayers;
     shower_bin[jj] = shower_profile->Integral(min,max)/integ;
     sum=sum+shower_bin[jj];
     //printf (" jj=%d, min=%f, max=%f, bin=%f, sum=%f\n",jj,min,max,shower_bin[jj],sum);
   }
     

   // use N as the mean total statistics for a given event
 
   Double_t N=1000;
   Double_t a=sqrt(N);
   Double_t b=0.;
   Double_t c=0;
   // Double_t b=0.4;
   // Double_t c=0.2;

   Double_t sigma[npts];
   Double_t sigma_th[npts];
   Double_t layer[npts];

   //
   xmin = 0;
   xmax = 10;
   ymin = 0.5;
   ymax = 1.5;

   // set up constant term so that it is b for Ntotal
   Double_t thresh = N;
   Double_t N_min=N/10;
   // Double_t N_min=0;  

   Double_t sigmab_N = sqrt(a*a/N + b*b);
   Double_t sigmac_N = sqrt(a*a/N + c*c*(N*N)/(N*N));

   for (int jj=1;jj<npts;jj++) {
   Double_t Sumw=0;
   Int_t events=0;
   Double_t Sumw_th=0;
      for (int j=0;j<jj;j++) {
        Double_t min = j*showermax/jj;
        Double_t max = (j+1)*showermax/jj;
        shower_bin[j] = shower_profile->Integral(min,max)/integ;
        Double_t N_layer = N*shower_bin[j];
        // printf ("jj=%d, j=%d, xmin=%f, xmax=%f, N_layer=%f\n",jj,j,min,max,N_layer);
        Double_t sigj2 = a*a/N_layer + b*b;
        Double_t sigj_th2 = a*a/N_layer + c*c*(N*N)/(N_layer*N_layer);
        if (N_layer > N_min) {
        
//           include in sum only if enough events

           Sumw = Sumw + 1/sigj2;
           Sumw_th = Sumw_th + 1/sigj_th2;
           events = events + N_layer;
           // printf(" Sumw=%f, Sumw_th=%f, events=%d\n",Sumw,Sumw_th,events);
       }

         // fill histogram

         if (jj == ndim) {
             Float_t cut = (integ/(max-min))*(N_min/N);
             Float_t x = (min+max)/2.;
             Float_t y = shower_bin[j]/(max-min);
             profile_bin->Fill(x,y);
             //printf ("ndim=%d, jj=%d, x=%f, y=%f cut=%f\n",ndim,jj,x,y,cut);
         }
             
      }
      layer[jj]=jj;
      if (sigmab_N>0 && sigmac_N>0 && Sumw > 0) {
         sigma[jj]=sqrt(1/Sumw)/sigmab_N;
         sigma_th[jj]=sqrt(1/Sumw_th)/sigmac_N;
         }
         else {
         sigma[jj] = 0.;
         sigma_th[jj] = 0.;
         }


      //printf (" Event total=%d\n\n",events);
   }

   TLine *line = new TLine(showermin,cut,showermax,cut);
   profile_bin->Draw("same");
   line->SetLineColor(2);
   line->Draw("same"); 

   for (jj=1;jj<npts;jj++){ 
      // printf ("jj=%d, sigmab_N=%f, sigmac_N=%f, layer=%f, sigma=%f, sigma_th=%f\n",jj,sigmab_N, sigmac_N, layer[jj], sigma[jj], sigma_th[jj]);
   }

   TCanvas *c1 = new TCanvas("c1","Timing systematics",200,10,700,700);

   c1->SetGridx();
   c1->SetGridy();
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   TGraph *gr = new TGraph (npts,layer,sigma);
   gr->GetXaxis()->SetRangeUser(xmin,xmax);
   gr->GetYaxis()->SetRangeUser(ymin,ymax);
   gr->GetXaxis()->SetTitleSize(0.05);
   gr->GetYaxis()->SetTitleSize(0.05);
   gr->GetXaxis()->SetTitle("Longitudinal Segmentation");
   gr->GetYaxis()->SetTitle("Timing resolution (relative)");
   gr->SetMarkerColor(2);
   gr->SetMarkerStyle(20);
   gr->SetTitle("");
   // gr->Fit("pol1","","",76,110);
   // TF1 *fit = gr->GetFunction("pol1");
   // Double_t p0 = fit->GetParameter(0); 
   // Double_t p1 = fit->GetParameter(1); 
   gr->Draw("AP");
   TGraph *gr2 = new TGraph (npts,layer,sigma_th);
   gr2->SetMarkerColor(4);
   gr2->SetMarkerStyle(21);
   gr2->Draw("same,P");

   // TF1 *calc = new TF1("calc","x-91",91,110);
   // calc->SetLineStyle(2);
   // calc->Draw("same");

   // sprintf (string,"#sigma = %.1f/sqrt(N)+ %.1f\n",a,b);
   sprintf (string,"#sigma = #frac{a}{#sqrt{N_{j}}} #oplus b #oplus c #frac{N}{N_{j}} \n");
   // sprintf (string,"#sigma = #frac{#sqrt{N}}{#sqrt{N_{j}}} #oplus b\n");
   t1 = new TLatex(0.2,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"N = %0.f, N_{min} = %0.f\n",N,N_min);
   t1 = new TLatex(0.2,0.70,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"#sigma = #frac{a}{#sqrt{N_{j}}} #oplus %.1f\n",b);
   // sprintf (string,"b = %.1f\n",b);
   t1 = new TLatex(0.6,0.40,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   // sprintf (string,"b = %.1f(#frac{N}{N_{j}})\n",b);
   sprintf (string,"#sigma = #frac{a}{#sqrt{N_{j}}} #oplus %.1f #frac{N}{N_{j}}\n",c);
   t1 = new TLatex(0.6,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c0->SaveAs("Shower_profile.eps");
   c0->SaveAs("Shower_profile.gif");
   c1->SaveAs("DeltaT_sys.eps");
   c1->SaveAs("DeltaT_sys.gif");
}
Double_t shower_profile (Double_t *x, Double_t *par) {
     Double_t a=par[0];
     Double_t b=par[1];
     Double_t t=x[0]; 

     char string[256];

     Double_t func=0;
    
     func = a + b*t;
     func = b * pow(b*t,a-1) * exp(-b*t)/TMath::Gamma(a);
     // func = b * (b*t)*(b*t) * exp(-b*t)/TMath::Gamma(a);
 
        // func = func *(TMath::Erf(erfarg1) - TMath::Erf(erfarg2)); 


        // sprintf (string,"a=%f, b=%f, func=%f\n",a,b,func);
        // printf ("string=%s",string);

     return func;
}