#include <iostream>
using namespace std;

void tof_calibration_read(void)
{
// 
// Input root tree generated by genr8_tree.C and histogram variables.
//

#include <TRandom.h>
#include <TTree.h>

gROOT->Reset();

gStyle->SetPalette(1,0);
// gStyle->SetOptStat(kFALSE);;
gStyle->SetOptStat(kTRUE);
gStyle->SetOptStat(1111111);
gStyle->SetOptFit(kTRUE);
gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
// gStyle->SetFillColor(0);
//

   char string[256];
   char filename[132];
   char name[10];
   char cut[132];
   Int_t j,jj;
   Double_t pi=3.14159;
   
   #define npts 55;
   Int_t nevents=10000;
   #define ntop2 8;
   Float_t content[8];

   // 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;

   // open output files

   Int_t eventno;

  // define histograms

   Double_t xmin=-120;
   Double_t xmax=120;
   Double_t ymin=-120;
   Double_t ymax=120;
   Double_t tmin=-10;
   Double_t tmax=40;
   Double_t emin=0;
   Double_t emax=20;
   Double_t vmin=0;
   Double_t vmax=2;
   Int_t nbinsx = 100;
   Int_t nbinsy = 100;
   Int_t nbinst = 200;
   Int_t nbinse = 100;
   Int_t nbinsv = 100;

  
   TH1F *t0 =   new TH1F("t0","Generated t0",nbinst,tmin-5,tmax+5);
   TH2F *genxy =   new TH2F("genxy","y vs x",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
   TH1F *dE1 =   new TH1F("dE1","Energy Loss Horizonal",nbinse,emin,emax);
   TH1F *dE2 =   new TH1F("dE2","Energy Loss Vertical",nbinse,emin,emax);

   // reconstructed quantities

   TH2F *t0H =   new TH2F("t0H","0.5*(t_R + t_L) vs t0",nbinst,tmin,tmax,nbinst,tmin,tmax);
   TH2F *t0V =   new TH2F("t0V","0.5*(t_R + t_L) vs t0",nbinst,tmin,tmax,nbinst,tmin,tmax);
   TH2F *xraw =   new TH2F("xraw","0.5*(t_R - t_L) vs x",nbinsx,xmin,xmax,nbinsx,xmin,xmax);
   TH2F *yraw =   new TH2F("yraw","0.5*(t_T - t_B) vs y",nbinsy,ymin,ymax,nbinsy,ymin,ymax);

   // same for times corrected for time-walk
  // Note: change title of "cor" -> "naive" to distinguish from Minuit fit. Leave variables alone.

   TH2F *t0corH =   new TH2F("t0naiveH","0.5*(tnaive_R + tnaive_L) vs t0",nbinst,tmin,tmax,nbinst,tmin,tmax);
   TH2F *t0corV =   new TH2F("t0naiveV","0.5*(tnaive_T + tnaive_B) vs t0",nbinst,tmin,tmax,nbinst,tmin,tmax);
   TH2F *xcor =   new TH2F("xnaive","0.5*(tnaive_R - tnaive_L) vs x",nbinsx,xmin,xmax,nbinsx,xmin,xmax);
   TH2F *ycor =   new TH2F("ynaive","0.5*(tnaive_T - tnaive_B) vs y",nbinsy,ymin,ymax,nbinsy,ymin,ymax);

   TH2F *wL = new TH2F("wL","t_L vs P_L",nbinsv,vmin,vmax,nbinst,5,15);
   TH2F *wR = new TH2F("wR","t_R vs P_R",nbinsv,vmin,vmax,nbinst,5,15);
   TH2F *wT = new TH2F("wT","t_T vs P_T",nbinsv,vmin,vmax,nbinst,5,15);
   TH2F *wB= new TH2F("wB","t_B vs P_B",nbinsv,vmin,vmax,nbinst,5,15);

   TH1F *sigtrawH = new TH1F("sigtraw","0.5*(t_R + t_L) - t0",nbinst,10,20);
   TH1F *sigtrawV = new TH1F("sigtraw","0.5*(t_T + t_B) - t0",nbinst,10,20);

   TH1F *sigteasyH = new TH1F("sigt naive cor","0.5*(tnaive_R + tnaive_L) - t0",nbinst,6,10);
   TH1F *sigteasyV = new TH1F("sigt naive cor","0.5*(tnaive_T + tnaive_B) - t0",nbinst,6,10);
   TH1F *sigtcorH = new TH1F("sigtcor","0.5*(tcor_R + tcor_L) - t0",nbinst,12,16);
   TH1F *sigtcorV = new TH1F("sigtcor","0.5*(tcor_T + tcor_B) - t0",nbinst,12,16);

   TH1F *sigtgenH = new TH1F("sigtgen","0.5*(tgen_R + tgen_L) - t0",nbinst,23,27);
   TH1F *sigtgenV = new TH1F("sigtgen","0.5*(tgen_T + tgen_B) - t0",nbinst,23,27);

  
   // read root tree 

   sprintf(filename,"tof_tree.root");
   printf ("filename = %s\n",filename);
   TFile *tfile = new TFile(filename);
   TTree *gentof = (TTree *) gROOT->FindObject("gentof");

   gentof->SetBranchAddress("Events",&Events.dE_V);
   
   Int_t nentries = gentof->GetEntries();
   for (Int_t j=0; j<nentries && j<nevents; j++){
      gentof->GetEntry(j);
      // printf ("nentries=%d, j=%d, Eb=%f, mg1g22=%f, mKg12=%f t=%f\n",nentries,j,Events.Eb,Events.mg1g22,Events.t);

         t0->Fill(Events.t0_H);
         t0->Fill(Events.t0_V);
         genxy->Fill(Events.x,Events.y);
         dE1->Fill(Events.dE_H);
         dE2->Fill(Events.dE_V);

         // compute derived quantities of interest

        t0H->Fill(Events.t0_H, 0.5*(Events.t_L+Events.t_R));
        t0V->Fill(Events.t0_V, 0.5*(Events.t_B+Events.t_T));
        xraw->Fill(Events.x, 0.5*(Events.t_L - Events.t_R)*Events.veff);
        yraw->Fill(Events.y, 0.5*(Events.t_B - Events.t_T)*Events.veff);

        t0corH->Fill(Events.t0_H, 0.5*(Events.tcor_L+Events.tcor_R));
        t0corV->Fill(Events.t0_V, 0.5*(Events.tcor_B+Events.tcor_T));
        xcor->Fill(Events.x, 0.5*(Events.tcor_L - Events.tcor_R)*Events.veff);
        ycor->Fill(Events.y, 0.5*(Events.tcor_B - Events.tcor_T)*Events.veff);

       sigtrawH->Fill(0.5*(Events.t_L+Events.t_R) - Events.t0_H);
       sigtrawV->Fill(0.5*(Events.t_T+Events.t_B) - Events.t0_V);
       sigteasyH->Fill(0.5*(Events.tcor_L+Events.tcor_R) - Events.t0_H);
       sigteasyV->Fill(0.5*(Events.tcor_T+Events.tcor_B) - Events.t0_V);

      Double_t xx[8]={6.95868,6.95868,6.95868,6.95868,1.36879,1.36879,1.36879,1.36879}; // constants at x0=0+/-3, y0=0+/-3 1000 event fit
      // Double_t xx[8]={6.94243,7.38598,-7.10287e+06,9.15008,0.924641,1.89558,-6.85158e+06,0.792373}; // constants at x0=0+/-120, y0=60+/-3 500 event fit
      // Double_t xx[8]={7.82661,10.724,8.36327,6.1211e+07,0.502474,4.67712,0.476978,4.32326e+07}; // constants at x0=0+/-120, y0=60+/-3 1000 event fit
      // Double_t xx[8]={7.01895,7.3892,9.0154,2080.63,0.688212,2.07727,0.751238,1629.07}; // constants at x0=0+/-120, y0=60+/-3 100 event fit
      // Double_t xx[8]={-2.84873,15.4672,9.85997,719.957,0.31379,0.86346,0.995294,724.988}; // constants at x0=0+/-3, y0=60+/-3
      // Double_t xx[8]={47.5807,-62.8133,-23.5626,610.715,0.311166,0.0969951,-0.012414,198.725};
      Double_t walk0L = xx[0];
      Double_t walk0R = xx[1];
      Double_t walk0B = xx[2];
      Double_t walk0T = xx[3];
      Double_t walk1L = xx[4];
      Double_t walk1R = xx[5];
      Double_t walk1B = xx[6];
      Double_t walk1T = xx[7];

      Double_t tLc = Events.t_L - walk0L/(1+walk1L*sqrt(Events.P_L));
      Double_t tRc = Events.t_R - walk0R/(1+walk1R*sqrt(Events.P_R));
      Double_t tBc = Events.t_B - walk0B/(1+walk1B*sqrt(Events.P_B));
      Double_t tTc = Events.t_T - walk0T/(1+walk1T*sqrt(Events.P_T));

      Double_t t0Hc = 0.5*(tLc+tRc);
      Double_t t0Vc = 0.5*(tBc+tTc);

       sigtcorH->Fill(t0Hc - Events.t0_H);
       sigtcorV->Fill(t0Vc - Events.t0_V);

      Double_t tRgen = Events.t_R + Events.tw_R;
      Double_t tLgen = Events.t_L + Events.tw_L;
      Double_t tTgen = Events.t_T + Events.tw_T;
      Double_t tBgen = Events.t_B + Events.tw_B;

      Double_t t0Hgen = 0.5*(tLgen+tRgen);
      Double_t t0Vgen = 0.5*(tBgen+tTgen);

       sigtgenH->Fill(t0Hgen - Events.t0_H);
       sigtgenV->Fill(t0Vgen - Events.t0_V);


       wL->Fill (Events.P_L,Events.t_L - 0.5*(Events.tcor_B+Events.tcor_T));
       wR->Fill (Events.P_R,Events.t_R - 0.5*(Events.tcor_B+Events.tcor_T));
       wT->Fill (Events.P_T,Events.t_T - 0.5*(Events.tcor_L+Events.tcor_R));
       wB->Fill (Events.P_B,Events.t_B - 0.5*(Events.tcor_L+Events.tcor_R));

      }
   gentof->Print();

   //
   TCanvas *c1 = new TCanvas("c1","c1 Tof_Calibration_Read Canvas",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);

   Float_t Eb=11;
   sprintf(string,"Eb=%.2f GeV\n",Eb);
   t1 = new TLatex(0.5,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   // t1->Draw();
   // t0->SetTitle(string);
   /*t0->GetXaxis()->SetRangeUser(xmin,xmax);
   t0->GetYaxis()->SetRangeUser(ymin,ymax);
   t0->GetXaxis()->SetTitleSize(0.05);
   t0->GetYaxis()->SetTitleSize(0.05);
   t0->GetYaxis()->SetTitleOffset(1.5);
   t0->GetXaxis()->SetTitle("Beam Energy (GeV)");
   t0->GetYaxis()->SetTitle("m_{g1g2}^{2}");*/
   t0->GetXaxis()->SetNdivisions(5);
   t0->SetMarkerColor(1);
   t0->SetMarkerStyle(21);
   t0->Draw();

   c1->cd(2); 
   c1_2->SetGridx();
   c1_2->SetGridy();
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);
   // c1_2->SetLogx();
   genxy->GetXaxis()->SetNdivisions(5);
   genxy->Draw();

   c1->cd(3);
   c1_3->SetGridx();
   c1_3->SetGridy();
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);
   dE1->GetXaxis()->SetNdivisions(5);
   dE1->Draw();

   c1->cd(4);
   c1_4->SetGridx();
   c1_4->SetGridy();
   c1_4->SetBorderMode(0);
   c1_4->SetFillColor(0);
   dE2->GetXaxis()->SetNdivisions(5);
   dE2->Draw();

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

   c2->cd(1);  
   c2_1->SetGridx();
   c2_1->SetGridy();
   c2_1->SetBorderMode(0);
   c2_1->SetFillColor(0);
   t0H->GetXaxis()->SetNdivisions(5);
   t0H->Draw();

   c2->cd(2);  
   c2_2->SetGridx();
   c2_2->SetGridy();
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);
   t0V->GetXaxis()->SetNdivisions(5);
   t0V->Draw();

   c2->cd(3);  
   c2_3->SetGridx();
   c2_3->SetGridy();
   c2_3->SetBorderMode(0);
   c2_3->SetFillColor(0);
   xraw->Draw("box");

   c2->cd(4);  
   c2_4->SetGridx();
   c2_4->SetGridy();
   c2_4->SetBorderMode(0);
   c2_4->SetFillColor(0);
   yraw->GetXaxis()->SetNdivisions(5);
   yraw->Draw("box");


   TCanvas *c3 = new TCanvas("c3","c3 Tof_Calibration_Read Canvas",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);

   t1 = new TLatex(0.5,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   // t1->Draw();
   // t0corH->SetTitle(string);
   t0corH->GetXaxis()->SetTitleSize(0.05);
   t0corH->GetYaxis()->SetTitleSize(0.05);
   t0corH->GetYaxis()->SetTitleOffset(1.5);
   // t0corH->GetXaxis()->SetTitle("Beam Energy (GeV)");
   // t0corH->GetYaxis()->SetTitle("m_{g1g2}^{2}");
   t0corH->GetXaxis()->SetNdivisions(5);
   // t0corH->SetMarkerColor(1);
   // t0corH->SetMarkerStyle(21);
   t0corH->Draw();

   c3->cd(2); 
   c3_2->SetGridx();
   c3_2->SetGridy();
   c3_2->SetBorderMode(0);
   c3_2->SetFillColor(0);
   t0corV->GetXaxis()->SetNdivisions(5);
   t0corV->Draw();

   c3->cd(3);
   c3_3->SetGridx();
   c3_3->SetGridy();
   c3_3->SetBorderMode(0);
   c3_3->SetFillColor(0);
   xcor->GetXaxis()->SetNdivisions(5);
   //xcor->Fit("gaus","","",0.4,0.8);
   xcor->Draw();

   c3->cd(4);
   c3_4->SetGridx();
   c3_4->SetGridy();
   c3_4->SetBorderMode(0);
   c3_4->SetFillColor(0);
   ycor->GetXaxis()->SetNdivisions(5);
   //ycor->Fit("gaus","","",0.1,0.5);
   ycor->Draw();


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

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

   // Define walk function

  Int_t npar=2;
  Double_t b0=2+6;
  Double_t b1=0.71;
  TF1 *walk_func = new TF1("walk_func",walk_func,vmin,vmax,npar);
   walk_func->SetParameters(b0,b1);
   
   t1 = new TLatex(0.5,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   // t1->Draw();
   // wL->SetTitle(string);
   wL->GetXaxis()->SetTitleSize(0.05);
   wL->GetYaxis()->SetTitleSize(0.05);
   wL->GetYaxis()->SetTitleOffset(1.5);
   wL->GetXaxis()->SetTitle("Pulse height (V)");
   // wL->GetYaxis()->SetTitle("Corrected time (ns)");
   wL->GetXaxis()->SetNdivisions(5);
   // wL->SetMarkerColor(1);
   // wL->SetMarkerStyle(21);
   wL->Fit(walk_func,"","",0.2,1.5);

   wL->Draw();

   c4->cd(2); 
   c4_2->SetGridx();
   c4_2->SetGridy();
   c4_2->SetBorderMode(0);
   c4_2->SetFillColor(0);
   wR->GetXaxis()->SetNdivisions(5);
  Double_t b0=2+6;
  Double_t b1=0.71;
   walk_func->SetParameters(b0,b1);
   wR->Fit(walk_func,"","",0.2,1.5);
   wR->Draw();

   c4->cd(3);
   c4_3->SetGridx();
   c4_3->SetGridy();
   c4_3->SetBorderMode(0);
   c4_3->SetFillColor(0);
   wT->GetXaxis()->SetNdivisions(5);
  Double_t b0=2+4;
  Double_t b1=0.71;
   walk_func->SetParameters(b0,b1);
   wT->Fit(walk_func,"","",0.3,1.5);
   wT->Draw();

   c4->cd(4);
   c4_4->SetGridx();
   c4_4->SetGridy();
   c4_4->SetBorderMode(0);
   c4_4->SetFillColor(0);
   wB->GetXaxis()->SetNdivisions(5);
  Double_t b0=2+3;
  Double_t b1=0.71;
   walk_func->SetParameters(b0,b1);
   wB->Fit(walk_func,"","",0.2,1.5);
   wB->Draw();

   TCanvas *c5 = new TCanvas("c5","c5 Tof_Calibration_Read Canvas",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);
   sigtrawH->GetXaxis()->SetNdivisions(5);
   sigtrawH->Draw();

   c5->cd(2); 
   c5_2->SetGridx();
   c5_2->SetGridy();
   c5_2->SetBorderMode(0);
   c5_2->SetFillColor(0);
   sigtrawV->GetXaxis()->SetNdivisions(5);
   sigtrawV->Draw();

   c5->cd(3);
   c5_3->SetGridx();
   c5_3->SetGridy();
   c5_3->SetBorderMode(0);
   c5_3->SetFillColor(0);
   sigtgenH->GetXaxis()->SetNdivisions(5);
   sigtgenH->Fit("gaus","","",20,30);
   sigtgenH->Draw();
   /*sigtcorH->GetXaxis()->SetNdivisions(5);
   sigtcorH->Fit("gaus","","",0.4,0.8);
   sigtcorH->Draw();*/

   c5->cd(4);
   c5_4->SetGridx();
   c5_4->SetGridy();
   c5_4->SetBorderMode(0);
   c5_4->SetFillColor(0);
   sigtgenV->GetXaxis()->SetNdivisions(5);
   sigtgenV->Fit("gaus","","",20,30);
   sigtgenV->Draw();
   /*sigtcorV->GetXaxis()->SetNdivisions(5);
   sigtcorV->Fit("gaus","","",0.1,0.5);
   sigtcorV->Draw();*/

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

   c6->cd(1);  
   c6_1->SetGridx();
   c6_1->SetGridy();
   c6_1->SetBorderMode(0);
   c6_1->SetFillColor(0);
   sigteasyH->GetXaxis()->SetNdivisions(5);
   sigteasyH->Fit("gaus","","",5,15);
   sigteasyH->Draw();

   c6->cd(2); 
   c6_2->SetGridx();
   c6_2->SetGridy();
   c6_2->SetBorderMode(0);
   c6_2->SetFillColor(0);
   sigteasyV->GetXaxis()->SetNdivisions(5);
   sigteasyV->Fit("gaus","","",5,15);
   sigteasyV->Draw();

   c6->cd(3);
   c6_3->SetGridx();
   c6_3->SetGridy();
   c6_3->SetBorderMode(0);
   c6_3->SetFillColor(0);
   sigtcorH->GetXaxis()->SetNdivisions(5);
   sigtcorH->Fit("gaus","","",10,20);
   sigtcorH->Draw();
   /*sigtgenH->GetXaxis()->SetNdivisions(5);
   sigtgenH->Fit("gaus","","",20,30);
   sigtgenH->Draw();*/

   c6->cd(4);
   c6_4->SetGridx();
   c6_4->SetGridy();
   c6_4->SetBorderMode(0);
   c6_4->SetFillColor(0);
   sigtcorV->GetXaxis()->SetNdivisions(5);
   sigtcorV->Fit("gaus","","",10,20);
   sigtcorV->Draw();
   /*sigtgenV->GetXaxis()->SetNdivisions(5);
   sigtgenV->Fit("gaus","","",20,30);
   sigtgenV->Draw();*/

      sprintf(filename,"tof_calibration_read_c1.eps");
      c1->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c1.png");
      c1->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c2.eps");
      c2->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c2.png");
      c2->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c3.eps");
      c3->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c3.png");
      c3->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c4.eps");
      c4->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c4.png");
      c4->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c5.eps");
      c5->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c5.png");
      c5->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c6.eps");
      c6->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c6.png");
      c6->SaveAs(filename);
      /*sprintf(filename,"tof_calibration_read_c7.eps");
      c7->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c7.png");
      c7->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c8.eps");
      c8->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c8.png");
      c8->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c9.eps");
      c9->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c9.png");
      c9->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c10.eps");
      c10->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c10.png");
      c10->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c11.eps");
      c11->SaveAs(filename);
      sprintf(filename,"tof_calibration_read_c11.png");
      c11->SaveAs(filename);*/

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

   Double_t b0=par[0];
   Double_t b1=par[1];
   Double_t x1=x[0];
   Double_t pi=3.14159;

   Int_t j;
   char string[256];

   Double_t correction = b0/(1+b1*sqrt(x1));
   // Double_t correction = b0/(1+b1*x1);

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

   return correction;
}