chookASbyH2.f

Go to the documentation of this file.

!
!     chookASbyH2:  This is a special user hook to compute
!                  AS generated by heavy primaries with
!                  generate='qas'.  This is a slave of chookASbyH.
!
      subroutine chookbgrunas
      implicit none
#include "Zmanagerp.h"
      include "ZASbyH.h"  ! this is in UserHook
      integer ios, i,j


!      open component disk file,
!      We assume that  over all information is in the first 
!      and 2nd record of the direct access file.
!      and record number =3 to NoOfCompS + 2  have component showers.
!
!   *************** record organization *****************
!
!    First rec:  NoOfCompS, NdpethComp, E0
!    Second rec: Depths
!    later rec:  first collision depth (vertical) and sizes, ages
!
!   **************************************************** 


!        
!        
!        
      open(13, file='cs', status='OLD')

!        read direct access file.
!        Number of  compoentent showers
!        Number of  depths
!        Depth list

      read(13, *) noofcomps, ndepthcomp, e0comp
     *      , firstcoldep
      if(noofcomps .gt. ncs) then
         write(0, *) 'too many comp showers '
         stop
      endif
      read(13, *)
     *      (depthcomp(i), i = 1, ndepthcomp)     

    do i=1,noofcomps
      read(13,*) (elecsizecomp(i,j),j=1,ndepthcomp)
      read(13,*) (elecagecomp(i,j),j=1,ndepthcomp)
    end do

    close(13)

      end

      subroutine chooknepinta
      implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Zobs.h"
#include  "Zobsv.h"
      include "ZASbyH.h"
!  #include  "Ztrackp.h"
      
      

!   MovedTrack is the particle that made interaction
!   Pwork contains produced particles.
!   Nproduced has the number of particles in Pwork
!   IntInfArray(ProcessNo) contains the type of interaction
!
!        IntInfArray(ProcessNo).process  will have
!             'col' or 'decay'

      integer i, recn, ios, j
      integer noofintn             ! no. of interacting nucleons
      type(ptcl):: ptcla(maxheavymassn)  ! to get nucleon inf.
      real*8 u     ! uniform random number
      real*8 esizetemp(ndepthcompm), depstep, depfirst, dep
      real*8 eagetemp(ndepthcompm)
      real*8 ans, zobs, t, cosz, length, clenbetween2h, eno
      real*8  clen2thick, age
    

!        inquire the interacting nucleons
      call cqinteptcl(ptcla, noofintn)


!          for each nucleon
      do i = 1, noofintn
!           select a component shower
         call rnde(u, 1)   ! get one unif. random. n.
!            
         recn = u* noofcomps + 1

!         you may need to verify that
!         the primary energy of the component shower and 
!          ptclA.fm.p(4) is almost the same for safety.

         depfirst = depthcomp(1)   !  first depth
         depstep = depthcomp(2) - depthcomp(1)  ! equi step depth


!          move to real*8  array for use of kintp3
         do j =1, ndepthcomp
            esizetemp(j) = elecsizecomp(recn,j)    ! log10( size )
        eagetemp(j)  = elecagecomp(recn, j)

         enddo
!              cos of c.s
         cosz = movedtrack.vec.coszenith
!          compute size at  
         do j = 1, noofassites
            zobs = asobssites(j).pos.depth 
!             current nucleon depth to obseration depth
!             (vertical depth)
            dep = zobs - movedtrack.pos.depth 
            if(dep .gt. 0.) then
                if(abs(cosz) .gt. 0.5) then
!                         amount of air
                   t= dep/abs(cosz)
                else
!                  slant lenngth between two depths
                   length = clenbetween2h(
     *                  movedtrack.pos.radiallen,
     *                  asobssites(j).pos.radiallen,
     *                  cosz)
!                         amount of air
                   t = clen2thick(
     *                movedtrack.pos.height,
     *                cosz,  length)
                endif
!                  current obs depth corresponds to dep in c.s
                dep = t + firstcoldep
!                 interpolate; use diff. one if DepthComp is not
!                  equisteped.
                call kintp3(esizetemp, 1, ndepthcomp,  depfirst,
     *              depstep, dep, ans)


!             *********************************
!                   for age
                call kintp3(eagetemp, 1, ndepthcomp,  depfirst,
     *             depstep, dep, age)
!             *********************************

!                   IMPORTANT:   assume that eSizeTemp is in log10
                eno =    10.**ans * movedtrack.wgt

                asobssites(j).esize =  asobssites(j).esize + eno

!             *********************************************
               asobssites(j).age = asobssites(j).age +
     *             eno * age
!             *********************************************
            endif
         enddo

      enddo
      end
!**********

!                   IMPORTANT:   assume that eSizeTemp is in log10
                eno =    10.**ans * movedtrack.wgt
                asobssites(j).esize =  asobssites(j).esize + eno

!             *********************************************
!               ASObsSites(j).age = ASObsSites(j).age +
!     *             eno * age
!             *********************************************
            endif
         enddo
      enddo
      end