csampPrimary.f

Go to the documentation of this file.

!       *******************************************************
!       *  csampPrimary: samples a primary particle; its type
!       *                 and energy
!       *******************************************************
!   Note:  Here primary spectrum information is contained in common
!       variable Prim (done by init routine). 
!       If we  give it to csampPrimary0 we get
!       one sampled priamary of which energy etc is given in
!       variable in Prim.  Only /ptcl/ information is returned
!       to the caller.
!
        subroutine csampprimary(p, fin)
!          p: /ptcl/  output.  energy, particle code, subcode,
!                              charge, mass are set.
!  

        implicit none
#include  "Zmanagerp.h"
#include  "Zptcl.h"
#include  "Zprimary.h"
#include  "Zprimaryv.h"

        type(ptcl)::  p
        integer fin ! output. if no more sim. -->1 else 0


        fin = 0
        if( desteventno(2) .lt. 0) then
           if(abs(desteventno(2)) .le.  prim%NoOfSamplings) then
              fin = 1
           endif
        endif
        if(fin .eq. 0) then
           call csampprm0(prim)


!          call cconv_prim_e(Prim)   ! to total energy
           call cconv_prim_e(prim%each(prim%label), prim%sampled_e, 
     *      prim%particle)         
           p = prim%particle
           prim%NoOfSamplings =  prim%NoOfSamplings  + 1 ! update counter.
           prim%NoOfSampComp(prim%label, 1) =
     *          prim%NoOfSampComp(prim%label, 1) +1
        endif
        end
!       **********************   
        subroutine cupdtprimc
        implicit none
#include  "Zptcl.h"
#include  "Zprimary.h"
#include  "Zprimaryv.h"

        prim%NoOfSampComp(prim%label, 2) =
     *       prim%NoOfSampComp(prim%label, 2) +1
        end

!       ************************************
        subroutine csampprm0(prm)
        implicit none

#include  "Zptcl.h"
#include  "Zprimary.h"
        type(primaries):: prm

        real*8 p_or_e
!
        integer i
        real*8 u
!
!          select one component
        call rndc(u)

        i = 1
        do while (u .gt. prm%cummInteFlux(i) )
           i = i + 1
        enddo
        prm%label = prm%each(i)%label
!
        call csampprm1(prm%each(prm%label), p_or_e)


        prm%sampled_e = p_or_e
        end
        subroutine csampprm1(each, p_or_e)
        implicit none

#include  "Zptcl.h"
#include  "Zprimary.h"
        type(component):: each
        real*8 p_or_e
!
        real*8 e1temp, ombeta, u
        integer i
!
        if( each%no_of_seg .eq. 0) then
           p_or_e = each%emin       ! = emax
        else
           call rndc(u)
           i = each%no_of_seg 
           do while (u .gt. each%norm_inte(i))
               i = i - 1
           enddo
!             use i-th segment
           if(each%diff_or_inte .eq. 'd') then    
               ombeta = (1.d0 - each%beta(i))           
               call rndc(u)
               if(abs(ombeta) .gt. 1.d-6) then
                  e1temp = each%energy(i)**ombeta
                  p_or_e =( u* (each%energy(i+1)**ombeta - e1temp) +
     *             e1temp )** (1.d0/ombeta)
               else    
                  p_or_e =( each%energy(i)/each%energy(i+1) )**u
     *                * each%energy(i+1)                  
               endif
           elseif(each%diff_or_inte .eq. 'i') then
               call rndc(u)
               p_or_e = each%energy(i)* (1.d0 - 
     *          u*(1.d0 - each%norm_inte(i+1)/each%norm_inte(i)))
     *          **(-1.d0/each%beta(i))
           else
               write(*, *) ' invlid diff_or_inte=',each%diff_or_inte,
     *         ' for primary=',each%symb
               stop 9999
           endif     
         endif  
         end
        subroutine cconv_prim_e(comp, e_or_p,  aPtcl)
        implicit none
!          given energy or momentum as given in primary 
!          specification is converted to total enregy.

