ctracking.f

Go to the documentation of this file.

      subroutine ctrackingall    
!     **************************************************
!       tracking particles given in the stack area.
!       During tracking, new particles may be produced
!       and pushed in the stack area.  They are all
!       treated here, until all particles are processed.
!     **************************************************
      implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include "Zincidentv.h"

      integer moreptcl, jold, icon
      real*8 smallAS/1./   ! electrons < 1GeV cannot make AS at all.
!c      real*8 u, iwgt  
      real*8 cosfromaxis
      external cblktracking
!//////
      logical show
      common /showshow/ show
!////////////
!
!      ***  do until no more ptcl in stack
      do while (.true.)
!           get one particle from stack
         call cpop(trackbefmove, moreptcl)
         if(moreptcl .eq. 0) goto 100   !  exit
         call rndsw(jold, 1)            !  specify 1st generator.
         call cifdead(trackbefmove, icon)  ! generator may be switched to 2
            
!                     icon=0; alive, 1 E<Emin, 2 long flight 3 h>H 4 h<L
!                             5 angle>limit
!            AS generation. 

         if(observeas .and. icon .le. 1 ) then
            if(trackbefmove%p%code .eq. kelec) then
               if(trackbefmove%asflag .eq. 0) then
                  if(trackbefmove%p%fm%p(4) .le. easwait) then
                     if(trackbefmove%p%fm%p(4) .gt. smallas) then ! 95/08/17
                        call cscalerprod(trackbefmove%vec%w,
     *                       dcatobsxyz, cosfromaxis)
!!!!!!!!!!!
                        write(0,*) ' cosfromaxis=',cosfromaxis
                        write(0, *) ' Track dc =',trackbefmove%vec%w
                        write(0,*) ' DcAtObsXyz=', dcatobsxyz
!!!!!!!!!!!!!!!                        
                        if(cosfromaxis .gt. 0.5 ) then
!                               as gneration for e with 60 deg or less from axis
                           call cobas(trackbefmove)
                        endif
                     endif
                     trackbefmove%asflag =1
                  endif
               endif
            elseif(trackbefmove%p%code .eq. kphoton) then
               if(icon .eq. 1) then
                  if(trackbefmove%asflag .eq. 0) then
                     if(trackbefmove%p%fm%p(4) .gt. smallas) then
                        icon = 0 ! follow until it becomes e-pair
                     endif
                  endif
               endif
            endif
         endif
!        ----------------------------
         if(icon .eq. 0) then
!              see if heavy and ssp is needed; icon =1 --> SSP so skip tracking here
            call cforcessp(trackbefmove, icon)
            if(icon .eq. 0) then
               call ctracking
            endif
         else
!             already dead ptcl has been considered in cifDead; nothing todo
         endif
      enddo
 100  continue
      end
! -------------------------------------tracking a ptcl
      subroutine ctracking
      use modxsecmedia, only: targetnucleonno, targetprotonno,
     * targetxs
      implicit none

#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include  "Zobs.h"
#include  "Zobsp.h"
#include  "Zobsv.h"
#include  "Zincidentv.h"
#include  "Zcode.h"
#include  "Zmanagerp.h"
#include  "Zstackv.h"
#include  "Zevhnv.h"
!                next is for DEBUG
!//////
      logical show
      common /showshow/ show
!/////////
#define DEBUG 0


      logical reset
      real*8 cosangle
      character*70 msg
      integer nextwhere, never, ptclidx
!c      integer icon

#if DEBUG == 1
      call checkstat('before fixModel')
#endif

      call cfixmodel( trackbefmove%p )      ! fix interaction model.
      call cfixmag            !   Mag field
!           sample interaction length.  The path may
!           be truncated to a shorter one than really sampled.
!           In that case, MoveStat == Truncated.
#if DEBUG == 1
      call checkstat('before csampIntl')
#endif
      call csampintl
#if DEBUG == 1
      call checkstat('after csampIntl')
#endif
!           Truncate the path if it is too long. Note, the truncated
!           path in the above process may be again truncated to
!           a shorter one. This happens, for example, if ptcl energy
!           is low.  Path inf is reset in IntInfArray.
      call ctruncpath

#if DEBUG == 1
      call checkstat('after ctruncPath')
#endif

!           compute path end inf. including scattering and
!           mag. deflection.  Path end inf is set in MovedTrack

      call ccomppathend

