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

#include <TRandom.h>

gROOT->Reset();

gStyle->SetPalette(1,0);
// gStyle->SetOptStat(kFALSE);;
gStyle->SetOptStat(kTRUE);
gStyle->SetOptStat(1111111);
TH1::SetDefaultSumw2(kTRUE);
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 histname[80];
   char filename[80];
   Double_t pi=3.14159;

   # define nevents 100000;

   //
   // define canvas 

   TCanvas *c0 = new TCanvas("c0","Shower Profile",200,10,700,500);
   c0->SetBorderMode(0);
   c0->SetFillColor(0);
   c0->SetGridx();
   c0->SetGridy();
   c0->Divide(2,1);

   c0->cd(1);  
   c0_1->SetGridx();
   c0_1->SetGridy();
   c0_1->SetBorderMode(0);
   c0_1->SetFillColor(0);
    
   // gStyle->SetPadRightMargin(0.35);
   // gStyle->SetPadLeftMargin(0.35);

   Double_t xmin=0;
   Double_t xmax = 40;
   Double_t showermin=0.0;
   Double_t showermax=40;
   Double_t Eshower = 0.06;
   Double_t theta=90;
   Double_t sinthe = sin(theta*pi/180.);
   Double_t ymin = 0;
   Double_t ymax=800;
   Double_t ztgt=65;
   Double_t Rcal=65;
   Double_t z_upstream=9;   // includes 8cm length of light guide
   Double_t z_downstream=415;  // includes 8 cm length of light guide

   // compute location of impact
   Double_t zimpact = ztgt + Rcal/tan(theta*pi/180.);
   Double_t pimpact = (zimpact-ztgt)/cos(theta*pi/180.);    // position along trajectory in cm
   printf ("Eshower=%f GeV, theta=%f deg , zimpact = %f cm, pimpact=%f\n",Eshower,theta,zimpact,pimpact);

   Double_t Eg=Eshower;
   Double_t Ec=0.01;
   Double_t Cg=0.5;
   Double_t bb=0.5;
   Double_t aa=1+bb*(log(Eg/Ec)+Cg);
    // take normalization roughly consistent with http://www.jlab.org/Hall-D/software/wiki/images/b/b2/Npe_Ngamma.pdf
   Double_t npe_per_GeV_per_side= 2200;
   printf ("Energy=%f, Shower parameter a=%f, b=%f, Npe/GeV/side=%f\n",Eg,aa,bb,npe_per_GeV_per_side);

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

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

   shower_profile->SetTitle("");
   // shower_profile->GetXaxis()->SetRangeUser(xmin,xmax);
   shower_profile->GetYaxis()->SetRangeUser(ymin,Eshower*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 (rad. leng.)");
   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,showermax);
   sprintf (string," %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();


   // radiation length (cm) in calorimeter from NIM article 
   Double_t X0=1.45;
   printf (" Rad Length=%f cm\n",X0);

   /* #define nlayers 10;
   #define npts nlayers+1;
   Double_t  height[nlayers]={2.06,2.06,2.06,2.06,2.06,2.06,2.46,2.46,2.46,2.46};*/

   /*#define nlayers 3;
   #define npts nlayers+1;
   Double_t  height[nlayers]={6.18,6.18,9.84};*/

   // 4-layer 1-2-3-4 configuration
   /*#define nlayers 4;
   #define npts nlayers+1;
   Double_t  height[nlayers]={2.06,4.12,6.18,9.84};*/

   // 4-layer 2-2-2-4 configuration
   #define nlayers 4;
   #define npts nlayers+1;
   Double_t  height[nlayers]={4.12,4.12,4.12,9.84};

   Double_t path[nlayers];    // distance in cm along particle trajectory within bin
   Double_t pathstep=0;
   for (int jj=0;jj<nlayers;jj++) {
        pathstep = pathstep + height[jj]/sinthe;
        path[jj]= pathstep;
        // printf ("jj=%d, path[jj]=%f\n",jj,path[jj]);
        }

   // 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 = nlayers;
   sprintf (string,"N_{pe}/GeV/side=%d, N_{pe}=%.1f\n",npe_per_GeV_per_side,shower_profile->Integral(0,showermax));
   t1 = new TLatex(0.2,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"E=%.2f GeV, #theta=%.0f Deg\n",Eshower,theta);
   t1 = new TLatex(0.2,0.60,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   Double_t min=0;
   Double_t max=0;
   for (int jj=0;jj<nlayers;jj++) {
     max += height[jj]/(X0*sinthe);
   }
   printf ("histogram boundaries: min=%f max=%f ndim=%d\n",min,max,ndim);

    // Book histograms

    Double_t rmin=0;
    Double_t rmax=10;
    Int_t nbins=100;
    TH1F *profile_bin = new TH1F ("profile_bin","Profile Binned",ndim,min, max);
    TH1F *profile_radial = new TH1F ("profile_radial","Radial profile",nbins,rmin, rmax);
    TH1F *shower_psi = new TH1F ("shower_psi","Shower angle psi",nbins,-20, 20);

    TRandom1 *r1 = new TRandom1();
    TRandom1 *r2 = new TRandom1();
    TRandom1 *rr = new TRandom1();
     Double_t tau1=0.2;
     Double_t tau2=1.5;
     Double_t radius=0;
     Double_t apar=1.124;
     Double_t bpar=-0.072;
     Double_t tt = 4;
     Double_t rtest=apar + bpar*tt;
     Double_t rhom=3.63;  // Moliere radii for the Bcal

    for (jj=0;jj<nevents;jj++) {

       // generate radial distribution from two exponential distribution
       Double_t r = rr->Uniform();
       if (r < rtest) {
           radius = r1->Exp(tau1)*rhom;
           }
       else {
           radius = r2->Exp(tau2)*rhom;
       }
       profile_radial->Fill(radius);
        }

    Double_t tmin=0;
    Double_t tmax=50;
    Int_t nbins=2000;
    TH1F *upstream_1 = new TH1F ("upstream_1","Upstream: Time in layer 1",nbins,tmin,tmax);
    TH1F *upstream_2 = new TH1F ("upstream_2","Upstream: Time in layer 2",nbins,tmin,tmax);
    TH1F *upstream_3 = new TH1F ("upstream_3","Upstream: Time in layer 3",nbins,tmin,tmax);
    TH1F *upstream_4 = new TH1F ("upstream_4","Upstream: Time in layer 4",nbins,tmin,tmax);
    TH1F *upstream_5 = new TH1F ("upstream_5","Upstream: Time in layer 5",nbins,tmin,tmax);
    TH1F *upstream_6 = new TH1F ("upstream_6","Upstream: Time in layer 6",nbins,tmin,tmax);
    TH1F *upstream_7 = new TH1F ("upstream_7","Upstream: Time in layer 7",nbins,tmin,tmax);
    TH1F *upstream_8 = new TH1F ("upstream_8","Upstream: Time in layer 8",nbins,tmin,tmax);
    TH1F *upstream_9 = new TH1F ("upstream_9","Upstream: Time in layer 9",nbins,tmin,tmax);
    TH1F *upstream_10 = new TH1F ("upstream_10","Upstream: Time in layer 10",nbins,tmin,tmax);

    TH1F *downstream_1 = new TH1F ("downstream_1","Downstream: Time in layer 1",nbins,tmin,tmax);
    TH1F *downstream_2 = new TH1F ("downstream_2","Downstream: Time in layer 2",nbins,tmin,tmax);
    TH1F *downstream_3 = new TH1F ("downstream_3","Downstream: Time in layer 3",nbins,tmin,tmax);
    TH1F *downstream_4 = new TH1F ("downstream_4","Downstream: Time in layer 4",nbins,tmin,tmax);
    TH1F *downstream_5 = new TH1F ("downstream_5","Downstream: Time in layer 5",nbins,tmin,tmax);
    TH1F *downstream_6 = new TH1F ("downstream_6","Downstream: Time in layer 6",nbins,tmin,tmax);
    TH1F *downstream_7 = new TH1F ("downstream_7","Downstream: Time in layer 7",nbins,tmin,tmax);
    TH1F *downstream_8 = new TH1F ("downstream_8","Downstream: Time in layer 8",nbins,tmin,tmax);
    TH1F *downstream_9 = new TH1F ("downstream_9","Downstream: Time in layer 9",nbins,tmin,tmax);
    TH1F *downstream_10 = new TH1F ("downstream_10","Downstream: Time in layer 10",nbins,tmin,tmax);

    tmin=-25;
    tmax=25;
    TH1F *diff_1 = new TH1F ("diff_1","Diff: Time in layer 1",nbins,tmin,tmax);
    TH1F *diff_2 = new TH1F ("diff_2","Diff: Time in layer 2",nbins,tmin,tmax);
    TH1F *diff_3 = new TH1F ("diff_3","Diff: Time in layer 3",nbins,tmin,tmax);
    TH1F *diff_4 = new TH1F ("diff_4","Diff: Time in layer 4",nbins,tmin,tmax);
    TH1F *diff_5 = new TH1F ("diff_5","Diff: Time in layer 5",nbins,tmin,tmax);
    TH1F *diff_6 = new TH1F ("diff_6","Diff: Time in layer 6",nbins,tmin,tmax);
    TH1F *diff_7 = new TH1F ("diff_7","Diff: Time in layer 7",nbins,tmin,tmax);
    TH1F *diff_8 = new TH1F ("diff_8","Diff: Time in layer 8",nbins,tmin,tmax);
    TH1F *diff_9 = new TH1F ("diff_9","Diff: Time in layer 9",nbins,tmin,tmax);
    TH1F *diff_10 = new TH1F ("diff_10","Diff: Time in layer 10",nbins,tmin,tmax);

    tmin=0;
    tmax=50;
    TH1F *sum_1 = new TH1F ("sum_1","Sum: Time in layer 1",nbins,tmin,tmax);
    TH1F *sum_2 = new TH1F ("sum_2","Sum: Time in layer 2",nbins,tmin,tmax);
    TH1F *sum_3 = new TH1F ("sum_3","Sum: Time in layer 3",nbins,tmin,tmax);
    TH1F *sum_4 = new TH1F ("sum_4","Sum: Time in layer 4",nbins,tmin,tmax);
    TH1F *sum_5 = new TH1F ("sum_5","Sum: Time in layer 5",nbins,tmin,tmax);
    TH1F *sum_6 = new TH1F ("sum_6","Sum: Time in layer 6",nbins,tmin,tmax);
    TH1F *sum_7 = new TH1F ("sum_7","Sum: Time in layer 7",nbins,tmin,tmax);
    TH1F *sum_8 = new TH1F ("sum_8","Sum: Time in layer 8",nbins,tmin,tmax);
    TH1F *sum_9 = new TH1F ("sum_9","Sum: Time in layer 9",nbins,tmin,tmax);
    TH1F *sum_10 = new TH1F ("sum_10","Sum: Time in layer 10",nbins,tmin,tmax);

   Double_t min=0;
   Double_t max=0;

   for (int jj=0;jj<nlayers;jj++) {
     max = min + height[jj]/(X0*sinthe);
     shower_bin[jj] = shower_profile->Integral(min,max)/(max-min);
     profile_bin -> Fill((min+max)/2,shower_bin[jj]);
     sum=sum+shower_bin[jj];
     printf (" jj=%d, min=%f, max=%f, bin=%f, sum=%f\n",jj,min,max,shower_bin[jj],sum);
     min = max;
   }

   profile_bin->Draw("histsame");

   c0->cd(2);  
   c0_2->SetGridx();
   c0_2->SetGridy();
   c0_2->SetBorderMode(0);
   c0_2->SetFillColor(0);
   c0_2->SetLogy();

   profile_radial->GetXaxis()->SetTitleSize(0.05);
   profile_radial->GetXaxis()->SetTitle("Transverse pos (cm)");
   profile_radial->Draw("hist");

   TCanvas *k0 = new TCanvas("k0","k0 geom_dt",200,10,700,700);
   k0->SetBorderMode(0);
   k0->SetFillColor(0);
   k0->Divide(2,3);

   /*k0->cd(1);  
   k0_1->SetGridx();
   k0_1->SetGridy();
   k0_1->SetBorderMode(0);
   k0_1->SetFillColor(0);
   k0_1->SetLogy();*/

   // shower_psi->GetXaxis()->SetRangeUser(xmin,xmax);
   // shower_psi->GetYaxis()->SetRangeUser(ymin,ymax);
   shower_psi->GetXaxis()->SetTitleSize(0.05);
   shower_psi->GetYaxis()->SetTitleSize(0.05);
   shower_psi->GetYaxis()->SetTitleOffset(1.5);
   shower_psi->GetXaxis()->SetTitle("Shower Angle #psi (deg)");
   // shower_psi->GetYaxis()->SetTitle("Npe");
   // shower_psi->GetXaxis()->SetNdivisions(5);
   shower_psi->SetMarkerColor(1);
   shower_psi->SetMarkerStyle(20);
   shower_psi->SetMarkerSize(0.5);
   shower_psi->Draw();


   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]);
   }

   // compute time distribution within each readout cell
   
  Double_t step=0.1*X0;   // step in cm
  Double_t pos=step;   // position along particle trajectory in cm
  Int_t Npe=0;
  Int_t Npetot=0;
  Double_t tpos=0;
  Double_t t_upstream=0;
  Double_t t_downstream=0;
  Double_t zpos=0;
  Double_t c=30;   // speed of light in cm/ns
  Double_t veff=16.2; // speed of light in fibers in cm/ns
  // Double_t c=veff;
  for (Int_t jj=0;jj<nlayers;jj++) {

         // compute distribution along trajectory with pos+ = step;
         while (pos <= path[jj]) {
               Npe = shower_profile->Eval(pos/X0)*step/X0 +0.6;  // use 0.6 so truncation matches sum
                Npetot += Npe;
               // printf ("jj=%d pos=%f, posX0=%f, Npe=%f\n",jj,pos,pos/X0,Npe);

               // make one entry per phoelectron

               for (Int_t kk=0; kk<Npe;kk++) {

                    Double_t rtest=apar + bpar*pos/X0;
                     // generate radial distribution from two exponential distribution
                    Double_t tempsin = 2;
                    while (tempsin > 1) {
                       Double_t r = rr->Uniform();
                       if (r < rtest) {
                           radius = r1->Exp(tau1)*rhom;
                           }
                        else {
                           radius = r2->Exp(tau2)*rhom;
                           }
                        tempsin=radius/pos;
                        }
                   Double_t psi=180*asin(tempsin)/pi;
                   if (rr->Uniform() > 0.5) psi = -psi;
                   // psi = 0;   // eliminate shower spread
                   shower_psi->Fill(psi);

                   // compute time along trajectory up to this point
                   zpos = zimpact + pos*cos((theta+psi)*pi/180.);
                   Double_t zadd = 0.5*pos*sin((theta+psi)*pi/180.);   // fudge factor 0.5 since particles are not traveling in straight lines.
                   tpos = (pimpact+pos)/c - pimpact/c + zadd/c;    // time relative to impact time
                   t_upstream = tpos + (zpos-z_upstream)/veff;
                   t_downstream = tpos + (z_downstream-zpos)/veff;

                  if (jj==2 && t_downstream < 21.5) {
                      // printf ("NEGATIVE zimpact=%f zpos=%f tpos=%f theta=%f psi=%f z_upstream=%f  t_upstream=%f t_downstream=%f\n",zimpact,zpos,tpos,theta,psi,z_upstream,t_upstream,t_downstream);
                      }
             
                   // fill histograms 
                  if (jj == 0) {upstream_1->Fill(t_upstream); downstream_1->Fill(t_downstream); }
                  else if (jj == 1) {upstream_2->Fill(t_upstream); downstream_2->Fill(t_downstream);}
                  else if (jj == 2) {upstream_3->Fill(t_upstream); downstream_3->Fill(t_downstream);}
                  else if (jj == 3)  {upstream_4->Fill(t_upstream); downstream_4->Fill(t_downstream);}
                  else if (jj == 4) {upstream_5->Fill(t_upstream); downstream_5->Fill(t_downstream);}
                  else if (jj == 5)  {upstream_6->Fill(t_upstream); downstream_6->Fill(t_downstream);}
                  else if (jj == 6) {upstream_7->Fill(t_upstream); downstream_7->Fill(t_downstream);}
                  else if (jj == 7) {upstream_8->Fill(t_upstream); downstream_8->Fill(t_downstream);}
                  else if (jj == 8) {upstream_9->Fill(t_upstream); downstream_9->Fill(t_downstream);}
                  else if (jj == 9) {upstream_10->Fill(t_upstream); downstream_10->Fill(t_downstream);}
                  else {
                     printf("*** switch illegal case jj=%d\n",jj);
                     } 

                  Double_t t_diff = 0.5*(t_upstream - t_downstream);
                  Double_t t_sum = 0.5*(t_upstream + t_downstream);
                  if (jj == 0)         {diff_1->Fill(t_diff); sum_1->Fill(t_sum); }
                  else if (jj == 1) {diff_2->Fill(t_diff); sum_2->Fill(t_sum);}
                  else if (jj == 2) {diff_3->Fill(t_diff); sum_3->Fill(t_sum);}
                  else if (jj == 3)  {diff_4->Fill(t_diff); sum_4->Fill(t_sum);}
                  else if (jj == 4) {diff_5->Fill(t_diff); sum_5->Fill(t_sum);}
                  else if (jj == 5)  {diff_6->Fill(t_diff); sum_6->Fill(t_sum);}
                  else if (jj == 6) {diff_7->Fill(t_diff); sum_7->Fill(t_sum);}
                  else if (jj == 7) {diff_8->Fill(t_diff); sum_8->Fill(t_sum);}
                  else if (jj == 8) {diff_9->Fill(t_diff); sum_9->Fill(t_sum);}
                  else if (jj == 9) {diff_10->Fill(t_diff); sum_10->Fill(t_sum);}
                  else {
                     printf("*** switch illegal case jj=%d\n",jj);
                     }

                   // printf ("psi=%f pos=%f zimpact=%f pimpact=%f zpos=%f tpos=%f tup=%f tdown=%f tdiff=%f tsum=%f\n",psi,pos,zimpact,pimpact,zpos,tpos,t_upstream,t_downstream,t_diff,t_sum);
                  }
               
               pos+ = step;
               }
         }
        printf ("Summed Npe=%d\n",Npetot);
               
        // Output histograms

   TCanvas *c1 = new TCanvas("c1","c1 geom_dt",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->Divide(2,2);

   c1->cd(1);  
   c1_1->SetGridx();
   c1_1->SetGridy();
   c1_1->SetBorderMode(0);
   c1_1->SetFillColor(0);

   // upstream_1->SetTitle("Upstream");
  Double_t deltax = 2;
   xmin = upstream_1->GetMean() - deltax;
   xmax = upstream_1->GetMean() + deltax;
   upstream_1->GetXaxis()->SetRangeUser(xmin,xmax);
   // upstream_1->GetYaxis()->SetRangeUser(ymin,ymax);
   upstream_1->GetXaxis()->SetTitleSize(0.05);
   upstream_1->GetYaxis()->SetTitleSize(0.05);
   upstream_1->GetYaxis()->SetTitleOffset(1.5);
   upstream_1->GetXaxis()->SetTitle("Time  (ns)");
   // upstream_1->GetYaxis()->SetTitle("Npe");
   upstream_1->GetXaxis()->SetNdivisions(5);
   upstream_1->SetMarkerColor(1);
   upstream_1->SetMarkerStyle(20);
   upstream_1->SetMarkerSize(0.5);
   upstream_1->Draw();

   // retrieving information from the histogram

   Double_t bin_thresh=5;
   Int_t jbins = upstream_1->GetNbinsX();
   Double_t xlo = upstream_1->GetBinLowEdge(1);
   Double_t width = upstream_1->GetBinWidth(1);
   printf ("jndim=%d, xlo=%f, width=%f\n",ndim,xlo,width);
   Double_t slope, Dslope, offset, Doffset,eps12;

   Double_t error_up_1;
   Double_t time_up_1;
   for(Int_t j=0;j<jbins;j++) {
      Double_t xbin = xlo + j*width;
      jj = (xbin-xlo)/width;

      Int_t nfitpts=3;
      Double_t datay[3], errory[3];
      Double_t datax[3];
      Double_t errorx[3]={0.025,0.025,0.025};
      Double_t binx1 = upstream_1->GetBinCenter(jj+1);
      datax[1] = 0;  // relative to binx1
      datay[1]= upstream_1->GetBinContent(jj+1) - bin_thresh;  // compute relative to (x-binx1,y-thresh)
      errory[1] = upstream_1->GetBinError(jj+1);
      // if (datay[1] > bin_thresh) {
      if (datay[1] > 0) {
         // calculate relative to middle bin
         datax[0] = upstream_1->GetBinCenter(jj) -binx1;
         datay[0] = upstream_1->GetBinContent(jj) - bin_thresh;
         errory[0] = upstream_1->GetBinError(jj);
         datax[2] = upstream_1->GetBinCenter(jj+2) -binx1;
         datay[2] = upstream_1->GetBinContent(jj+2) - bin_thresh;
         errory[2] = upstream_1->GetBinError(jj+2);

         // compute slope using data points before and after

          printf("datax %f %f %f\n",datax[0],datax[1],datax[2]);
          printf("datay %f %f %f\n",datay[0],datay[1],datay[2]);
          printf("errory %f %f %f\n",errory[0],errory[1],errory[2]);
          TGraphErrors *data = new TGraphErrors (nfitpts,datax,datay,errorx,errory);
          TFitResultPtr ptr = data->Fit("pol1","S","",-0.05,0.05);
          TMatrixDSym cov = ptr->GetCovarianceMatrix();       // access to covariance matrix
          cov.Print();
          Double_t determinant = cov.Determinant();
          // data->Draw("AP");

         if (determinant  != 0) {
              Double_t row1[2], row2[2];
              cov.ExtractRow(0,0,row1);
              cov.ExtractRow(1,0,row2);
              printf ("matrix=%f %f %f %f, determinant=%f\n",row1[0],row1[1],row2[0],row2[1],determinant);
              offset = ptr->Value(0);
              slope = ptr->Value(1);
              Doffset= sqrt(row1[0]);
              Dslope = sqrt(row2[1]);
              eps12 = row1[1];
              if (abs(slope) > 1.e-6) {
                 time_up_1 = -offset/slope + binx1;
                 error_up_1 = sqrt(Doffset*Doffset + (offset/slope)*(offset/slope)*Dslope*Dslope -(offset/slope)*eps12)/slope;
                 }
                 else {
                 time_up_1 = binx1;
                 error_up_1 = 0.025;
                 }
              printf (" slope =%f Dslope=%f offset=%f Doffset=%f binx1=%f time_up_1 =%f error_up_1 =%f\n",slope,Dslope,offset,Doffset,binx1,time_up_1,error_up_1);
              }
              else {
              printf ("*** Delta = 0\n",Delta);
              slope = 0;
              error_up_1 = 0;
              time_up_1 = 0;
              }
          j = jbins;
          }
   }

   sprintf(string,"Eshower=%.2f GeV, #theta=%.0f deg\n",Eshower,theta);
   t1 = new TLatex(0.18,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Time=%.3f, Delta x=%.3f ns\n",time_up_1,error_up_1);
   t1 = new TLatex(0.18,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Thresh=%.1f p.e.\n",bin_thresh);
   t1 = new TLatex(0.18,0.70,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c1->cd(2); 
   c1_2->SetGridx();
   c1_2->SetGridy();
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);
   // c1_2->SetLogx();
   xmin = upstream_2->GetMean() - deltax;
   xmax = upstream_2->GetMean() + deltax;
   upstream_2->GetXaxis()->SetRangeUser(xmin,xmax);
   upstream_2->GetXaxis()->SetTitleSize(0.05);
   upstream_2->GetYaxis()->SetTitleSize(0.05);
   upstream_2->GetYaxis()->SetTitleOffset(1.5);
   upstream_2->GetXaxis()->SetTitle("Time (ns)");
   upstream_2->GetXaxis()->SetNdivisions(5);
   upstream_2->SetMarkerStyle(20);
   upstream_2->SetMarkerSize(0.5);
   upstream_2->Draw();


   // retrieving information from the histogram

   Int_t jbins = upstream_2->GetNbinsX();
   Double_t xlo = upstream_2->GetBinLowEdge(1);
   Double_t width = upstream_2->GetBinWidth(1);
   printf ("jndim=%d, xlo=%f, width=%f\n",ndim,xlo,width);

   Double_t error_up_2;
   Double_t time_up_2;
   for(Int_t j=0;j<jbins;j++) {
      Double_t xbin = xlo + j*width;
      jj = (xbin-xlo)/width;

      Int_t nfitpts=3;
      Double_t datay[3], errory[3];
      Double_t datax[3];
      Double_t errorx[3]={0.025,0.025,0.025};
      Double_t binx1 = upstream_2->GetBinCenter(jj+1);
      datax[1] = 0;  // relative to binx1
      datay[1]= upstream_2->GetBinContent(jj+1) - bin_thresh;  // compute relative to (x-binx1,y-thresh)
      errory[1] = upstream_2->GetBinError(jj+1);
      // if (datay[1] > bin_thresh) {
      if (datay[1] > 0) {
         // calculate relative to middle bin
         datax[0] = upstream_2->GetBinCenter(jj) -binx1;
         datay[0] = upstream_2->GetBinContent(jj) - bin_thresh;
         errory[0] = upstream_2->GetBinError(jj);
         datax[2] = upstream_2->GetBinCenter(jj+2) -binx1;
         datay[2] = upstream_2->GetBinContent(jj+2) - bin_thresh;
         errory[2] = upstream_2->GetBinError(jj+2);

         // compute slope using data points before and after

          printf("datax %f %f %f\n",datax[0],datax[1],datax[2]);
          printf("datay %f %f %f\n",datay[0],datay[1],datay[2]);
          printf("errory %f %f %f\n",errory[0],errory[1],errory[2]);
          TGraphErrors *data = new TGraphErrors (nfitpts,datax,datay,errorx,errory);
          TFitResultPtr ptr = data->Fit("pol1","S","",-0.05,0.05);
          TMatrixDSym cov = ptr->GetCovarianceMatrix();       // access to covariance matrix
          cov.Print();
          Double_t determinant = cov.Determinant();
          // data->Draw("AP");

         if (determinant  != 0) {
              Double_t row1[2], row2[2];
              cov.ExtractRow(0,0,row1);
              cov.ExtractRow(1,0,row2);
              printf ("matrix=%f %f %f %f, determinant=%f\n",row1[0],row1[1],row2[0],row2[1],determinant);
              offset = ptr->Value(0);
              slope = ptr->Value(1);
              Doffset= sqrt(row1[0]);
              Dslope = sqrt(row2[1]);
              eps12 = row1[1];
              if (abs(slope) > 1.e-6) {
                 time_up_2 = -offset/slope + binx1;
                 error_up_2 = sqrt(Doffset*Doffset + (offset/slope)*(offset/slope)*Dslope*Dslope -(offset/slope)*eps12)/slope;
                 }
                 else {
                 time_up_2 = binx1;
                 error_up_2 = 0.025;
                 }
              printf (" slope =%f Dslope=%f offset=%f Doffset=%f binx1=%f time_up_2 =%f error_up_2 =%f\n",slope,Dslope,offset,Doffset,binx1,time_up_2,error_up_2);
              }
              else {
              printf ("*** Delta = 0\n",Delta);
              slope = 0;
              error_up_2 = 0;
              time_up_2 = 0;
              }
          j = jbins;
          }
   }

   sprintf(string,"Time=%.3f, Delta x=%.3f ns\n",time_up_2,error_up_2);
   t1 = new TLatex(0.18,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();


   c1->cd(3);
   c1_3->SetGridx();
   c1_3->SetGridy();
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);
   xmin = upstream_3->GetMean() - deltax;
   xmax = upstream_3->GetMean() + deltax;
   upstream_3->GetXaxis()->SetRangeUser(xmin,xmax);
   upstream_3->GetXaxis()->SetTitleSize(0.05);
   upstream_3->GetYaxis()->SetTitleSize(0.05);
   upstream_3->GetYaxis()->SetTitleOffset(1.5);
   upstream_3->GetXaxis()->SetTitle("Time (ns)");
   upstream_3->GetXaxis()->SetNdivisions(5);
   upstream_3->SetMarkerStyle(20);
   upstream_3->SetMarkerSize(0.5);
   upstream_3->Draw();

   // retrieving information from the histogram

   Int_t jbins = upstream_3->GetNbinsX();
   Double_t xlo = upstream_3->GetBinLowEdge(1);
   Double_t width = upstream_3->GetBinWidth(1);
   printf ("jndim=%d, xlo=%f, width=%f\n",ndim,xlo,width);

   Double_t error_up_3;
   Double_t time_up_3;
   for(Int_t j=0;j<jbins;j++) {
      Double_t xbin = xlo + j*width;
      jj = (xbin-xlo)/width;

      Int_t nfitpts=3;
      Double_t datay[3], errory[3];
      Double_t datax[3];
      Double_t errorx[3]={0.025,0.025,0.025};
      Double_t binx1 = upstream_3->GetBinCenter(jj+1);
      datax[1] = 0;  // relative to binx1
      datay[1]= upstream_3->GetBinContent(jj+1) - bin_thresh;  // compute relative to (x-binx1,y-thresh)
      errory[1] = upstream_3->GetBinError(jj+1);
      // if (datay[1] > bin_thresh) {
      if (datay[1] > 0) {
         // calculate relative to middle bin
         datax[0] = upstream_3->GetBinCenter(jj) -binx1;
         datay[0] = upstream_3->GetBinContent(jj) - bin_thresh;
         errory[0] = upstream_3->GetBinError(jj);
         datax[2] = upstream_3->GetBinCenter(jj+2) -binx1;
         datay[2] = upstream_3->GetBinContent(jj+2) - bin_thresh;
         errory[2] = upstream_3->GetBinError(jj+2);

         // compute slope using data points before and after

          printf("datax %f %f %f\n",datax[0],datax[1],datax[2]);
          printf("datay %f %f %f\n",datay[0],datay[1],datay[2]);
          printf("errory %f %f %f\n",errory[0],errory[1],errory[2]);
          TGraphErrors *data = new TGraphErrors (nfitpts,datax,datay,errorx,errory);
          TFitResultPtr ptr = data->Fit("pol1","S","",-0.05,0.05);
          TMatrixDSym cov = ptr->GetCovarianceMatrix();       // access to covariance matrix
          cov.Print();
          Double_t determinant = cov.Determinant();
          // data->Draw("AP");

         if (determinant  != 0) {
              Double_t row1[2], row2[2];
              cov.ExtractRow(0,0,row1);
              cov.ExtractRow(1,0,row2);
              printf ("matrix=%f %f %f %f, determinant=%f\n",row1[0],row1[1],row2[0],row2[1],determinant);
              offset = ptr->Value(0);
              slope = ptr->Value(1);
              Doffset= sqrt(row1[0]);
              Dslope = sqrt(row2[1]);
              eps12 = row1[1];
              if (abs(slope) > 1.e-6) {
                 time_up_3 = -offset/slope + binx1;
                 error_up_3 = sqrt(Doffset*Doffset + (offset/slope)*(offset/slope)*Dslope*Dslope -(offset/slope)*eps12)/slope;
                 }
                 else {
                 time_up_3 = binx1;
                 error_up_3 = 0.025;
                 }
              printf (" slope =%f Dslope=%f offset=%f Doffset=%f binx1=%f time_up_3 =%f error_up_3 =%f\n",slope,Dslope,offset,Doffset,binx1,time_up_3,error_up_3);
              }
              else {
              printf ("*** Delta = 0\n",Delta);
              slope = 0;
              error_up_3 = 0;
              time_up_3 = 0;
              }
          j = jbins;
          }
   }

   sprintf(string,"Time=%.3f, Delta x=%.3f ns\n",time_up_3,error_up_3);
   t1 = new TLatex(0.18,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();



   c1->cd(4);
   c1_4->SetGridx();
   c1_4->SetGridy();
   c1_4->SetBorderMode(0);
   c1_4->SetFillColor(0);
   xmin = upstream_4->GetMean() - deltax;
   xmax = upstream_4->GetMean() + deltax;
   upstream_4->GetXaxis()->SetRangeUser(xmin,xmax);
   upstream_4->GetXaxis()->SetTitleSize(0.05);
   upstream_4->GetYaxis()->SetTitleSize(0.05);
   upstream_4->GetYaxis()->SetTitleOffset(1.5);
   upstream_4->GetXaxis()->SetTitle("Time (ns)");
   upstream_4->GetXaxis()->SetNdivisions(5);
   upstream_4->SetMarkerStyle(20);
   upstream_4->SetMarkerSize(0.5);
   upstream_4->Draw();


   // retrieving information from the histogram

   Int_t jbins = upstream_4->GetNbinsX();
   Double_t xlo = upstream_4->GetBinLowEdge(1);
   Double_t width = upstream_4->GetBinWidth(1);
   printf ("jndim=%d, xlo=%f, width=%f\n",ndim,xlo,width);



   Double_t error_up_4;
   Double_t time_up_4;
   for(Int_t j=0;j<jbins;j++) {
      Double_t xbin = xlo + j*width;
      jj = (xbin-xlo)/width;

      Int_t nfitpts=3;
      Double_t datay[3], errory[3];
      Double_t datax[3];
      Double_t errorx[3]={0.025,0.025,0.025};
      Double_t binx1 = upstream_4->GetBinCenter(jj+1);
      datax[1] = 0;  // relative to binx1
      datay[1]= upstream_4->GetBinContent(jj+1) - bin_thresh;  // compute relative to (x-binx1,y-thresh)
      errory[1] = upstream_4->GetBinError(jj+1);
      // if (datay[1] > bin_thresh) {
      if (datay[1] > 0) {
         // calculate relative to middle bin
         datax[0] = upstream_4->GetBinCenter(jj) -binx1;
         datay[0] = upstream_4->GetBinContent(jj) - bin_thresh;
         errory[0] = upstream_4->GetBinError(jj);
         datax[2] = upstream_4->GetBinCenter(jj+2) -binx1;
         datay[2] = upstream_4->GetBinContent(jj+2) - bin_thresh;
         errory[2] = upstream_4->GetBinError(jj+2);

         // compute slope using data points before and after

          printf("datax %f %f %f\n",datax[0],datax[1],datax[2]);
          printf("datay %f %f %f\n",datay[0],datay[1],datay[2]);
          printf("errory %f %f %f\n",errory[0],errory[1],errory[2]);
          TGraphErrors *data = new TGraphErrors (nfitpts,datax,datay,errorx,errory);
          TFitResultPtr ptr = data->Fit("pol1","S","",-0.05,0.05);
          TMatrixDSym cov = ptr->GetCovarianceMatrix();       // access to covariance matrix
          cov.Print();
          Double_t determinant = cov.Determinant();
          printf (" Determinant=%f\n",determinant);
          // data->Draw("AP");

         if (determinant  != 0) {
              Double_t row1[2], row2[2];
              cov.ExtractRow(0,0,row1);
              cov.ExtractRow(1,0,row2);
              printf ("matrix=%f %f %f %f, determinant=%f\n",row1[0],row1[1],row2[0],row2[1],determinant);
              offset = ptr->Value(0);
              slope = ptr->Value(1);
              Doffset= sqrt(row1[0]);
              Dslope = sqrt(row2[1]);
              eps12 = row1[1];
              if (abs(slope) > 1.e-6) {
                 time_up_4 = -offset/slope + binx1;
                 error_up_4 = sqrt(Doffset*Doffset + (offset/slope)*(offset/slope)*Dslope*Dslope -(offset/slope)*eps12)/slope;
                 }
                 else {
                 time_up_4 = binx1;
                 error_up_4  = 0.025;
                 }
              printf (" slope =%f Dslope=%f offset=%f Doffset=%f binx1=%f time_up_4 =%f error_up_4 =%f\n",slope,Dslope,offset,Doffset,binx1,time_up_4,error_up_4);
              }
              else {
              printf ("*** Delta = 0\n",Delta);
              slope = 0;
              error_up_4 = 0;
              time_up_4 = 0;
              }
          j = jbins;
          }
   }


   sprintf(string,"Time=%.3f, Delta x=%.3f ns\n",time_up_4,error_up_4);
   t1 = new TLatex(0.18,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();



   /*c1->cd(5);
   c1_5->SetGridx();
   c1_5->SetGridy();
   c1_5->SetBorderMode(0);
   c1_5->SetFillColor(0);
   xmin = upstream_5->GetMean() - deltax;
   xmax = upstream_5->GetMean() + deltax;
   upstream_5->GetXaxis()->SetRangeUser(xmin,xmax);
   upstream_5->GetXaxis()->SetTitleSize(0.05);
   upstream_5->GetYaxis()->SetTitleSize(0.05);
   upstream_5->GetYaxis()->SetTitleOffset(1.5);
   upstream_5->GetXaxis()->SetTitle("Time (ns)");
   upstream_5->GetXaxis()->SetNdivisions(5);
   upstream_5->SetMarkerStyle(20);
   upstream_5->SetMarkerSize(0.5);
   upstream_5->Draw();

   c1->cd(6);
   c1_6->SetGridx();
   c1_6->SetGridy();
   c1_6->SetBorderMode(0);
   c1_6->SetFillColor(0);
   xmin = upstream_6->GetMean() - deltax;
   xmax = upstream_6->GetMean() + deltax;
   upstream_6->GetXaxis()->SetRangeUser(xmin,xmax);
   upstream_6->GetXaxis()->SetTitleSize(0.05);
   upstream_6->GetYaxis()->SetTitleSize(0.05);
   upstream_6->GetYaxis()->SetTitleOffset(1.5);
   upstream_6->GetXaxis()->SetTitle("Time (ns)");
   upstream_6->GetXaxis()->SetNdivisions(5);
   upstream_6->SetMarkerStyle(20);
   upstream_6->SetMarkerSize(0.5);
   upstream_6->Draw();*/

   /*TCanvas *c2 = new TCanvas("c2","c2 geom_dt",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   c2->Divide(2,3);

   c2->cd(1);  
   c2_1->SetGridx();
   c2_1->SetGridy();
   c2_1->SetBorderMode(0);
   c2_1->SetFillColor(0);


   xmin = upstream_7->GetMean() - deltax;
   xmax = upstream_7->GetMean() + deltax;
   upstream_7->GetXaxis()->SetRangeUser(xmin,xmax);
   upstream_7->GetXaxis()->SetTitleSize(0.05);
   upstream_7->GetYaxis()->SetTitleSize(0.05);
   upstream_7->GetYaxis()->SetTitleOffset(1.5);
   upstream_7->GetXaxis()->SetTitle("Time (ns)");
   upstream_7->GetXaxis()->SetNdivisions(5);
   upstream_7->SetMarkerStyle(20);
   upstream_7->SetMarkerSize(0.5);
   upstream_7->SetMarkerColor(1);
   upstream_7->Draw();

   sprintf(string,"Eshower=%.2f GeV, #theta=%.0f deg\n",Eshower,theta);
   t1 = new TLatex(0.18,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c2->cd(2); 
   c2_2->SetGridx();
   c2_2->SetGridy();
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);
   xmin = upstream_8->GetMean() - deltax;
   xmax = upstream_8->GetMean() + deltax;
   upstream_8->GetXaxis()->SetRangeUser(xmin,xmax);
   upstream_8->GetXaxis()->SetTitleSize(0.05);
   upstream_8->GetYaxis()->SetTitleSize(0.05);
   upstream_8->GetYaxis()->SetTitleOffset(1.5);
   upstream_8->GetXaxis()->SetTitle("Time (ns)");
   upstream_8->GetXaxis()->SetNdivisions(5);
   upstream_8->SetMarkerStyle(20);
   upstream_8->SetMarkerSize(0.5);
   upstream_8->Draw();

   c2->cd(3);
   c2_3->SetGridx();
   c2_3->SetGridy();
   c2_3->SetBorderMode(0);
   c2_3->SetFillColor(0);
   xmin = upstream_9->GetMean() - deltax;
   xmax = upstream_9->GetMean() + deltax;
   upstream_9->GetXaxis()->SetRangeUser(xmin,xmax);
   upstream_9->GetXaxis()->SetTitleSize(0.05);
   upstream_9->GetYaxis()->SetTitleSize(0.05);
   upstream_9->GetYaxis()->SetTitleOffset(1.5);
   upstream_9->GetXaxis()->SetTitle("Time (ns)");
   upstream_9->GetXaxis()->SetNdivisions(5);
   upstream_9->SetMarkerStyle(20);
   upstream_9->SetMarkerSize(0.5);
   upstream_9->Draw();

   c2->cd(4);
   c2_4->SetGridx();
   c2_4->SetGridy();
   c2_4->SetBorderMode(0);
   c2_4->SetFillColor(0);
   xmin = upstream_10->GetMean() - deltax;
   xmax = upstream_10->GetMean() + deltax;
   upstream_10->GetXaxis()->SetRangeUser(xmin,xmax);
   upstream_10->GetXaxis()->SetTitleSize(0.05);
   upstream_10->GetYaxis()->SetTitleSize(0.05);
   upstream_10->GetYaxis()->SetTitleOffset(1.5);
   upstream_10->GetXaxis()->SetTitle("Time (ns)");
   upstream_10->GetXaxis()->SetNdivisions(5);
   upstream_10->SetMarkerStyle(20);
   upstream_10->SetMarkerSize(0.5);
   upstream_10->Draw();*/


   TCanvas *c3 = new TCanvas("c3","c3 geom_dt",200,10,700,700);
   c3->SetBorderMode(0);
   c3->SetFillColor(0);
   c3->Divide(2,2);

   c3->cd(1);  
   c3_1->SetGridx();
   c3_1->SetGridy();
   c3_1->SetBorderMode(0);
   c3_1->SetFillColor(0);

   // downstream_1->SetTitle("Downstream");
   xmin = downstream_1->GetMean() - deltax;
   xmax = downstream_1->GetMean() + deltax;
   downstream_1->GetXaxis()->SetRangeUser(xmin,xmax);
   // downstream_1->GetYaxis()->SetRangeUser(ymin,ymax);
   downstream_1->GetXaxis()->SetTitleSize(0.05);
   downstream_1->GetYaxis()->SetTitleSize(0.05);
   downstream_1->GetYaxis()->SetTitleOffset(1.5);
   downstream_1->GetXaxis()->SetTitle("Time (ns)");
   // downstream_1->GetYaxis()->SetTitle("Npe");
   downstream_1->GetXaxis()->SetNdivisions(5);
   downstream_1->SetMarkerStyle(20);
   downstream_1->SetMarkerSize(0.5);
   downstream_1->SetMarkerColor(1);
   downstream_1->Draw();

   // retrieving information from the histogram

   Int_t jbins = downstream_1->GetNbinsX();
   Double_t xlo = downstream_1->GetBinLowEdge(1);
   Double_t width = downstream_1->GetBinWidth(1);
   printf ("jndim=%d, xlo=%f, width=%f\n",ndim,xlo,width);

   Double_t error_dn_1;
   Double_t time_dn_1;
   for(Int_t j=0;j<jbins;j++) {
      Double_t xbin = xlo + j*width;
      jj = (xbin-xlo)/width;

      Int_t nfitpts=3;
      Double_t datay[3], errory[3];
      Double_t datax[3];
      Double_t errorx[3]={0.025,0.025,0.025};
      Double_t binx1 = downstream_1->GetBinCenter(jj+1);
      datax[1] = 0;  // relative to binx1
      datay[1]= downstream_1->GetBinContent(jj+1) - bin_thresh;  // compute relative to (x-binx1,y-thresh)
      errory[1] = downstream_1->GetBinError(jj+1);
      // if (datay[1] > bin_thresh) {
      if (datay[1] > 0) {
         // calculate relative to middle bin
         datax[0] = downstream_1->GetBinCenter(jj) -binx1;
         datay[0] = downstream_1->GetBinContent(jj) - bin_thresh;
         errory[0] = downstream_1->GetBinError(jj);
         datax[2] = downstream_1->GetBinCenter(jj+2) -binx1;
         datay[2] = downstream_1->GetBinContent(jj+2) - bin_thresh;
         errory[2] = downstream_1->GetBinError(jj+2);

         // compute slope using data points before and after

          printf("datax %f %f %f\n",datax[0],datax[1],datax[2]);
          printf("datay %f %f %f\n",datay[0],datay[1],datay[2]);
          printf("errory %f %f %f\n",errory[0],errory[1],errory[2]);
          TGraphErrors *data = new TGraphErrors (nfitpts,datax,datay,errorx,errory);
          TFitResultPtr ptr = data->Fit("pol1","S","",-0.05,0.05);
          TMatrixDSym cov = ptr->GetCovarianceMatrix();       // access to covariance matrix
          cov.Print();
          Double_t determinant = cov.Determinant();
          // data->Draw("AP");

         if (determinant  != 0) {
              Double_t row1[2], row2[2];
              cov.ExtractRow(0,0,row1);
              cov.ExtractRow(1,0,row2);
              printf ("matrix=%f %f %f %f, determinant=%f\n",row1[0],row1[1],row2[0],row2[1],determinant);
              offset = ptr->Value(0);
              slope = ptr->Value(1);
              Doffset= sqrt(row1[0]);
              Dslope = sqrt(row2[1]);
              eps12 = row1[1];
              if (abs(slope) > 1.e-6) {
                 time_dn_1 = -offset/slope + binx1;
                 error_dn_1 = sqrt(Doffset*Doffset + (offset/slope)*(offset/slope)*Dslope*Dslope -(offset/slope)*eps12)/slope;
                 }
                 else {
                 time_dn_1 = binx1;
                 error_dn_1  = 0.025;
                 }
              printf (" slope =%f Dslope=%f offset=%f Doffset=%f binx1=%f time_dn_1 =%f error_dn_1 =%f\n",slope,Dslope,offset,Doffset,binx1,time_dn_1,error_dn_1);
              }
              else {
              printf ("*** Delta = 0\n",Delta);
              slope = 0;
              error_dn_1 = 0;
              time_dn_1 = 0;
              }
          j = jbins;
          }
   }

   sprintf(string,"Time=%.3f, Delta x=%.3f ns\n",time_dn_1,error_dn_1);
   t1 = new TLatex(0.18,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Thresh=%.1f p.e.\n",bin_thresh);
   t1 = new TLatex(0.18,0.70,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();


   sprintf(string,"Eshower=%.2f GeV, #theta=%.0f deg\n",Eshower,theta);
   t1 = new TLatex(0.18,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c3->cd(2); 
   c3_2->SetGridx();
   c3_2->SetGridy();
   c3_2->SetBorderMode(0);
   c3_2->SetFillColor(0);
   // c3_2->SetLogx();
   xmin = downstream_2->GetMean() - deltax;
   xmax = downstream_2->GetMean() + deltax;
   downstream_2->GetXaxis()->SetRangeUser(xmin,xmax);
   downstream_2->GetXaxis()->SetTitleSize(0.05);
   downstream_2->GetYaxis()->SetTitleSize(0.05);
   downstream_2->GetYaxis()->SetTitleOffset(1.5);
   downstream_2->GetXaxis()->SetTitle("Time (ns)");
   downstream_2->GetXaxis()->SetNdivisions(5);
   downstream_2->SetMarkerStyle(20);
   downstream_2->SetMarkerSize(0.5);
   downstream_2->Draw();

   // retrieving information from the histogram

   Int_t jbins = downstream_2->GetNbinsX();
   Double_t xlo = downstream_2->GetBinLowEdge(1);
   Double_t width = downstream_2->GetBinWidth(1);
   printf ("jndim=%d, xlo=%f, width=%f\n",ndim,xlo,width);

   Double_t error_dn_2;
   Double_t time_dn_2;
   for(Int_t j=0;j<jbins;j++) {
      Double_t xbin = xlo + j*width;
      jj = (xbin-xlo)/width;

      Int_t nfitpts=3;
      Double_t datay[3], errory[3];
      Double_t datax[3];
      Double_t errorx[3]={0.025,0.025,0.025};
      Double_t binx1 = downstream_2->GetBinCenter(jj+1);
      datax[1] = 0;  // relative to binx1
      datay[1]= downstream_2->GetBinContent(jj+1) - bin_thresh;  // compute relative to (x-binx1,y-thresh)
      errory[1] = downstream_2->GetBinError(jj+1);
      // if (datay[1] > bin_thresh) {
      if (datay[1] > 0) {
         // calculate relative to middle bin
         datax[0] = downstream_2->GetBinCenter(jj) -binx1;
         datay[0] = downstream_2->GetBinContent(jj) - bin_thresh;
         errory[0] = downstream_2->GetBinError(jj);
         datax[2] = downstream_2->GetBinCenter(jj+2) -binx1;
         datay[2] = downstream_2->GetBinContent(jj+2) - bin_thresh;
         errory[2] = downstream_2->GetBinError(jj+2);

         // compute slope using data points before and after

          printf("datax %f %f %f\n",datax[0],datax[1],datax[2]);
          printf("datay %f %f %f\n",datay[0],datay[1],datay[2]);
          printf("errory %f %f %f\n",errory[0],errory[1],errory[2]);
          TGraphErrors *data = new TGraphErrors (nfitpts,datax,datay,errorx,errory);
          TFitResultPtr ptr = data->Fit("pol1","S","",-0.05,0.05);
          TMatrixDSym cov = ptr->GetCovarianceMatrix();       // access to covariance matrix
          cov.Print();
          Double_t determinant = cov.Determinant();
          // data->Draw("AP");

         if (determinant  != 0) {
              Double_t row1[2], row2[2];
              cov.ExtractRow(0,0,row1);
              cov.ExtractRow(1,0,row2);
              printf ("matrix=%f %f %f %f, determinant=%f\n",row1[0],row1[1],row2[0],row2[1],determinant);
              offset = ptr->Value(0);
              slope = ptr->Value(1);
              Doffset= sqrt(row1[0]);
              Dslope = sqrt(row2[1]);
              eps12 = row1[1];
              if (abs(slope) > 1.e-6) {
                 time_dn_2 = -offset/slope + binx1;
                 error_dn_2 = sqrt(Doffset*Doffset + (offset/slope)*(offset/slope)*Dslope*Dslope -(offset/slope)*eps12)/slope;
                 }
                 else {
                 time_dn_2 = binx1;
                 error_dn_2  = 0.025;
                 }
              printf (" slope =%f Dslope=%f offset=%f Doffset=%f binx1=%f time_dn_2 =%f error_dn_2 =%f\n",slope,Dslope,offset,Doffset,binx1,time_dn_2,error_dn_2);
              }
              else {
              printf ("*** Delta = 0\n",Delta);
              slope = 0;
              error_dn_2 = 0;
              time_dn_2 = 0;
              }
          j = jbins;
          }
   }



   sprintf(string,"Time=%.3f, Delta x=%.3f ns\n",time_dn_2,error_dn_2);
   t1 = new TLatex(0.18,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();


   c3->cd(3);
   c3_3->SetGridx();
   c3_3->SetGridy();
   c3_3->SetBorderMode(0);
   c3_3->SetFillColor(0); 
   xmin = downstream_3->GetMean() - deltax;
   xmax = downstream_3->GetMean() + deltax;
   downstream_3->GetXaxis()->SetRangeUser(xmin,xmax);
   downstream_3->GetXaxis()->SetTitleSize(0.05);
   downstream_3->GetYaxis()->SetTitleSize(0.05);
   downstream_3->GetYaxis()->SetTitleOffset(1.5);
   downstream_3->GetXaxis()->SetTitle("Time (ns)");
   downstream_3->GetXaxis()->SetNdivisions(5);
   downstream_3->SetMarkerStyle(20);
   downstream_3->SetMarkerSize(0.5);
   downstream_3->Draw();


   // retrieving information from the histogram

   Int_t jbins = downstream_3->GetNbinsX();
   Double_t xlo = downstream_3->GetBinLowEdge(1);
   Double_t width = downstream_3->GetBinWidth(1);
   printf ("jndim=%d, xlo=%f, width=%f\n",ndim,xlo,width);

   Double_t error_dn_3;
   Double_t time_dn_3;
   for(Int_t j=0;j<jbins;j++) {
      Double_t xbin = xlo + j*width;
      jj = (xbin-xlo)/width;

      Int_t nfitpts=3;
      Double_t datay[3], errory[3];
      Double_t datax[3];
      Double_t errorx[3]={0.025,0.025,0.025};
      Double_t binx1 = downstream_3->GetBinCenter(jj+1);
      datax[1] = 0;  // relative to binx1
      datay[1]= downstream_3->GetBinContent(jj+1) - bin_thresh;  // compute relative to (x-binx1,y-thresh)
      errory[1] = downstream_3->GetBinError(jj+1);
      // if (datay[1] > bin_thresh) {
      if (datay[1] > 0) {
         // calculate relative to middle bin
         datax[0] = downstream_3->GetBinCenter(jj) -binx1;
         datay[0] = downstream_3->GetBinContent(jj) - bin_thresh;
         errory[0] = downstream_3->GetBinError(jj);
         datax[2] = downstream_3->GetBinCenter(jj+2) -binx1;
         datay[2] = downstream_3->GetBinContent(jj+2) - bin_thresh;
         errory[2] = downstream_3->GetBinError(jj+2);

         // compute slope using data points before and after

          printf("datax %f %f %f\n",datax[0],datax[1],datax[2]);
          printf("datay %f %f %f\n",datay[0],datay[1],datay[2]);
          printf("errory %f %f %f\n",errory[0],errory[1],errory[2]);
          TGraphErrors *data = new TGraphErrors (nfitpts,datax,datay,errorx,errory);
          TFitResultPtr ptr = data->Fit("pol1","S","",-0.05,0.05);
          TMatrixDSym cov = ptr->GetCovarianceMatrix();       // access to covariance matrix
          cov.Print();
          Double_t determinant = cov.Determinant();
          // data->Draw("AP");

         if (determinant  != 0) {
              Double_t row1[2], row2[2];
              cov.ExtractRow(0,0,row1);
              cov.ExtractRow(1,0,row2);
              printf ("matrix=%f %f %f %f, determinant=%f\n",row1[0],row1[1],row2[0],row2[1],determinant);
              offset = ptr->Value(0);
              slope = ptr->Value(1);
              Doffset= sqrt(row1[0]);
              Dslope = sqrt(row2[1]);
              eps12 = row1[1];
              if (abs(slope) > 1.e-6) {
                 time_dn_3 = -offset/slope + binx1;
                 error_dn_3 = sqrt(Doffset*Doffset + (offset/slope)*(offset/slope)*Dslope*Dslope -(offset/slope)*eps12)/slope;
                 }
                 else {
                 time_dn_3 = binx1;
                 error_dn_3  = 0.025;
                 }
              printf (" slope =%f Dslope=%f offset=%f Doffset=%f binx1=%f time_dn_3 =%f error_dn_3 =%f\n",slope,Dslope,offset,Doffset,binx1,time_dn_3,error_dn_3);
              }
              else {
              printf ("*** Delta = 0\n",Delta);
              slope = 0;
              error_dn_3 = 0;
              time_dn_3 = 0;
              }
          j = jbins;
          }
   }


   sprintf(string,"Time=%.3f, Delta x=%.3f ns\n",time_dn_3,error_dn_3);
   t1 = new TLatex(0.18,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();


   c3->cd(4);
   c3_4->SetGridx();
   c3_4->SetGridy();
   c3_4->SetBorderMode(0);
   c3_4->SetFillColor(0);
   xmin = downstream_4->GetMean() - deltax;
   xmax = downstream_4->GetMean() + deltax;
   downstream_4->GetXaxis()->SetRangeUser(xmin,xmax);
   downstream_4->GetXaxis()->SetTitleSize(0.05);
   downstream_4->GetYaxis()->SetTitleSize(0.05);
   downstream_4->GetYaxis()->SetTitleOffset(1.5);
   downstream_4->GetXaxis()->SetTitle("Time (ns)");
   downstream_4->GetXaxis()->SetNdivisions(5);
   downstream_4->SetMarkerStyle(20);
   downstream_4->SetMarkerSize(0.5);
   downstream_4->Draw();


   // retrieving information from the histogram

   Int_t jbins = downstream_4->GetNbinsX();
   Double_t xlo = downstream_4->GetBinLowEdge(1);
   Double_t width = downstream_4->GetBinWidth(1);
   printf ("jndim=%d, xlo=%f, width=%f\n",ndim,xlo,width);


   Double_t error_dn_4;
   Double_t time_dn_4;
   for(Int_t j=0;j<jbins;j++) {
      Double_t xbin = xlo + j*width;
      jj = (xbin-xlo)/width;

      Int_t nfitpts=3;
      Double_t datay[3], errory[3];
      Double_t datax[3];
      Double_t errorx[3]={0.025,0.025,0.025};
      Double_t binx1 = downstream_4->GetBinCenter(jj+1);
      datax[1] = 0;  // relative to binx1
      datay[1]= downstream_4->GetBinContent(jj+1) - bin_thresh;  // compute relative to (x-binx1,y-thresh)
      errory[1] = downstream_4->GetBinError(jj+1);
      // if (datay[1] > bin_thresh) {
      if (datay[1] > 0) {
         // calculate relative to middle bin
         datax[0] = downstream_4->GetBinCenter(jj) -binx1;
         datay[0] = downstream_4->GetBinContent(jj) - bin_thresh;
         errory[0] = downstream_4->GetBinError(jj);
         datax[2] = downstream_4->GetBinCenter(jj+2) -binx1;
         datay[2] = downstream_4->GetBinContent(jj+2) - bin_thresh;
         errory[2] = downstream_4->GetBinError(jj+2);

         // compute slope using data points before and after

          printf("datax %f %f %f\n",datax[0],datax[1],datax[2]);
          printf("datay %f %f %f\n",datay[0],datay[1],datay[2]);
          printf("errory %f %f %f\n",errory[0],errory[1],errory[2]);
          TGraphErrors *data = new TGraphErrors (nfitpts,datax,datay,errorx,errory);
          TFitResultPtr ptr = data->Fit("pol1","S","",-0.05,0.05);
          TMatrixDSym cov = ptr->GetCovarianceMatrix();       // access to covariance matrix
          cov.Print();
          Double_t determinant = cov.Determinant();
          // data->Draw("AP");

         if (determinant  != 0) {
              Double_t row1[2], row2[2];
              cov.ExtractRow(0,0,row1);
              cov.ExtractRow(1,0,row2);
              printf ("matrix=%f %f %f %f, determinant=%f\n",row1[0],row1[1],row2[0],row2[1],determinant);
              offset = ptr->Value(0);
              slope = ptr->Value(1);
              Doffset= sqrt(row1[0]);
              Dslope = sqrt(row2[1]);
              eps12 = row1[1];
              if (abs(slope) > 1.e-6) {
                 time_dn_4 = -offset/slope + binx1;
                 error_dn_4 = sqrt(Doffset*Doffset + (offset/slope)*(offset/slope)*Dslope*Dslope -(offset/slope)*eps12)/slope;
                 }
                 else {
                 time_dn_4 = binx1;
                 error_dn_4  = 0.025;
                 }
              printf (" slope =%f Dslope=%f offset=%f Doffset=%f binx1=%f time_dn_4 =%f error_dn_4 =%f\n",slope,Dslope,offset,Doffset,binx1,time_dn_4,error_dn_4);
              }
              else {
              printf ("*** Delta = 0\n",Delta);
              slope = 0;
              error_dn_4 = 0;
              time_dn_4 = 0;
              }
          j = jbins;
          }
   }

   sprintf(string,"Time=%.3f, Delta x=%.3f ns\n",time_dn_4,error_dn_4);
   t1 = new TLatex(0.18,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();


   /*c3->cd(5);
   c3_5->SetGridx();
   c3_5->SetGridy();
   c3_5->SetBorderMode(0);
   c3_5->SetFillColor(0);
   xmin = downstream_5->GetMean() - deltax;
   xmax = downstream_5->GetMean() + deltax;
   downstream_5->GetXaxis()->SetRangeUser(xmin,xmax);
   downstream_5->GetXaxis()->SetTitleSize(0.05);
   downstream_5->GetYaxis()->SetTitleSize(0.05);
   downstream_5->GetYaxis()->SetTitleOffset(1.5);
   downstream_5->GetXaxis()->SetTitle("Time (ns)");
   downstream_5->GetXaxis()->SetNdivisions(5);
   downstream_5->SetMarkerStyle(20);
   downstream_5->SetMarkerSize(0.5);
   downstream_5->Draw();

   c3->cd(6);
   c3_6->SetGridx();
   c3_6->SetGridy();
   c3_6->SetBorderMode(0);
   c3_6->SetFillColor(0);
   xmin = downstream_6->GetMean() - deltax;
   xmax = downstream_6->GetMean() + deltax;
   downstream_6->GetXaxis()->SetRangeUser(xmin,xmax);
   downstream_6->GetXaxis()->SetTitleSize(0.05);
   downstream_6->GetYaxis()->SetTitleSize(0.05);
   downstream_6->GetYaxis()->SetTitleOffset(1.5);
   downstream_6->GetXaxis()->SetTitle("Time (ns)");
   downstream_6->GetXaxis()->SetNdivisions(5);
   downstream_6->SetMarkerStyle(20);
   downstream_6->SetMarkerSize(0.5);
   downstream_6->Draw();*/

   /*TCanvas *c4 = new TCanvas("c4","c4 geom_dt",200,10,700,700);
   c4->SetBorderMode(0);
   c4->SetFillColor(0);
   c4->Divide(2,3);

   c4->cd(1);  
   c4_1->SetGridx();
   c4_1->SetGridy();
   c4_1->SetBorderMode(0);
   c4_1->SetFillColor(0);


   xmin = downstream_7->GetMean() - deltax;
   xmax = downstream_7->GetMean() + deltax;
   downstream_7->GetXaxis()->SetRangeUser(xmin,xmax);
   downstream_7->GetXaxis()->SetTitleSize(0.05);
   downstream_7->GetYaxis()->SetTitleSize(0.05);
   downstream_7->GetYaxis()->SetTitleOffset(1.5);
   downstream_7->GetXaxis()->SetTitle("Time (ns)");
   downstream_7->GetXaxis()->SetNdivisions(5);
   downstream_7->SetMarkerStyle(20);
   downstream_7->SetMarkerSize(0.5);
   downstream_7->SetMarkerColor(1);
   downstream_7->Draw();

   sprintf(string,"Eshower=%.2f GeV, #theta=%.0f deg\n",Eshower,theta);
   t1 = new TLatex(0.18,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c4->cd(2); 
   c4_2->SetGridx();
   c4_2->SetGridy();
   c4_2->SetBorderMode(0);
   c4_2->SetFillColor(0);
   xmin = downstream_8->GetMean() - deltax;
   xmax = downstream_8->GetMean() + deltax;
   downstream_8->GetXaxis()->SetRangeUser(xmin,xmax);
   downstream_8->GetXaxis()->SetTitleSize(0.05);
   downstream_8->GetYaxis()->SetTitleSize(0.05);
   downstream_8->GetYaxis()->SetTitleOffset(1.5);
   downstream_8->GetXaxis()->SetTitle("Time (ns)");
   downstream_8->GetXaxis()->SetNdivisions(5);
   downstream_8->SetMarkerStyle(20);
   downstream_8->SetMarkerSize(0.5);
   downstream_8->Draw();

   c4->cd(3);
   c4_3->SetGridx();
   c4_3->SetGridy();
   c4_3->SetBorderMode(0);
   c4_3->SetFillColor(0);
   xmin = downstream_9->GetMean() - deltax;
   xmax = downstream_9->GetMean() + deltax;
   downstream_9->GetXaxis()->SetRangeUser(xmin,xmax);
   downstream_9->GetXaxis()->SetTitleSize(0.05);
   downstream_9->GetYaxis()->SetTitleSize(0.05);
   downstream_9->GetYaxis()->SetTitleOffset(1.5);
   downstream_9->GetXaxis()->SetTitle("Time (ns)");
   downstream_9->GetXaxis()->SetNdivisions(5);
   downstream_9->SetMarkerStyle(20);
   downstream_9->SetMarkerSize(0.5);
   downstream_9->Draw();

   c4->cd(4);
   c4_4->SetGridx();
   c4_4->SetGridy();
   c4_4->SetBorderMode(0);
   c4_4->SetFillColor(0);
   xmin = downstream_10->GetMean() - deltax;
   xmax = downstream_10->GetMean() + deltax;
   downstream_10->GetXaxis()->SetRangeUser(xmin,xmax);
   downstream_10->GetXaxis()->SetTitleSize(0.05);
   downstream_10->GetYaxis()->SetTitleSize(0.05);
   downstream_10->GetYaxis()->SetTitleOffset(1.5);
   downstream_10->GetXaxis()->SetTitle("Time (ns)");
   downstream_10->GetXaxis()->SetNdivisions(5);
   downstream_10->SetMarkerStyle(20);
   downstream_10->SetMarkerSize(0.5);
   downstream_10->Draw();*/


   TCanvas *c5 = new TCanvas("c5","c5 geom_dt",200,10,700,700);
   c5->SetBorderMode(0);
   c5->SetFillColor(0);
   c5->Divide(2,2);

   c5->cd(1);  
   c5_1->SetGridx();
   c5_1->SetGridy();
   c5_1->SetBorderMode(0);
   c5_1->SetFillColor(0);

   // diff_1->SetTitle("Diff");
   xmin = diff_1->GetMean() - deltax;
   xmax = diff_1->GetMean() + deltax;
   diff_1->GetXaxis()->SetRangeUser(xmin,xmax);
   // diff_1->GetYaxis()->SetRangeUser(ymin,ymax);
   diff_1->GetXaxis()->SetTitleSize(0.05);
   diff_1->GetYaxis()->SetTitleSize(0.05);
   diff_1->GetYaxis()->SetTitleOffset(1.5);
   diff_1->GetXaxis()->SetTitle("Time (ns)");
   // diff_1->GetYaxis()->SetTitle("Npe");
   diff_1->GetXaxis()->SetNdivisions(5);
   diff_1->SetMarkerStyle(20);
   diff_1->SetMarkerSize(0.5);
   diff_1->SetMarkerColor(1);
   diff_1->Draw();

   sprintf(string,"Eshower=%.2f GeV, #theta=%.0f deg\n",Eshower,theta);
   t1 = new TLatex(0.18,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c5->cd(2); 
   c5_2->SetGridx();
   c5_2->SetGridy();
   c5_2->SetBorderMode(0);
   c5_2->SetFillColor(0);
   // c5_2->SetLogx();
   xmin = diff_2->GetMean() - deltax;
   xmax = diff_2->GetMean() + deltax;
   diff_2->GetXaxis()->SetRangeUser(xmin,xmax);
   diff_2->GetXaxis()->SetTitleSize(0.05);
   diff_2->GetYaxis()->SetTitleSize(0.05);
   diff_2->GetYaxis()->SetTitleOffset(1.5);
   diff_2->GetXaxis()->SetTitle("Time (ns)");
   diff_2->GetXaxis()->SetNdivisions(5);
   diff_2->SetMarkerStyle(20);
   diff_2->SetMarkerSize(0.5);
   diff_2->Draw();

   c5->cd(3);
   c5_3->SetGridx();
   c5_3->SetGridy();
   c5_3->SetBorderMode(0);
   c5_3->SetFillColor(0);
   xmin = diff_3->GetMean() - deltax;
   xmax = diff_3->GetMean() + deltax;
   diff_3->GetXaxis()->SetRangeUser(xmin,xmax);
   diff_3->GetXaxis()->SetTitleSize(0.05);
   diff_3->GetYaxis()->SetTitleSize(0.05);
   diff_3->GetYaxis()->SetTitleOffset(1.5);
   diff_3->GetXaxis()->SetTitle("Time (ns)");
   diff_3->GetXaxis()->SetNdivisions(5);
   diff_3->SetMarkerStyle(20);
   diff_3->SetMarkerSize(0.5);
   diff_3->Draw();

   c5->cd(4);
   c5_4->SetGridx();
   c5_4->SetGridy();
   c5_4->SetBorderMode(0);
   c5_4->SetFillColor(0);
   xmin = diff_4->GetMean() - deltax;
   xmax = diff_4->GetMean() + deltax;
   diff_4->GetXaxis()->SetRangeUser(xmin,xmax);
   diff_4->GetXaxis()->SetTitleSize(0.05);
   diff_4->GetYaxis()->SetTitleSize(0.05);
   diff_4->GetYaxis()->SetTitleOffset(1.5);
   diff_4->GetXaxis()->SetTitle("Time (ns)");
   diff_4->GetXaxis()->SetNdivisions(5);
   diff_4->SetMarkerStyle(20);
   diff_4->SetMarkerSize(0.5);
   diff_4->Draw();

   /*c5->cd(5);
   c5_5->SetGridx();
   c5_5->SetGridy();
   c5_5->SetBorderMode(0);
   c5_5->SetFillColor(0);
   xmin = diff_5->GetMean() - deltax;
   xmax = diff_5->GetMean() + deltax;
   diff_5->GetXaxis()->SetRangeUser(xmin,xmax);
   diff_5->GetXaxis()->SetTitleSize(0.05);
   diff_5->GetYaxis()->SetTitleSize(0.05);
   diff_5->GetYaxis()->SetTitleOffset(1.5);
   diff_5->GetXaxis()->SetTitle("Time (ns)");
   diff_5->GetXaxis()->SetNdivisions(5);
   diff_5->SetMarkerStyle(20);
   diff_5->SetMarkerSize(0.5);
   diff_5->Draw();

   c5->cd(6);
   c5_6->SetGridx();
   c5_6->SetGridy();
   c5_6->SetBorderMode(0);
   c5_6->SetFillColor(0);
   xmin = diff_6->GetMean() - deltax;
   xmax = diff_6->GetMean() + deltax;
   diff_6->GetXaxis()->SetRangeUser(xmin,xmax);
   diff_6->GetXaxis()->SetTitleSize(0.05);
   diff_6->GetYaxis()->SetTitleSize(0.05);
   diff_6->GetYaxis()->SetTitleOffset(1.5);
   diff_6->GetXaxis()->SetTitle("Time (ns)");
   diff_6->GetXaxis()->SetNdivisions(5);
   diff_6->SetMarkerStyle(20);
   diff_6->SetMarkerSize(0.5);
   diff_6->Draw();*/

   /*TCanvas *c6 = new TCanvas("c6","c6 geom_dt",200,10,700,700);
   c6->SetBorderMode(0);
   c6->SetFillColor(0);
   c6->Divide(2,3);

   c6->cd(1);  
   c6_1->SetGridx();
   c6_1->SetGridy();
   c6_1->SetBorderMode(0);
   c6_1->SetFillColor(0);

   xmin = diff_7->GetMean() - deltax;
   xmax = diff_7->GetMean() + deltax;
   diff_7->GetXaxis()->SetRangeUser(xmin,xmax);
   diff_7->GetXaxis()->SetTitleSize(0.05);
   diff_7->GetYaxis()->SetTitleSize(0.05);
   diff_7->GetYaxis()->SetTitleOffset(1.5);
   diff_7->GetXaxis()->SetTitle("Time (ns)");
   diff_7->GetXaxis()->SetNdivisions(5);
   diff_7->SetMarkerStyle(20);
   diff_7->SetMarkerSize(0.5);
   diff_7->SetMarkerColor(1);
   diff_7->Draw();

   sprintf(string,"Eshower=%.2f GeV, #theta=%.0f deg\n",Eshower,theta);
   t1 = new TLatex(0.18,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c6->cd(2); 
   c6_2->SetGridx();
   c6_2->SetGridy();
   c6_2->SetBorderMode(0);
   c6_2->SetFillColor(0);
   xmin = diff_8->GetMean() - deltax;
   xmax = diff_8->GetMean() + deltax;
   diff_8->GetXaxis()->SetRangeUser(xmin,xmax);
   diff_8->GetXaxis()->SetTitleSize(0.05);
   diff_8->GetYaxis()->SetTitleSize(0.05);
   diff_8->GetYaxis()->SetTitleOffset(1.5);
   diff_8->GetXaxis()->SetTitle("Time (ns)");
   diff_8->GetXaxis()->SetNdivisions(5);
   diff_8->SetMarkerStyle(20);
   diff_8->SetMarkerSize(0.5);
   diff_8->Draw();

   c6->cd(3);
   c6_3->SetGridx();
   c6_3->SetGridy();
   c6_3->SetBorderMode(0);
   c6_3->SetFillColor(0);
   xmin = diff_9->GetMean() - deltax;
   xmax = diff_9->GetMean() + deltax;
   diff_9->GetXaxis()->SetRangeUser(xmin,xmax);
   diff_9->GetXaxis()->SetTitleSize(0.05);
   diff_9->GetYaxis()->SetTitleSize(0.05);
   diff_9->GetYaxis()->SetTitleOffset(1.5);
   diff_9->GetXaxis()->SetTitle("Time (ns)");
   diff_9->GetXaxis()->SetNdivisions(5);
   diff_9->SetMarkerStyle(20);
   diff_9->SetMarkerSize(0.5);
   diff_9->Draw();

   c6->cd(4);
   c6_4->SetGridx();
   c6_4->SetGridy();
   c6_4->SetBorderMode(0);
   c6_4->SetFillColor(0);
   xmin = diff_10->GetMean() - deltax;
   xmax = diff_10->GetMean() + deltax;
   diff_10->GetXaxis()->SetRangeUser(xmin,xmax);
   diff_10->GetXaxis()->SetTitleSize(0.05);
   diff_10->GetYaxis()->SetTitleSize(0.05);
   diff_10->GetYaxis()->SetTitleOffset(1.5);
   diff_10->GetXaxis()->SetTitle("Time (ns)");
   diff_10->GetXaxis()->SetNdivisions(5);
   diff_10->SetMarkerStyle(20);
   diff_10->SetMarkerSize(0.5);
   diff_10->Draw();*/


   TCanvas *c7 = new TCanvas("c7","c7 geom_dt",200,10,700,700);
   c7->SetBorderMode(0);
   c7->SetFillColor(0);
   c7->Divide(2,2);

   c7->cd(1);  
   c7_1->SetGridx();
   c7_1->SetGridy();
   c7_1->SetBorderMode(0);
   c7_1->SetFillColor(0);

   // sum_1->SetTitle("Sum");
   xmin = sum_1->GetMean() - deltax;
   xmax = sum_1->GetMean() + deltax;
   sum_1->GetXaxis()->SetRangeUser(xmin,xmax);
   // sum_1->GetYaxis()->SetRangeUser(ymin,ymax);
   sum_1->GetXaxis()->SetTitleSize(0.05);
   sum_1->GetYaxis()->SetTitleSize(0.05);
   sum_1->GetYaxis()->SetTitleOffset(1.5);
   sum_1->GetXaxis()->SetTitle("Time (ns)");
   // sum_1->GetYaxis()->SetTitle("Npe");
   sum_1->GetXaxis()->SetNdivisions(5);
   sum_1->SetMarkerStyle(20);
   sum_1->SetMarkerSize(0.5);
   sum_1->SetMarkerColor(1);
   sum_1->Draw();

   sprintf(string,"Eshower=%.2f GeV, #theta=%.0f deg\n",Eshower,theta);
   t1 = new TLatex(0.18,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c7->cd(2); 
   c7_2->SetGridx();
   c7_2->SetGridy();
   c7_2->SetBorderMode(0);
   c7_2->SetFillColor(0);
   // c7_2->SetLogx();
   xmin = sum_2->GetMean() - deltax;
   xmax = sum_2->GetMean() + deltax;
   sum_2->GetXaxis()->SetRangeUser(xmin,xmax);
   sum_2->GetXaxis()->SetTitleSize(0.05);
   sum_2->GetYaxis()->SetTitleSize(0.05);
   sum_2->GetYaxis()->SetTitleOffset(1.5);
   sum_2->GetXaxis()->SetTitle("Time (ns)");
   sum_2->GetXaxis()->SetNdivisions(5);
   sum_2->SetMarkerStyle(20);
   sum_2->SetMarkerSize(0.5);
   sum_2->Draw();

   c7->cd(3);
   c7_3->SetGridx();
   c7_3->SetGridy();
   c7_3->SetBorderMode(0);
   c7_3->SetFillColor(0);
   xmin = sum_3->GetMean() - deltax;
   xmax = sum_3->GetMean() + deltax;
   sum_3->GetXaxis()->SetRangeUser(xmin,xmax);
   sum_3->GetXaxis()->SetTitleSize(0.05);
   sum_3->GetYaxis()->SetTitleSize(0.05);
   sum_3->GetYaxis()->SetTitleOffset(1.5);
   sum_3->GetXaxis()->SetTitle("Time (ns)");
   sum_3->GetXaxis()->SetNdivisions(5);
   sum_3->SetMarkerStyle(20);
   sum_3->SetMarkerSize(0.5);
   sum_3->Draw();

   c7->cd(4);
   c7_4->SetGridx();
   c7_4->SetGridy();
   c7_4->SetBorderMode(0);
   c7_4->SetFillColor(0);
   xmin = sum_4->GetMean() - deltax;
   xmax = sum_4->GetMean() + deltax;
   sum_4->GetXaxis()->SetRangeUser(xmin,xmax);
   sum_4->GetXaxis()->SetTitleSize(0.05);
   sum_4->GetYaxis()->SetTitleSize(0.05);
   sum_4->GetYaxis()->SetTitleOffset(1.5);
   sum_4->GetXaxis()->SetTitle("Time (ns)");
   sum_4->GetXaxis()->SetNdivisions(5);
   sum_4->SetMarkerStyle(20);
   sum_4->SetMarkerSize(0.5);
   sum_4->Draw();

   /*c7->cd(5);
   c7_5->SetGridx();
   c7_5->SetGridy();
   c7_5->SetBorderMode(0);
   c7_5->SetFillColor(0);
   xmin = sum_5->GetMean() - deltax;
   xmax = sum_5->GetMean() + deltax;
   sum_5->GetXaxis()->SetRangeUser(xmin,xmax);
   sum_5->GetXaxis()->SetTitleSize(0.05);
   sum_5->GetYaxis()->SetTitleSize(0.05);
   sum_5->GetYaxis()->SetTitleOffset(1.5);
   sum_5->GetXaxis()->SetTitle("Time (ns)");
   sum_5->GetXaxis()->SetNdivisions(5);
   sum_5->SetMarkerStyle(20);
   sum_5->SetMarkerSize(0.5);
   sum_5->Draw();

   c7->cd(6);
   c7_6->SetGridx();
   c7_6->SetGridy();
   c7_6->SetBorderMode(0);
   c7_6->SetFillColor(0);
   xmin = sum_6->GetMean() - deltax;
   xmax = sum_6->GetMean() + deltax;
   sum_6->GetXaxis()->SetRangeUser(xmin,xmax);
   sum_6->GetXaxis()->SetTitleSize(0.05);
   sum_6->GetYaxis()->SetTitleSize(0.05);
   sum_6->GetYaxis()->SetTitleOffset(1.5);
   sum_6->GetXaxis()->SetTitle("Time (ns)");
   sum_6->GetXaxis()->SetNdivisions(5);
   sum_6->SetMarkerStyle(20);
   sum_6->SetMarkerSize(0.5);
   sum_6->Draw();*/

   /*TCanvas *c8 = new TCanvas("c8","c8 geom_dt",200,10,700,700);
   c8->SetBorderMode(0);
   c8->SetFillColor(0);
   c8->Divide(2,3);

   c8->cd(1);  
   c8_1->SetGridx();
   c8_1->SetGridy();
   c8_1->SetBorderMode(0);
   c8_1->SetFillColor(0);


   xmin = sum_7->GetMean() - deltax;
   xmax = sum_7->GetMean() + deltax;
   sum_7->GetXaxis()->SetRangeUser(xmin,xmax);
   sum_7->GetXaxis()->SetTitleSize(0.05);
   sum_7->GetYaxis()->SetTitleSize(0.05);
   sum_7->GetYaxis()->SetTitleOffset(1.5);
   sum_7->GetXaxis()->SetTitle("Time (ns)");
   sum_7->GetXaxis()->SetNdivisions(5);
   sum_7->SetMarkerStyle(20);
   sum_7->SetMarkerSize(0.5);
   sum_7->SetMarkerColor(1);
   sum_7->Draw();

   sprintf(string,"Eshower=%.2f GeV, #theta=%.0f deg\n",Eshower,theta);
   t1 = new TLatex(0.18,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c8->cd(2); 
   c8_2->SetGridx();
   c8_2->SetGridy();
   c8_2->SetBorderMode(0);
   c8_2->SetFillColor(0);
   xmin = sum_8->GetMean() - deltax;
   xmax = sum_8->GetMean() + deltax;
   sum_8->GetXaxis()->SetRangeUser(xmin,xmax);
   sum_8->GetXaxis()->SetTitleSize(0.05);
   sum_8->GetYaxis()->SetTitleSize(0.05);
   sum_8->GetYaxis()->SetTitleOffset(1.5);
   sum_8->GetXaxis()->SetTitle("Time (ns)");
   sum_8->GetXaxis()->SetNdivisions(5);
   sum_8->SetMarkerStyle(20);
   sum_8->SetMarkerSize(0.5);
   sum_8->Draw();

   c8->cd(3);
   c8_3->SetGridx();
   c8_3->SetGridy();
   c8_3->SetBorderMode(0);
   c8_3->SetFillColor(0);
   xmin = sum_9->GetMean() - deltax;
   xmax = sum_9->GetMean() + deltax;
   sum_9->GetXaxis()->SetRangeUser(xmin,xmax);
   sum_9->GetXaxis()->SetTitleSize(0.05);
   sum_9->GetYaxis()->SetTitleSize(0.05);
   sum_9->GetYaxis()->SetTitleOffset(1.5);
   sum_9->GetXaxis()->SetTitle("Time (ns)");
   sum_9->GetXaxis()->SetNdivisions(5);
   sum_9->SetMarkerStyle(20);
   sum_9->SetMarkerSize(0.5);
   sum_9->Draw();

   c8->cd(4);
   c8_4->SetGridx();
   c8_4->SetGridy();
   c8_4->SetBorderMode(0);
   c8_4->SetFillColor(0);
   xmin = sum_10->GetMean() - deltax;
   xmax = sum_10->GetMean() + deltax;
   sum_10->GetXaxis()->SetRangeUser(xmin,xmax);
   sum_10->GetXaxis()->SetTitleSize(0.05);
   sum_10->GetYaxis()->SetTitleSize(0.05);
   sum_10->GetYaxis()->SetTitleOffset(1.5);
   sum_10->GetXaxis()->SetTitle("Time (ns)");
   sum_10->GetXaxis()->SetNdivisions(5);
   sum_10->SetMarkerStyle(20);
   sum_10->SetMarkerSize(0.5);
   sum_10->Draw();*/


  // open output file

  Double_t vfactor = 1.00;  // empirical factor to account for spread of shower effectively reducing velocity time
  sprintf (filename,"geom_dt_2-%d_%.0f_%.0f_%.0f.txt",nlayers,Eshower*1000,theta,bin_thresh);
  FILE *out;
  out = fopen(filename,"w+");

  // add offset to sum
   Double_t offset_sum = 0;
   Double_t zoffset=2.06/2;

  Double_t tsum_1 = 0.5*(time_up_1+time_dn_1) + offset_sum;
  Double_t tsum_err_1 =  0.5*sqrt(error_up_1*error_up_1 + error_dn_1*error_dn_1);
  Double_t tdiff_1=0.5*(time_up_1-time_dn_1)*veff;
  Double_t tdiff_err_1 = 0.5*sqrt(error_up_1*error_up_1 + error_dn_1*error_dn_1)*veff;
  Double_t tpos_1 = (z_downstream-z_upstream)/(2*veff*vfactor);   // time relative to impact time
  Double_t zpos_1 = zimpact + (zoffset)/(tan(theta*pi/180.)) - (z_downstream+z_upstream)/2.;
  fprintf (out,"Upstream 1: %.3f  +/- %.3f, Downstream 1: %.3f +/- %.3f, Mean 1: %.3f +/- %.3f, Tpos=%.3f ns, Half Diff 1: %.3f +/- %.3f, zpos: %.3f cm, Dmean: %.3f, Dpos: %.2f\n",time_up_1,error_up_1,time_dn_1,error_dn_1,tsum_1,tsum_err_1,tpos_1,tdiff_1,tdiff_err_1,zpos_1,tsum_1-tpos_1,tdiff_1-zpos_1);
  fprintf (out,"%.3f   %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f %.3f %.2f\n",time_up_1,error_up_1,time_dn_1,error_dn_1,tsum_1,tsum_err_1,tpos_1,tdiff_1,tdiff_err_1,zpos_1,tsum_1-tpos_1,tdiff_1-zpos_1);

  Double_t tsum_2 = 0.5*(time_up_2+time_dn_2) + offset_sum;
  Double_t tsum_err_2 =  0.5*sqrt(error_up_2*error_up_2 + error_dn_2*error_dn_2);
  Double_t tdiff_2=0.5*(time_up_2-time_dn_2)*veff;
  Double_t tdiff_err_2 = 0.5*sqrt(error_up_2*error_up_2 + error_dn_2*error_dn_2)*veff;
  Double_t tpos_2 = path[0]/c + (z_downstream-z_upstream)/(2*veff*vfactor);   // time relative to impact time
  Double_t zpos_2 = zimpact  + (zoffset+height[0])/(tan(theta*pi/180.))  - (z_downstream+z_upstream)/2.;
  fprintf (out,"%.3f   %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f %.3f %.2f\n",time_up_2,error_up_2,time_dn_2,error_dn_2,tsum_2,tsum_err_2,tpos_2,tdiff_2,tdiff_err_2,zpos_2,tsum_2-tpos_2,tdiff_2-zpos_2);
  // fprintf (out,"Upstream 2: %.3f  +/- %.3f, Downstream 2: %.3f +/- %.3f, Mean 2: %.3f +/- %.3f, Tpos=%.3f ns, Half Diff 2: %.3f +/- %.3f, zpos: %.3f cm\n",time_up_2,error_up_2,time_dn_2,error_dn_2,tsum_2,tsum_err_2,tpos_2,tdiff_2,tdiff_err_2,zpos_2);

  Double_t tsum_3 = 0.5*(time_up_3+time_dn_3) + offset_sum;
  Double_t tsum_err_3 =  0.5*sqrt(error_up_3*error_up_3 + error_dn_3*error_dn_3);
  Double_t tdiff_3=0.5*(time_up_3-time_dn_3)*veff;
  Double_t tdiff_err_3 = 0.5*sqrt(error_up_3*error_up_3 + error_dn_3*error_dn_3)*veff;
  Double_t tpos_3 =  path[1]/c +  (z_downstream-z_upstream)/(2*veff*vfactor);   // time relative to impact time
  Double_t zpos_3 = zimpact  +  (zoffset+height[0]+height[1])/(tan(theta*pi/180.)) - (z_downstream+z_upstream)/2.;
  fprintf (out,"%.3f   %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f %.3f %.2f\n",time_up_3,error_up_3,time_dn_3,error_dn_3,tsum_3,tsum_err_3,tpos_3,tdiff_3,tdiff_err_3,zpos_3,tsum_3-tpos_3,tdiff_3-zpos_3);
  // fprintf (out,"Upstream 3: %.3f  +/- %.3f, Downstream 3: %.3f +/- %.3f, Mean 3: %.3f +/- %.3f, Tpos=%.3f ns, Half Diff 3: %.3f +/- %.3f, zpos: %.3f cm\n",time_up_3,error_up_3,time_dn_3,error_dn_3,tsum_3,tsum_err_3,tpos_3,tdiff_3,tdiff_err_3,zpos_3);

  Double_t tsum_4 = 0.5*(time_up_4+time_dn_4) + offset_sum;
  Double_t tsum_err_4 =  0.5*sqrt(error_up_4*error_up_4 + error_dn_4*error_dn_4);
  Double_t tdiff_4=0.5*(time_up_4-time_dn_4)*veff;
  Double_t tdiff_err_4 = 0.5*sqrt(error_up_4*error_up_4 + error_dn_4*error_dn_4)*veff;
  Double_t tpos_4 =  path[2]/c + (z_downstream-z_upstream)/(2*veff*vfactor);   // time relative to impact time
  Double_t zpos_4 = zimpact  + (zoffset+height[0]+height[1]+height[2])/(tan(theta*pi/180.))  - (z_downstream+z_upstream)/2.;
  fprintf (out,"%.3f   %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f %.3f %.2f\n",time_up_4,error_up_4,time_dn_4,error_dn_4,tsum_4,tsum_err_4,tpos_4,tdiff_4,tdiff_err_4,zpos_4,tsum_4-tpos_4,tdiff_4-zpos_4);
  // fprintf (out,"Upstream 4: %.3f  +/- %.3f, Downstream 4: %.3f +/- %.3f, Mean 4: %.3f +/- %.3f, Tpos=%.3f ns, Half Diff 4: %.3f +/- %.3f, zpos: %.3f cm\n",time_up_4,error_up_4,time_dn_4,error_dn_4,tsum_4,tsum_err_4,tpos_4,tdiff_4,tdiff_err_4,zpos_4);
  fclose (out);


   // output canvases to file
   char filename1[80];
   char filename2[80];

   c0->SaveAs("Shower_profile.pdf");
   c0->SaveAs("Shower_profile.png");
   sprintf (filename1,"geom_dt_2-%d_%.0f_%.0f_%.0f_k0.pdf",nlayers,Eshower*1000,theta,bin_thresh);
   sprintf (filename2,"geom_dt_2-%d_%.0f_%.0f_%.0f_k0.png",nlayers,Eshower*1000,theta,bin_thresh);
   k0->SaveAs(filename1);
   k0->SaveAs(filename2);
   sprintf (filename1,"geom_dt_2-%d_%.0f_%.0f_%.0f_c1.pdf",nlayers,Eshower*1000,theta,bin_thresh);
   sprintf (filename2,"geom_dt_2-%d_%.0f_%.0f_%.0f_c1.png",nlayers,Eshower*1000,theta,bin_thresh);
   c1->SaveAs(filename1);
   c1->SaveAs(filename2);
   sprintf (filename1,"geom_dt_2-%d_%.0f_%.0f_%.0f_c2.pdf",nlayers,Eshower*1000,theta,bin_thresh);
   sprintf (filename2,"geom_dt_2-%d_%.0f_%.0f_%.0f_c2.png",nlayers,Eshower*1000,theta,bin_thresh);
   /*c2->SaveAs(filename1);
   c2->SaveAs(filename2);*/
   sprintf (filename1,"geom_dt_2-%d_%.0f_%.0f_%.0f_c3.pdf",nlayers,Eshower*1000,theta,bin_thresh);
   sprintf (filename2,"geom_dt_2-%d_%.0f_%.0f_%.0f_c3.png",nlayers,Eshower*1000,theta,bin_thresh);
   c3->SaveAs(filename1);
   c3->SaveAs(filename2);
   sprintf (filename1,"geom_dt_2-%d_%.0f_%.0f_%.0f_c4.pdf",nlayers,Eshower*1000,theta,bin_thresh);
   sprintf (filename2,"geom_dt_2-%d_%.0f_%.0f_%.0f_c4.png",nlayers,Eshower*1000,theta,bin_thresh);
   /*c4->SaveAs(filename1);
   c4->SaveAs(filename2);*/
   sprintf (filename1,"geom_dt_2-%d_%.0f_%.0f_%.0f_c5.pdf",nlayers,Eshower*1000,theta,bin_thresh);
   sprintf (filename2,"geom_dt_2-%d_%.0f_%.0f_%.0f_c5.png",nlayers,Eshower*1000,theta,bin_thresh);
   c5->SaveAs(filename1);
   c5->SaveAs(filename2);
   sprintf (filename1,"geom_dt_2-%d_%.0f_%.0f_%.0f_c6.pdf",nlayers,Eshower*1000,theta,bin_thresh);
   sprintf (filename2,"geom_dt_2-%d_%.0f_%.0f_%.0f_c6.png",nlayers,Eshower*1000,theta,bin_thresh);
   /*c6->SaveAs(filename1);
   c6->SaveAs(filename2);*/
   sprintf (filename1,"geom_dt_2-%d_%.0f_%.0f_%.0f_c7.pdf",nlayers,Eshower*1000,theta,bin_thresh);
   sprintf (filename2,"geom_dt_2-%d_%.0f_%.0f_%.0f_c7.png",nlayers,Eshower*1000,theta,bin_thresh);
   c7->SaveAs(filename1);
   c7->SaveAs(filename2);
   sprintf (filename1,"geom_dt_2-%d_%.0f_%.0f_%.0f_c8.pdf",nlayers,Eshower*1000,theta,bin_thresh);
   sprintf (filename2,"geom_dt_2-%d_%.0f_%.0f_%.0f_c8.png",nlayers,Eshower*1000,theta,bin_thresh);
   /*c8->SaveAs(filename1);
   c8->SaveAs(filename2);*/
}
Double_t shower_profile (Double_t *x, Double_t *par) {

     // return shower profile in number of photoelectrons/GeV/side in the Bcal

     Double_t a=par[0];
     Double_t b=par[1];
     Double_t npe_per_GeV_per_side=par[2];
     Double_t t=x[0]; 

     char string[256];

     Double_t func=0;
    
     func = a + b*t;
     func = npe_per_GeV_per_side * 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;
}