/* -------------------------------------------------------------------
   Monte Carlo Simulation for the Measurement of
   Charged Pion Life Time

   Detailed description can be found at:
    http://www1.gantep.edu.tr/~bingul/seminar/pion-lifetime/

  Ahmet Bingul, March 2006
  --------------------------------------------------------------------- */


#define N0 10000    // # of pion generated at t=0 
#define NP 10       // # of points to take data 
#define T0 2.602e-8 // mean life time (s) #define P0
#define P0 50.0     // pion beam momentum (MeV/c)
#define M0 139.5658 // pion rest mass (MeV/c^2)
#define c 3.0e+8    // speed of light
#include <stdlib.h>

 void generate(Double_t PIC[])
 //--------------------------------------------------
 // Generates N0 pions each has decay length DL and
 // Counts # of living pions from 1m to NP=10m.
 //--------------------------------------------------
 {
     Int_t i;
     Double_t Lambda,Gamma,Beta,DL,R;

     for(i=0;i<NP;i++) PIC[i] = 0;

     srand(time(NULL));

     for(i=1;i<=N0;i++)
     {
       Gamma = sqrt(P0*P0 + M0*M0)/M0;
       Beta  = sqrt(1.0-1.0/Gamma/Gamma);
       Lambda= Gamma*Beta*T0*c;

       R  = (double) rand()/RAND_MAX;
       DL = -Lambda*log(R);

       if(DL>=1.)  PIC[0]++;
       if(DL>=2.)  PIC[1]++;
       if(DL>=3.)  PIC[2]++;
       if(DL>=4.)  PIC[3]++;
       if(DL>=5.)  PIC[4]++;
       if(DL>=6.)  PIC[5]++;
       if(DL>=7.)  PIC[6]++;
       if(DL>=8.)  PIC[7]++;
       if(DL>=9.)  PIC[8]++;
       if(DL>=10.) PIC[9]++;
     }
 }


 void results(Double_t PIC[],Double_t p1, Double_t p2)
 //--------------------------------------------------
 // Prints the results to the user screen
 //--------------------------------------------------
 {
    Double_t T,Gamma,Beta;

    printf("==================================================\n");
    printf("\n");
    printf("*** Monte Carlo Simulation for the **********\n");
    printf("*** Measurement of Charged Pion Life Time ***\n");
    printf("\n");
    printf("---------------------------------------------\n");
    printf("Pion mom. generated (MeV/c): %10.3e\n",P0);
    printf("Pion life time      (s)    : %10.3e\n",T0);
    printf("Number of Pion generated   : %8d\n",N0);
    printf("Number of counter positon  : %8d\n",NP);
    printf("---------------------------------------------\n");
    printf("Values obtained for each position:\n");
    for(Int_t i=0;i<NP;i++)
        printf("%2d %8d %15.3e\n",i+1,PIC[i],log(PIC[i]));
    printf("---------------------------------------------\n");
    printf("Fitting results:\n");
    printf("Slope,     A = %15.4e\n",p2);
    printf("Intercept, B = %15.4e\n",p1);
    printf("---------------------------------------------\n");

    Gamma = sqrt(P0*P0 + M0*M0)/M0;
    Beta  = sqrt(1.0-1.0/Gamma/Gamma);
    T     = -1.0/(p2*Gamma*Beta*c);

    printf("\n");
    printf("Measured lifetime: %15.4e\n",T);
    printf("\n");
    printf("==================================================\n");
 }


void plt(){

   Int_t     i;
   Double_t  x[NP],PIC[NP];
   Double_t  p1,p2;

   TCanvas *c1 = new TCanvas("c1","Fitting Demo",10,10,800,500);
   TH1F    *hi = new TH1F("histo","Decay rate ve Distance",NP,0.5,10.5);
   lf          = new TF1("fitFun","[0]+x*[1]",1,10);

// generate data
   generate(PIC);

// Fill histogram
   for(i=0; i<NP; i++){
     x[i] = (i+1)*1.0;
     hi->SetBinContent(i+1,log(PIC[i]));
   }

// generate graph
   TGraph  *gr = new TGraph(NP,x,PIC);

// Settings
   lf->SetParName(0,"Intercept");
   lf->SetParName(1,"Slope");
   gr->SetLineColor(1);
   gr->SetLineWidth(2);
   gr->SetMarkerColor(1);
   gr->SetMarkerStyle(20);
   gr->SetMarkerSize(1.0);
   gr->SetTitle("Decay Rate vs Counter Distance");
   gr->GetXaxis()->SetTitle("x(m)");
   gr->GetYaxis()->SetTitle("N");
   gr->GetXaxis()->CenterTitle();
   gr->GetYaxis()->CenterTitle();
   gr->GetYaxis()->SetTitleOffset(2);
   hi->SetMarkerStyle(20);
   hi->SetMarkerSize(1.0);
   hi->SetStats(0);
   hi->SetMinimum(5);
   hi->SetMaximum(9.5);
   hi->SetXTitle("x(m)");
   hi->SetYTitle("log(N)");
   hi->GetXaxis()->CenterTitle();
   hi->GetYaxis()->CenterTitle();
   hi->GetYaxis()->SetTitleOffset(2);
   fitFun->SetLineColor(kBlue);

// divide canvas
   c1->Divide(2,1);

   c1->cd(1);
   gr->Draw("ACP");
   TLatex *t = new TLatex();
   t->SetNDC();
   t->SetTextFont(62);
   t->SetTextColor(1);
   t->SetTextSize(0.04);
   t->SetTextAlign(12);
   t->SetTextColor(1); t->DrawLatex(0.5,0.80,"N_{0} = 10,000");
   t->SetTextColor(1); t->DrawLatex(0.5,0.75,"P_{0} = 50 MeV/c");


   c1->cd(2);
// Perform fitting
   fitFun->SetParameters(9.0,-0.3);
   hi->Fit("fitFun","V+","p");

// get fit parameters
   p1 = fitFun->GetParameter(0);
   p2 = fitFun->GetParameter(1);

   results(PIC, p1, p2);

}