#if DEBUG == 1
      call checkstat('after ccompPathEnd')
#endif
!           see if MovedTrack crosses an observation site.
!         (or reaches BorderHeightL)  If so, reset MovedTrack. 
      call cifxobssite(nextwhere)


      if(trackbefmove%p%charge .ne. 0 .and. eabsorb(1) .ne. 0) then
!               4th arg. not used now
         call chookeabsorb(trackbefmove, movedtrack, energyloss, 0)
      endif
#if DEBUG == 1
      call checkstat('after cifXObsSite')
#endif
      if(trace .ne. 0) then
         if( mod(trace, 10) .le. 2) then
            call cputtrinfo     ! put trace info.
         elseif(mod(trace,10) .le. 4) then
            if( trackbefmove%pos%height .le. heightlist(1)) then
               call cputtrinfo  ! put trace info.
            endif
         endif
      endif
#if DETAILED_TRACKING >= 3
!            the user may kill the ptcl
      call cobservation(8)
#endif

!     
      if(movestat .eq. truncated) then
!        &&&&&&&&&&&&&& some may lose energy by dE/dx etc. 
!                       This should be recorded
         if(job .eq. 'newskel' .and. 
     *        (movedtrack%p%fm%p(4) - movedtrack%p%mass) .lt. 
     *        keminobs(1))  then
            never = -1
            call chookeint(never)
         endif
!        &&&&&&&&&&&&
!              stack current ptcl
         call cpush(movedtrack)
      elseif(movestat .eq. tointeract ) then
         if(zfirst%pos%depth .eq. 0.) then
!            if(  (MovedTrack.p.code .ge. kpion  .and. 
!     *            MovedTrack.p.code .le. knuc ) .or.
!     *         MovedTrack.p.code .eq. kgnuc ) then
!            from  v.7.10
            if(  (incidentcopy%p%code .ge. kpion  .and. 
     *            incidentcopy%p%code .le. knuc ) .or.
     *            incidentcopy%p%code .eq. kgnuc ) then

!
!                   from  v6.31 1ry is judged by Stack_pos
!                  (even if this is not used, almost no problem)
!                     ! if no stack, the particle should be 1ry (since Zfirst=0)

               if( stack_pos .eq. 0 .and.
     *             intinfarray(processno)%process .eq. 'coll') then
                  reset = .true.
               else
                  reset = .false.
               endif
            else
               reset = .true.
            endif
            if(reset) then
!             reset minimum time to reach the obs level (time from this
!             hight to the obs. level), if no mag. exists until
!                  the first collision point. or ptcl goes streight
!                  comment out from v6.31
!               call cresetTimer(MovedTrack)

               zfirst = movedtrack
               firstcola = targetnucleonno
               firstcolz = targetprotonno
               firstcolxs = targetxs
            endif
         endif
#if DETAILED_TRACKING >= 1 
         call cobservation(4)
#endif

#if DEBUG == 1
         call checkstat('before cinteraction')
#endif
         call cinteraction
#if DEBUG == 1
         call checkstat('after cinteraction')
#endif
      elseif(movestat .eq. tobeobserved ) then
         ptclidx = movedtrack%p%code  ! make 4 byte integer
         ptclidx = min(ptclidx, 8)   ! this is only for using KEminObs
         if(movedtrack%p%fm%p(4)-movedtrack%p%mass 
     *        .ge. keminobs(ptclidx) ) then    ! .gt. for <= uv6.00
