PROGRAM PION_LIFETIME
!---------------------------------------------------------------------
!  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
!---------------------------------------------------------------------
IMPLICIT NONE
INTEGER, PARAMETER :: N0 = 10000      ! # of pion generated at t=0
INTEGER, PARAMETER :: NP = 10         ! # of points to take data
REAL,    PARAMETER :: T0 = 2.602E-8   ! mean life time (s)
REAL,    PARAMETER :: P0 = 50.0       ! pion beam momentum (MeV/c)
REAL,    PARAMETER :: M0 = 139.5658   ! pion rest mass (MeV/c^2)
REAL,    PARAMETER :: c  = 3.0E+8     ! speed of light
INTEGER :: I,PIC(NP)                  ! PIC: pion counter
REAL    :: Lambda,Gamma,Beta,DL,R,A,B,T


CALL RANDOM_SEED
CALL Generate
CALL Fitting
CALL Results

CONTAINS


   SUBROUTINE Results
   !--------------------------------------------------
   ! Prints the results to the user screen
   !--------------------------------------------------
    WRITE(*,*)
    WRITE(*,12) "*** Monte Carlo Simulation for the **********"
    WRITE(*,12) "*** Measurement of Charged Pion Life Time ***"
    WRITE(*,*)
    WRITE(*,12) "---------------------------------------------"
    WRITE(*,10) "Pion mom. generated (MeV/c): ",P0
    WRITE(*,10) "Pion life time      (s)    : ",T0
    WRITE(*,11) "Number of Pion generated   : ",N0
    WRITE(*,11) "Number of counter positon  : ",NP
    WRITE(*,12) "---------------------------------------------"
    WRITE(*,12) "Values obtained for each position:"
    DO I=1,NP
      PRINT 13,I,PIC(I),LOG(PIC(I)*1.0)
    END DO
    WRITE(*,12) "---------------------------------------------"
    WRITE(*,12) "Fitting results:"
    WRITE(*,14)  "Slope,     A = ",A
    WRITE(*,14)  "Intercept, B = ",B
    WRITE(*,12) "---------------------------------------------"

    t = -1.0/(A*Gamma*Beta*c)

    WRITE(*,*)
    WRITE(*,14) "Measured lifetime: ",t
    WRITE(*,*)

    10 FORMAT(A,ES10.3)
    11 FORMAT(A,I5)
    12 FORMAT(A)
    13 FORMAT(I2,I8,ES15.3)
    14 FORMAT(A,ES15.4)

   END SUBROUTINE Results

   SUBROUTINE Generate
   !--------------------------------------------------
   ! Generates N0 pions each has decay length DL and
   ! Counts # of living pions from 1m to NP=10m.
   !--------------------------------------------------
     PIC = 0
     DO I=1,N0
       Gamma = SQRT(P0**2+M0**2)/M0
       Beta  = SQRT(1.0-1.0/Gamma**2)
       Lambda= Gamma*Beta*T0*c

       CALL RANDOM_NUMBER(R)
       DL    =-Lambda*LOG(R)

       IF(DL>=1.)  PIC(1) = PIC(1)  + 1
       IF(DL>=2.)  PIC(2) = PIC(2)  + 1
       IF(DL>=3.)  PIC(3) = PIC(3)  + 1
       IF(DL>=4.)  PIC(4) = PIC(4)  + 1
       IF(DL>=5.)  PIC(5) = PIC(5)  + 1
       IF(DL>=6.)  PIC(6) = PIC(6)  + 1
       IF(DL>=7.)  PIC(7) = PIC(7)  + 1
       IF(DL>=8.)  PIC(8) = PIC(8)  + 1
       IF(DL>=9.)  PIC(9) = PIC(9)  + 1
       IF(DL>=10.) PIC(10)= PIC(10) + 1
     END DO
   END SUBROUTINE Generate

   SUBROUTINE Fitting
   !--------------------------------------------------
   ! Performs least square fitting for the data
   ! generated from the decay rates vs distance plot.
   !--------------------------------------------------
   REAL :: x(NP),y(NP)

    x = (/ (I*1.0, I=1,NP) /)
    y = LOG(REAL(PIC))

    A = (NP*SUM(x*y) - SUM(x)*SUM(y)) / &         ! slope
        (NP*SUM(x*x) - SUM(x)*SUM(x))
    B = (SUM(x*x)*SUM(y) - SUM(x)*SUM(x*y)) / &   ! intercept
        (NP*SUM(x*x)     - SUM(x)*SUM(x))
   END SUBROUTINE Fitting


END PROGRAM