#include  "Zptcl.h"
#include  "Zcode.h"
#include  "Zprimary.h"
         type(component):: comp  ! input. primary component
         real*8  e_or_p           ! input. energy or momentum as
                                  !        specified in primary data file
         type(ptcl)::aPtcl       ! output. 'partcle' of the primary 
         integer  code, massn
!       
         real*8 inGeV  ! p_or_e in  Gev
         real(8),parameter::hbarc=197.0e6*1e-15/1.e-9 ! hbarc in eV nm
         real(8),parameter::hc=hbarc*2*3.14159  ! hc/wl = eV
!
         ingev= e_or_p *   comp%togev
         call cmkptc(comp%code,
     *               comp%subcode,
     *               comp%charge,
     *               aptcl)        

         if(comp%code .eq. kmuon) then
!            give polarization for muon primary
            call csetmuonpol(0.d0)
         endif

         if(comp%etype .eq. 'e/p' .or.
     *      comp%etype .eq. 'e') then
            aptcl%fm%p(4) = ingev
     *                  
         elseif(comp%etype .eq.'p/p' .or.
     *      comp%etype .eq. 'p')then
            aptcl%fm%p(4) = sqrt(aptcl%mass**2 +
     *             ingev**2)
         elseif(comp%etype .eq. 'ke/p' .or.
     *      comp%etype .eq. 'ke') then
            aptcl%fm%p(4) = ingev + aptcl%mass
         elseif(comp%etype .eq. 'e/n') then
            code = comp%code
            if(code .eq. kgnuc) then
               aptcl%fm%p(4) = aptcl%subcode * ingev
            elseif(code .ge. kalfa .and. 
     *         code .le. khvymax) then
               call cghvm(code, massn)
               aptcl%fm%p(4) = massn * ingev
            else
               aptcl%fm%p(4) =  ingev
            endif
         elseif(comp%etype .eq. 'ke/n') then
            code = comp%code
            if(code .eq. kgnuc) then
               aptcl%fm%p(4) = aptcl%subcode * ingev +
     *            aptcl%mass
            elseif(code .ge. kalfa .and. 
     *         code .le. khvymax) then
               call cghvm(code, massn)
               aptcl%fm%p(4) = massn * ingev +
     *            aptcl%mass
            else
               aptcl%fm%p(4) = ingev +
     *            aptcl%mass
            endif
     elseif(comp%etype .eq. 'p/n') then
            code = comp%code
            if(code .eq. kgnuc) then
               aptcl%fm%p(4) =
     *            sqrt( (aptcl%subcode * ingev) **2 + 
     *            aptcl%mass **2)
            elseif(code .ge. kalfa .and. 
     *         code .le. khvymax) then
               call cghvm(code, massn)
               aptcl%fm%p(4) =
     *            sqrt( (massn * ingev) **2 + 
     *            aptcl%mass **2)
            else
               aptcl%fm%p(4) = sqrt(ingev**2 +
     *            aptcl%mass**2)
            endif
         elseif(comp%etype .eq. 'rig') then
            aptcl%fm%p(4) =sqrt(
     *          ( ingev*aptcl%charge)**2 + aptcl%mass**2 )
         elseif(comp%etype == 'nm' ) then
            if(comp%code /= klight ) then
               write(0,*) ' nm unit can be used only for light'
               stop
            endif
            aptcl%fm%p(4) =  hc/e_or_p  ! hc/wl --> eV 
         else
            write(*, *) ' energy type=', comp%etype,
     *                  ' invalid. label=', comp%label, ' symb =',
     *                  comp%symb
            stop 9999
         endif
         if(comp%code == klight .and. comp%etype /= 'nm' ) then
            !  energy should be in eV so GeV to eV  conversion
            aptcl%fm%p(4) = aptcl%fm%p(4)*1.d9
         endif
        end
!       *******************************************************
        subroutine cconv_prim_e2(comp, E, e_or_p)
        implicit none
!         inverse of cconv_prim_e

#include  "Zptcl.h"
#include  "Zcode.h"
#include  "Zprimary.h"
         type(component):: comp  ! input. primary component
         real*8  E                ! input. total energy of a particle
                                  !   of 'comp' primary
         real*8 e_or_p            ! output. e or p (/n, etc) as given 
                                  !      in primary paticle file.
         integer  code, massn