!             call only for high energy particles
!                  (at AS generation, some ptcl may have < KEminObs
!                  or some may loose energy by dE/dx etc.
            call cobservation(2)

            if(job .eq. 'newskel' .and.
     *          endlevel .lt. noofsites .and. 
     *          movedtrack%where .eq. endlevel ) then
!                in this case, even E> KEminObs must be recorded
!                 for skeleton making for kahanshin
!                  where should be +1 at flesh time

               movedtrack%where = movedtrack%where + 1
               never = -3
               call chookeint(never)
               movedtrack%where = movedtrack%where - 1 ! for safty
            endif
         elseif(job .eq. 'newskel'  .and. 
     *          endlevel .eq. noofsites .and.
     *          movedtrack%where .lt. endlevel )  then
!            &&&&&&  this is to be recorded
!                where should be + 1 at flesh time
            movedtrack%where = movedtrack%where + 1
            never = -2
            call chookeint(never)
            movedtrack%where = movedtrack%where - 1
!            &&&&& 
         endif
!               update observation place
!         if(abs(ObsPlane) .eq. perpendicular) then
         movedtrack%where = nextwhere
!                  if more obs-site, stack current ptcl.
         if(.not. upgoing) then ! incident is downgoing
            if(movedtrack%where .gt. endlevel) then
                                ! no need to stack; discard ptcl.
            else
               movedtrack%where =max( movedtrack%where*1, 1)
               call cpush(movedtrack)
            endif
         else                   !incident is  Upgoing
            if(movedtrack%where .lt.1 ) then
                                ! no need to stack
            else
               movedtrack%where =
     *              min(movedtrack%where*1, endlevel)
               call cpush(movedtrack)
            endif
         endif
!        else
!            MovedTrack.where = nextwhere
!            call cpush(MovedTrack)
!         endif
      elseif(movestat .eq. borderl) then
         call cobservation(3)
         call cpush(movedtrack)
      elseif(movestat .eq. borderh) then
         call cobservation(1)
         call cpush(movedtrack)
      else
         write(msg, *) ' movestat=',movestat,' invalid'
         call cerrormsg(msg, 0)
      endif

      end
      subroutine cresettimer(aTrack)
      implicit none
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include  "Zobs.h"
#include  "Zobsp.h"
#include  "Zobsv.h"
#include  "Zincidentv.h"
#include  "Zcode.h"
#include  "Zmanagerp.h"
      type(track)::aTrack
!             reset minimum time to reach the obs level (time from this
!             hight to the obs. level), if no mag. exists until
!                  the first collision point. or ptcl goes streight
      if(mod(howgeomag, 2) .eq. 1 .or. 
     *     incidentcopy%p%charge .eq. 0 .and.
     *     obsplane .ne. 3) then
         call csetmintime(atrack)
      else
!         MovedTrack.t = 0.     ! bug bef uv6.09 ;;; but why commentout ??? 
      endif
      end
!           if you activate the next routine, Epics will have link problem  b

      subroutine checkstat(str)
      implicit none
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include  "Zobs.h"
#include  "Zobsv.h"
#include  "Zevhnv.h"

      logical deb
      common /cdebug/deb

      character*(*) str
      integer i
      integer nc
      data nc/0/

      nc = nc + 1
!      if(nc .gt. 750) then
!      if( deb) then
            write(0,*) ' ActiveModel=',activemdl
            write(0,*) str, ' code=',trackbefmove%p%code,
     *         ' sub',trackbefmove%p%subcode,
     *         ' chg',trackbefmove%p%charge
            write(0,*) 'Px,y,z e=', trackbefmove%p%fm%p(1:4)
            write(0,*)  'Move stat=',movestat
            write(0, *) ' where', trackbefmove%where
            write(0,*)  ' h==', trackbefmove%pos%height,
     *       ' d=', trackbefmove%pos%depth
            write(0,*) ' cos=',trackbefmove%vec%coszenith
            do i = 1, numberofinte
               write(0, *) ' process=',intinfarray(i)%process
               write(0, *) ' thickness=',intinfarray(i)%thickness
               write(0, *) ' length=',intinfarray(i)%length
            enddo
            write(0,*)'  ProcessNo=', processno
            write(0,*)  '--------------'
!         endif
!         endif
      end
!cccc  #endif
!-------------------------------------------------------
      subroutine cfixmag
      implicit none
#include "Ztrack.h"
#include "Ztrackv.h"
#include "Ztrackp.h"
#include "Zobs.h"
#include "Zobsp.h"
#include "Zobsv.h"

      real*8 x/0./, y/0./, z/0./
      save x, y, z
      real*8 distant
      integer icon
      real*8 u

      if(howgeomag .le. 2 .or. howgeomag .eq. 31) then
!             distant; change of B is < 1 % (for default
!             MagChgDist =20 km; you can change it)
         call clengsmallbc(trackbefmove, distant)
         if( (trackbefmove%pos%xyz%r(1) - x)**2 + 
     *        (trackbefmove%pos%xyz%r(2) - y)**2 + 
     *        (trackbefmove%pos%xyz%r(3) - z)**2 
     *        .gt. distant**2) then
!              if more than MagChgDist m from previous mag set.,
!               reset mag field.
!
            call cgeomag(yearofgeomag,  trackbefmove%pos%xyz, 
     *                mag, icon)

            call ctransmagto('xyz', trackbefmove%pos%xyz, 
     *        mag, mag)
            x = trackbefmove%pos%xyz%r(1)
            y = trackbefmove%pos%xyz%r(2)
            z = trackbefmove%pos%xyz%r(3)

         endif
      else
         mag = magfieldxyz
      endif
      end

!     --------------------------------------------------
      subroutine csampintl
      use modxsecmedia
!     **************************************************
!          
!          sample interaction length 
!          and interaction type.
!       Method.  Sample interaction lengths for all possible 
!                interactions (for collisions, in kg/m2, for decay
!                in m).  They are stored in a record /intinf/ IntInfArray
!                given in Ztrackv.h;  
!            In the routine,  cfixProc, 
!                we take the minimum of values given in kg/m2, and
!                convert it in real length (m).  In this process,
!                path truncation may occur (if the particle is
!                going upward, and there is very thin air there,
!                then the given thickness may not be realized.
!                Or if the length is too large and accuracy 
!                of convesion is not enough due to the earch
!                curverture)
!                If the decay process exists, we compare the length
!                given by the above treatment, and take shorter one.
!                If decay length is shorter,  we assume the decay
!                takes place, else some collision takes place if
!                the path is not truncated. In the latter case,
!                if the path is truncated, MoveStat == Truncated
!                is set.
      implicit none

#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"

      real*8 leng/1.d50/   ! will  be truncated  by ctruncPath
      integer:: jcon
!/////////
      logical::show
      common /showshow/ show
!//////////

!     **************************************************
      call cinismpintl   ! init for interaction length sampling
!
      if(reverse .ne. 0) then
         call csetintinf(leng, .true., 'none')
         processno = 1
      else
         if(trackbefmove%p%code .eq. kelec) then
            call csampeintl
         elseif(trackbefmove%p%code .eq. kphoton) then
            call csampgintl
         else
            call csampnepintl   ! non Electron Photon  Interaction.
         endif
         call cfixproc

         if( intinfarray(processno)%process == 'coll' ) then
!               some of the code can not deal with coll.of sigma 
!     etc.
            call cseecolpossible( trackbefmove%p, jcon)

            if( jcon == -1) then
!                  then force to decay
               call cresetintinf
            else
               call cfixtarget( media(1) )
            endif
         elseif( intinfarray(processno)%process == 'photop') then
            call  cfixtarget( media(1) )
         elseif( intinfarray(processno)%process == 'munuci') then
            call cfixtargetmuni(media(1))
         endif

         if(.not. freec .and. zfirst%pos%depth .eq. 0.) then
            intinfarray(processno)%length = 0.
            intinfarray(processno)%thickness = 0.
         endif
      endif

      end
!     **************************** cfixProc **********
      subroutine cfixproc
      implicit none

#include  "Zglobalc.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Zearth.h"


      real*8 h, leng, t, minlen, clen2thick, pcut
!     
      integer i, jcut
!     
       minlen = infty + infty
!              get vertical height
       h = trackbefmove%pos%height
       do i = 1, numberofinte
         if(.not. intinfarray(i)%decay) then
!              convert kg/m2 into length
            if(numberofinte .eq. 1 .or.
     *        intinfarray(i)%thickness .ne. infty) then
               call cthick2len(trackbefmove,
     *         intinfarray(i)%thickness, leng, t, jcut) 

               if( jcut .ne. 0) then
                  intinfarray(i)%thickness = t
               endif
!                  some work around for insufficient accuracy 
               if(leng .lt. 0.) then
                  leng = 1.d-3
                  intinfarray(i)%thickness = clen2thick(h, 
     *              trackbefmove%vec%coszenith, 
     *              intinfarray(i)%length )
               endif
            else
               leng = infty
               jcut = 0
            endif
         else
            leng = intinfarray(i)%length
            jcut = 0
         endif
         if(leng .lt. minlen) then
            if(jcut .ne. 0) then
               movestat = truncated
            else
               movestat = tointeract
            endif
            processno = i
            intinfarray(i)%length = leng
            minlen = leng
         endif
      enddo
!      if(leng .ge. minlen) then  !   bug in 6.00
!      if(leng .ge. Infty) then  
!         ProcessNo = 1
!         IntInfArray(1).length = 1.e5
!      endif         
      if(intinfarray(processno)%decay ) then
!         In the case  of muon, if individual knockon process  
!       is neglected (by parameter setting or with high Emin)
!       the length could be quite large (say, 6000  km).
!       and results in error in the next call.
!       To avoid that, we cut the path here
!
         if( trackbefmove%vec%coszenith .lt. 0.) then
            pcut = 300.d3
         else
            pcut = 30.d3
         endif
         if(intinfarray(processno)%length .gt. pcut) then
            movestat = truncated
            intinfarray(processno)%length  = pcut
         endif
         intinfarray(processno)%thickness = clen2thick(h, 
     *      trackbefmove%vec%coszenith, 
     *      intinfarray(processno)%length )
      endif
      end
!     ***********************
!          truncate a path of the particle, if the path
!        given in InfIntArray is too long, or the path traverses
!        an observation depth.
      subroutine ctruncpath
      implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"

      real*8 leng, thick
      real(8),parameter::veryshort=0.1 ! 10cm
      real(8),parameter::verylow=2d-3 !  2MeV
!
!//////////////////
      logical show
      common /showshow/ show
!////////////
      call cmaxmovlen(leng, thick)
      if(leng .lt. intinfarray(processno)%length ) then
         movestat = truncated
         intinfarray(processno)%length = leng
         intinfarray(processno)%thickness  = thick
         if(leng < veryshort ) then
!            some very unhappy case, B and dE/dx forces
!             muon (typically) path very short and the
!             muon cannot decay, while E loss is very small
!             and muon KE is still > 0 after the path.
!             This is repeated and KE becomes smaller and
!             smaller (say, until 10^-13 GeV); it may need
!             huge steps for traveling only a few cm and   
!             cpu time more than 5 sec.  This is avoided
!             by next 
            if( movedtrack%p%code == kmuon .or.
     *          movedtrack%p%code == kpion .or.
     *          movedtrack%p%code == kkaon ) then
               if( movedtrack%p%fm%p(4)- movedtrack%p%mass 
     *            < verylow ) then
!                  force to decay; 
                  movestat = tointeract
               endif
            endif
         endif
      endif
      end
!      *************** for same model, we must use simple super position for heavy
!      If so, current particle is decomposed into nucelons and put into stack
!      and icon =1  is made. else icon = 0
      subroutine  cforcessp(atrack, icon)
#include "Ztrack.h"
#include "Zevhnv.h"
#include "Zevhnp.h"
#include "Zcode.h"
#include "Zmass.h"

      integer icon
      type(track)::atrack  ! input current projectile.

      type(track)::pj
      integer i, j
      real*8  p, po, pr

      icon = 0
      if(atrack%p%code .eq. kgnuc) then
         call cfixmodel(atrack%p)
         if(activemdl .eq. 'dpmjet3' .and. 
     *      atrack%p%subcode .gt. 18) then
!/////            icon = 1 ///// 
         elseif( activemdl .eq. 'incdpm3') then
            icon = 1
         endif
      endif
      if( icon  .eq. 1 ) then
         pj = atrack
         pj%p%fm%p(4) = atrack%p%fm%p(4)/atrack%p%subcode
         pj%p%mass = masp
         pj%p%code = knuc
         pj%p%subcode = -1
         pj%p%charge = 1

         p = sqrt(pj%p%fm%p(4)**2 - masp**2)      
!          next one leads to seg. violation on opteron and ifc64
!         po = atrack.p.fm.p(1)**2 + atrack.p.fm.p(2)**2 +
!     *           atrack.p.fm.p(3)**2 
!       (  still by version 9. )

         po = atrack%p%fm%p(1)**2 + atrack%p%fm%p(2)**2
         po = po +  atrack%p%fm%p(3)**2 

         po = sqrt(po)
         pr = p/po
         do j = 1, 3
            pj%p%fm%p(j) = atrack%p%fm%p(j)*pr
         enddo
         do i = 1, atrack%p%charge
            call cpush(pj)
         enddo

         pj%p%mass = masn
         p = sqrt(pj%p%fm%p(4)**2 - masn**2)      
         pr = p/po
         pj%p%charge = 0
         do j = 1, 3
            pj%p%fm%p(j) = atrack%p%fm%p(j)*pr
         enddo
         do i = 1, atrack%p%subcode - atrack%p%charge
            call cpush(pj)
         enddo
      endif
      end