#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

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

void generate();      // .
void fitting();       // function prototypes
void results();       // .


int    i,PIC[NP];     // PIC: pion counter
double Lambda,Gamma,Beta,DL,R,A,B,T;


int main(void)
{

  srand( time(NULL) );

  generate();
  fitting();
  results();

  return 0;
}


 void generate()
 //--------------------------------------------------
 // Generates N0 pions each has decay length DL and
 // Counts # of living pions from 1m to NP=10m.
 //--------------------------------------------------
 {

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

     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 fitting()
 //--------------------------------------------------
 // Performs least square fitting for the data
 // generated from the decay rates vs distance plot.
 //--------------------------------------------------
 {
   double x[NP],y[NP];
   double sumx=0.0,sumy=0.0,sumxx=0.0,sumxy=0.0;

    for(i=0;i<NP;i++){
      x[i] = (double) i;
      y[i] = log(PIC[i]);

      sumx  += x[i];
      sumy  += y[i];
      sumxx += x[i]*x[i];
      sumxy += x[i]*y[i];
    }

    A = (NP*sumxy   - sumx*sumy )/(NP*sumxx - sumx*sumx); // slope
    B = (sumxx*sumy - sumx*sumxy)/(NP*sumxx - sumx*sumx); // intercept
 }


 void results()
 //--------------------------------------------------
 // Prints the results to the user screen
 //--------------------------------------------------
 {
    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(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",A);
    printf("Intercept, B = %15.4e\n",B);
    printf("---------------------------------------------\n");


    T = -1.0/(A*Gamma*Beta*c);

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