!       
         real*8 Ex         !  e_or_p in GeV

         type(ptcl)::aPtcl

         call cmkptc(comp%code,
     *               comp%subcode,
     *               comp%charge,
     *               aptcl)        


         if(comp%etype .eq. 'e/p' .or.
     *      comp%etype .eq. 'e') then
            ex = e
         elseif(comp%etype .eq.'p/p' .or.
     *      comp%etype .eq. 'p')then
            ex = sqrt(e**2 - aptcl%mass**2)

         elseif(comp%etype .eq. 'ke/p' .or.
     *      comp%etype .eq. 'ke') then
            ex = e - aptcl%mass

         elseif(comp%etype .eq. 'e/n') then
            code = comp%code
            if(code .eq. kgnuc) then
               ex = e/aptcl%subcode
            elseif(code .ge. kalfa .and. 
     *         code .le. khvymax) then
               call cghvm(code, massn)
               ex = e/massn
            else
               ex = e
            endif
         elseif(comp%etype .eq. 'ke/n') then
            code = comp%code
            if( code .eq. kgnuc ) then
               ex =( e -   aptcl%mass)/aptcl%subcode
            elseif(code .ge. kalfa .and. 
     *         code .le. khvymax) then
               call cghvm(code, massn)
               ex =( e -   aptcl%mass)/massn
            else
               ex = e - aptcl%mass
            endif
     elseif(comp%etype .eq. 'p/n') then
            code = comp%code
            if( code .eq. kgnuc) then
               ex =sqrt( e**2 -   aptcl%mass **2)/aptcl%subcode
            elseif(code .ge. kalfa .and. 
     *         code .le. khvymax) then
               call cghvm(code, massn)
               ex =sqrt( e**2 -   aptcl%mass **2)/massn
            else
               ex = sqrt(e**2 - aptcl%mass**2)
            endif
         elseif(comp%etype .eq. 'rig') then
            ex = sqrt( e**2 - aptcl%mass**2) /aptcl%charge
         else
            write(*, *) ' energy type=', comp%etype,
     *                  ' invalid. label=', comp%label, ' symb =',
     *                  comp%symb
            stop 9999
         endif
        e_or_p = ex/ comp%togev
        end
!       *******************************************************
!       *  cqPrimE: inquire sampled primary energy or p or rigidity
!       *                    as it is
!       *******************************************************
!
        subroutine cqprime(p_or_e)
!  
        implicit none

#include  "Zptcl.h"
#include  "Zprimary.h"
#include  "Zprimaryv.h"
        real*8 p_or_e
        call cqprime0(prim, p_or_e)
        end
!       *******************************************************
!       *  cqPrimE0: inquire sampled primary energy or p or rigidity
!       *                     as it is
!       *******************************************************
!
        subroutine cqprime0(prm, p_or_e)
!  
        implicit none

#include  "Zptcl.h"
#include  "Zprimary.h"
        real*8 p_or_e
        type(primaries):: prm
!
        p_or_e = prm%sampled_e
!
        end
!       *******************************************************
!       *  cqPrimLabel: inquire sampled primary label

!       *******************************************************
!
        subroutine cqprimlabel(label)
!  
        implicit none

#include  "Zptcl.h"
#include  "Zprimary.h"
#include  "Zprimaryv.h"
        integer label
        call cqprimlabel0(prim, label)
        end
!       ************************  inquire the number of primaries sampled
        subroutine cqnoofprim(no)
!       ************************
        implicit none
#include  "Zmanagerp.h"
        integer no   ! output.  no. of sampled primaries so far.
        no = eventno
        end

!       *******************************************************
!
        subroutine cqprimlabel0(prm,label)
!  
        implicit none
#include  "Zptcl.h"
#include  "Zprimary.h"
        integer label
        type(primaries):: prm

        label = prm%label
        end
!       *******************************************************
!           inquire all about current primaries
        subroutine cqprimary(prm)
!         prm /primaires/ output.
!  
        implicit none
#include  "Zptcl.h"
#include  "Zprimary.h"
#include  "Zprimaryv.h"

        type(primaries):: prm
        prm = prim
        end