void tof_calibration_tree(void)
{
// 
// Develop a toy model for calibrating the tof counters with particles
//

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

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

   char string[256];
   char filename[80];
   Int_t j,jj;
   
   #define npts 3;   // number of points in dummy plot
   #define np 1001;  // number of steps for pulse parameterization
   #define nevents 1000;    // number of events to generate and process


   // define structure

   struct Events_t {

      // generated values
      Float_t dE_V;
      Float_t dE_H;
      Float_t t0_H;
      Float_t t0_V;
      Float_t x;
      Float_t y;
      Float_t tw_R;
      Float_t tw_L;
      Float_t tw_T;
      Float_t tw_B;
      Float_t t_offset;
      Float_t sigma_t;
      Float_t veff;
      Float_t lambda;
      Float_t Vmax;
      Float_t Vthr;

      // reconstructed values
      Float_t P_R;
      Float_t P_L;
      Float_t P_T;
      Float_t P_B;
      Float_t t_R;
      Float_t t_L;
      Float_t t_T;
      Float_t t_B;
      Float_t tcor_R;
      Float_t tcor_L;
      Float_t tcor_T;
      Float_t tcor_B;
      };
   Events_t Events;

   // create a tree 

   sprintf(filename,"tof_tree.root");
   TFile *tfile = new TFile(filename,"RECREATE");
   TTree *gentof = new TTree("gentof","generate pulses in tof counters");
   gentof->Branch("Events",&Events.dE_V,
           "dE_V/F:dE_H:t0_H:t0_V:x:y:tw_R:tw_L:tw_T:tw_B:t_offset:sigma_t:veff:lambda:Vmax:Vthr:P_R:P_L:P_T:P_B:t_R:t_L:t_T:t_B:tcor_R:tcor_L:tcor_T:tcor_B");

   // prepare pulse parameters

    Double_t pulse_t[np];      // time axis
    Double_t pulse_V[np];    // voltage scale

   Double_t ptmin = -20;    // begin range of pulse
   Double_t ptmax= 120;   // end range of pulse
   TF1 *pulse = new TF1("pulse_func",pulse_func,ptmin,ptmax,4);
   // Double_t mu=1/8.;    // unis of 1/x
   Double_t mu=1/20.;    // unis of 1/x
   // Double_t sigma=6;     // units of x  
   Double_t sigma=2;     // units of x  
   Double_t norm=20;  // normalization 
   Double_t npe=100;    // dummy for now
   pulse->SetParameters(mu,sigma,norm,npe);  
  
   Events.Vmax = 2;    // max voltage on ADC

   Double_t t_max = -1;
   Double_t maxpulse = -1;
   Double_t maxj = -1;
   for (j=0;j<np;j++) {
         Double_t npt = ptmin + (ptmax-ptmin)*j/(np-1);
         pulse_t[j] = npt;
         pulse_V[j] = pulse->Eval(npt);
         if (pulse->Eval(npt) > maxpulse) {
             maxj = j;
             maxpulse = pulse->Eval(npt);
             t_max = npt;
             }
         // printf ("j=%d t=%f V=%f\n",j,pulse_t[j],pulse_V[j]);
         }
         printf ("maxj=%d t_max=%f maxpulse=%f\n",maxj,t_max,maxpulse);

   // normalize pulse
   norm = Events.Vmax*norm/maxpulse;   // set maximum to 2 V for 20 MeV energy deposition.
   pulse->SetParameters(mu,sigma,norm,npe); 

   // find time-walk from peak to threshold

   Double_t t_thresh = -100;
   Double_t t_walk = -1;
   Double_t V_thresh = 0.020;    // take 20 mV threshold
   for (j=0;j<maxj+1;j++) {
         Double_t npt = ptmin + (ptmax-ptmin)*j/(np-1);
         if (pulse->Eval(npt) > V_thresh && t_thresh== -100) {
             t_thresh = npt;
             }
         // printf ("j=%d t=%f V=%f\n",j,pulse_t[j],pulse_V[j]);
         }
     t_walk = t_max - t_thresh;     // t(corrected) = t(measured) + t_walk

     // printf ("t_thresh=%f, t_max=%f, t_walk=%f\n",t_thresh,t_max,t_walk); 

   // input data set first point (actually -100) to 100 so that it is not plotted.

   // use offset on x-axis to make data visible and give negative y-values for reference (not plotted)

   Double_t dummyx[npts]={-150,0,150};
   Double_t dummyy[npts]={-150,-150,-150};

   //
   TCanvas *c1 = new TCanvas("c1","c1 tof calibration",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   // c1->Divide(2,2);
   //c1->SetGridx();
   //c1->SetGridy();
   // c1->SetLogy();
   // c1->cd(1);
   // c1_1->SetBorderMode(0);
   // c1_1->SetFillColor(0);

   // dummy to draw axes

   Double_t ymin=-120;
   Double_t ymax=120;
   Double_t xmin=-120;
   Double_t xmax=120;

   TGraph *dummy = new TGraph (npts,dummyx,dummyy);
   TLegend *leg = new TLegend(0.35,0.86,0.85,0.9);
   dummy->SetMarkerColor(1);
   dummy->SetMarkerStyle(21);
   t1 = new TLatex(0.15,0.91,"Hamamatsu: Min pulse of max channel");
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw(); 
   dummy->SetTitle("");
   dummy->GetXaxis()->SetRangeUser(xmin,xmax);
   dummy->GetYaxis()->SetRangeUser(ymin,ymax);
   dummy->GetXaxis()->SetTitleSize(0.04);
   dummy->GetYaxis()->SetTitleSize(0.04);
   dummy->GetYaxis()->SetTitleOffset(1.5);
   dummy->GetXaxis()->SetTitle("Time (ns)");
   dummy->GetYaxis()->SetTitle("Amplitude (mV)");
   dummy->Draw("Ap");
   dummy->GetXaxis()->SetNdivisions(505);

   // use energy deposited in cell with max energy.

   Events.lambda=300.;
   Double_t y0=0;   // 60 cm is quarter point
   Double_t x0=0;     // 0 cm is at center
   Double_t width=6;
   Double_t xwidth=240;
   // Double_t xwidth=6;
   Double_t ycmin  = y0 - width/2;
   Double_t ycmax  = y0 + width/2;
   Double_t xcmin= x0 - xwidth/2;
   Double_t xcmax= x0 + xwidth/2;;
   Double_t tmin=0;     // units are ns and cm
   Double_t tmax=0;
   Events.Vthr=0.020; // LE disc threshold
   Events.Vmax=2; // Max voltage at ADC
   Events.veff=16;   // 16 cm/ns conversion of time to distance
   Events.sigma_t=sqrt(2)*0.080;   // assume 80 ps sigma for counter average
   Events.t_offset=25;   // use time offset to move times of "data"
   Double_t emax = 20;   // take 20 MeV for max energy

   Double_t peakE=2.5*2;   // Most probable energy deposition (mep) 2.5 cm * 2 MeV/g/cm3, with rho=1 g/cm3
   Double_t sigma_E=1.;

   // generate "data"
   TRandom1 *r = new TRandom1();

   for (j=0;j<nevents;j++) {
        Events.x = xcmin + r->Uniform()*(xcmax-xcmin);
        Events.y  = ycmin + r->Uniform()*(ycmax-ycmin);
        Float_t t0 = tmin +  r->Uniform()*(tmax-tmin);
        Events.t0_H = t0;
        Events.t0_V = t0;

      // generate distribution of energy deposited

      Events.dE_V = r->Landau(peakE,sigma_E);   
      // Events.dE_V = 10; 
      Events.dE_H = r->Landau(peakE,sigma_E);   
      // Events.dE_H = 10;
        Events.P_R = Events.Vmax*(Events.dE_H/emax)*exp(Events.x/Events.lambda);    // normalize to 2 V for emax
        Events.P_L = Events.Vmax*(Events.dE_H/emax)*exp(-Events.x/Events.lambda);     
        Events.P_T = Events.Vmax*(Events.dE_V/emax)*exp(Events.y/Events.lambda);
        Events.P_B = Events.Vmax*(Events.dE_V/emax)*exp(-Events.y/Events.lambda);


      // normalize to energy deposition
      norm = emax*Events.P_R/maxpulse;   // set maximum to 2 V for 20 MeV energy deposition;
      // norm = 0.5*emax/maxpulse;   // set maximum to 2 V for 20 MeV energy deposition;
      pulse->SetParameters(mu,sigma,norm,npe); 

      // generate time-walk and LE times relative to "peak time"

         Double_t Vj;
         Double_t Vj1=0;
         Double_t t_thresh = -100;
         Double_t t_walk = -1;
         for (jj=0;jj<maxj+1;jj++) {
              Double_t npt = ptmin + (ptmax-ptmin)*jj/(np-1);
              Double_t Vj =  pulse->Eval(npt) ;
              if (Vj > Events.Vthr && t_thresh== -100) {
                  t_thresh = npt - ((Vj-Events.Vthr)/(Vj-Vj1))*(ptmax-ptmin)/(np-1);
                  jj = maxj;
                  }
              // printf ("jj=%d t=%f Vj=%f\n",jj,npt,Vj);
              Vj1 = Vj;
              }
         // if (t_thresh == -100) t_thresh = ptmin + (ptmax-ptmin)*maxj/(np-1);
         Events.tw_R= t_max - t_thresh;     // t(corrected) = t(measured) + t_walk;
         // printf ("dE=%f, emax=%f, t_thresh=%f, t_max=%f, tw_R=%f\n",Events.dE_H,emax,t_thresh,t_max,Events.tw_R); 

        // "measured times" including time-walk

        Events.t_R = r->Gaus(Events.t0_H - Events.x/Events.veff,Events.sigma_t) - Events.tw_R+ Events.t_offset;
        // ;
        Double_t b0=10.03;
        Double_t b1=0.73;
        Double_t b0=11.9;
        Double_t b1=0.22;
        Events.tcor_R = Events.t_R - b0/(1+b1*sqrt(Events.P_R));

      norm = emax*Events.P_L/maxpulse;   // set maximum to 2 V for 20 MeV energy deposition;
      // norm = 0.5*emax/maxpulse;   // set maximum to 2 V for 20 MeV energy deposition;
      pulse->SetParameters(mu,sigma,norm,npe); 

      // generate time-walk and LE times relative to "peak time"

         Double_t Vj;
         Double_t Vj1=0;
         Double_t t_thresh = -100;
         Double_t t_walk = -1;
         Double_t V_thresh = 0.020;    // take 20 mV threshold
         for (jj=0;jj<maxj+1;jj++) {
              Double_t npt = ptmin + (ptmax-ptmin)*jj/(np-1);
              Double_t Vj =  pulse->Eval(npt) ;
              if (Vj > Events.Vthr && t_thresh== -100) {
                  t_thresh = npt - ((Vj-Events.Vthr)/(Vj-Vj1))*(ptmax-ptmin)/(np-1);
                  jj = maxj;
                  }
              // printf ("jj=%d t=%f Vj=%f\n",jj,npt,Vj);
              Vj1 = Vj;
              }
         // if (t_thresh == -100) t_thresh = ptmin + (ptmax-ptmin)*maxj/(np-1);
         Events.tw_L = t_max - t_thresh;     // t(corrected) = t(measured) + t_walk;
         // printf ("dE=%f, emax=%f, t_thresh=%f, t_max=%f, tw_L=%f\n",Events.dE_H,emax,t_thresh,t_max,Events.tw_L); 


        Events.t_L = r->Gaus(Events.t0_H + Events.x/Events.veff,Events.sigma_t) - Events.tw_L + Events.t_offset;
        Events.tcor_L = Events.t_L + Events.tw_L;
        Double_t b0=10.06;
        Double_t b1=0.73;
        Double_t b0=11.9;
        Double_t b1=0.22;
        Events.tcor_L = Events.t_L - b0/(1+b1*sqrt(Events.P_L));

      norm = emax*Events.P_T/maxpulse;   // set maximum to 2 V for 20 MeV energy deposition;
      // norm = 0.5*emax/maxpulse;   // set maximum to 2 V for 20 MeV energy deposition;
      pulse->SetParameters(mu,sigma,norm,npe); 

      // generate time-walk and LE times relative to "peak time"

         Double_t Vj;
         Double_t Vj1=0;
         Double_t t_thresh = -100;
         Double_t t_walk = -1;
         for (jj=0;jj<maxj+1;jj++) {
              Double_t npt = ptmin + (ptmax-ptmin)*jj/(np-1);
              Double_t Vj =  pulse->Eval(npt) ;
              if (Vj > Events.Vthr && t_thresh== -100) {
                  t_thresh = npt - ((Vj-Events.Vthr)/(Vj-Vj1))*(ptmax-ptmin)/(np-1);
                  jj = maxj;
                  }
              // printf ("jj=%d t=%f Vj=%f\n",jj,npt,Vj);
              Vj1 = Vj;
              }
         // if (t_thresh == -100) t_thresh = ptmin + (ptmax-ptmin)*maxj/(np-1);
         Events.tw_T = t_max - t_thresh;     // t(corrected) = t(measured) + t_walk;
         // printf ("dE=%f, emax=%f, t_thresh=%f, t_max=%f, tw_T=%f\n",Events.dE_V,emax,t_thresh,t_max,Events.tw_T); 


        Events.t_T = r->Gaus(Events.t0_V - Events.y/Events.veff,Events.sigma_t) - Events.tw_T + Events.t_offset;
        Events.tcor_T = Events.t_T + Events.tw_T;
        Double_t b0=13.53;
        Double_t b1=5.25;
        Double_t b0=11.9;
        Double_t b1=0.22;
        Events.tcor_T = Events.t_T - b0/(1+b1*sqrt(Events.P_T));

      norm = emax*Events.P_B/maxpulse;   // set maximum to 2 V for 20 MeV energy deposition;
      // norm = 0.5*emax/maxpulse;   // set maximum to 2 V for 20 MeV energy deposition;
      pulse->SetParameters(mu,sigma,norm,npe); 

      // generate time-walk and LE times relative to "peak time"

         Double_t Vj;
         Double_t Vj1=0;
         Double_t t_thresh = -100;
         Double_t t_walk = -1;
         for (jj=0;jj<maxj+1;jj++) {
              Double_t npt = ptmin + (ptmax-ptmin)*jj/(np-1);
              Double_t Vj =  pulse->Eval(npt) ;
              if (Vj > Events.Vthr && t_thresh== -100) {
                  t_thresh = npt - ((Vj-Events.Vthr)/(Vj-Vj1))*(ptmax-ptmin)/(np-1);
                  jj = maxj;
                  }
              // printf ("jj=%d t=%f Vj=%f\n",jj,npt,Vj);
              Vj1 = Vj;
              }
         // if (t_thresh == -100) t_thresh = ptmin + (ptmax-ptmin)*maxj/(np-1);
         Events.tw_B  = t_max - t_thresh;     // t(corrected) = t(measured) + t_walk;
         // printf ("dE=%f, emax=%f, t_thresh=%f, t_max=%f, tw_B=%f\n\n",Events.dE_V,emax,t_thresh,t_max,Events.tw_B); 

        Events.t_B = r->Gaus(Events.t0_V + Events.y/Events.veff,Events.sigma_t) - Events.tw_B + Events.t_offset;
        Events.tcor_B = Events.t_B + Events.tw_B;
        Double_t b0=10.18;
        Double_t b1=0.83;
        Double_t b0=11.9;
        Double_t b1=0.22;
        Events.tcor_B = Events.t_B - b0/(1+b1*sqrt(Events.P_B));

       // Fill root tree
      
      gentof->Fill();

       // printf ("x=%f, y=%f, tR=%f, tL=%f, tT=%f, tB=%f, xH=%f, yV=%f, t_H=%f, t_V=%f\n",x,y,t0,t_R,t_L,t_T,t_B,x_H,y_V,t_H,t_V);
       // printf("\ndE_V=%f dE_H=%f, PR=%f, PL=%f, PT=%f, PB=%f\n",dE_H,dE_H,P_R,P_L,P_T,P_B);
       }
 
   gentof->Print();


    TGraph *curve = new TGraph(np,pulse_t,pulse_V);
   curve->SetMarkerColor(1);
   curve->SetMarkerStyle(21);
   curve->SetMarkerSize(0.3);
   curve->SetTitle("");
   curve->GetXaxis()->SetRangeUser(xmin,xmax);
   //curve->GetYaxis()->SetRangeUser(xmin,xmax);
   curve->GetXaxis()->SetTitle("Time (ns)");
   curve->GetYaxis()->SetTitle("Pulse Height (arb. units)");
   curve->GetYaxis()->SetTitleOffset(1.5);
   curve->Draw("Ap");


   sprintf(filename,"tof_calibration_tree.eps");
   c1->SaveAs(filename);
   sprintf(filename,"tof_calibration_tree.png");
   c1->SaveAs(filename);

   tfile->Write();


}
Double_t pulse_func (Double_t *x, Double_t *par) {

   Double_t timeoffset=0;
   Double_t mu=par[0];
   Double_t sigma=par[1];
   Double_t norm=par[2];
   Double_t npe=par[3];
   Double_t x1=x[0]-timeoffset;
   Double_t pi=3.14159;

   Int_t j;
   char string[256];

   Double_t arg = (mu*sigma*sigma-x1)/(sqrt(2)*sigma);
   Double_t amplitude = (mu/2)* exp(-mu*x1 +(mu*sigma)*(mu*sigma)/2) *  TMath::Erfc(arg);

   // convert to voltage scale (mV).
   
   amplitude = norm*amplitude ;

   // sprintf (string,"amplitude=%f mu=%f sigma=%f arg=%f, npe=%f\n",amplitude,mu,sigma,arg,npe);
   // printf ("string=%s",string);

   return amplitude;
}