ctracking.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	ctrackingall
subroutine	ctracking
subroutine	cresettimer (aTrack)
subroutine	checkstat (str)
subroutine	cfixmag
subroutine	csampintl
subroutine	cfixproc
subroutine	ctruncpath
subroutine	cforcessp (atrack, icon)

Function/Subroutine Documentation

cfixmag()

subroutine cfixmag

(

)

Definition at line 369 of file ctracking.f.

References cgeomag(), clengsmallbc(), and ctransmagto().

Referenced by ctracking().

      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

Here is the call graph for this function:

Here is the caller graph for this function:

cfixproc()

subroutine cfixproc

(

)

Definition at line 492 of file ctracking.f.

References cthick2len(), d, d3, tointeract, and truncated.

Referenced by csampintl().

      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

Here is the call graph for this function:

Here is the caller graph for this function:

cforcessp()

subroutine cforcessp	(	type(track)	atrack,
		integer	icon
	)

Definition at line 619 of file ctracking.f.

References cfixmodel(), cpush(), kgnuc, knuc, masn, and masp.

Referenced by ctrackingall().

#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

Here is the call graph for this function:

Here is the caller graph for this function:

checkstat()

subroutine checkstat

(

character*(*)

str

)

Definition at line 327 of file ctracking.f.

Referenced by ctracking().

      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

Here is the caller graph for this function:

cresettimer()

subroutine cresettimer

(

type(track)

aTrack

)

Definition at line 302 of file ctracking.f.

References csetmintime().

Referenced by cbegin1ev().

      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

Here is the call graph for this function:

Here is the caller graph for this function:

csampintl()

subroutine csampintl

(

)

Definition at line 411 of file ctracking.f.

References cfixproc(), cfixtarget(), cfixtargetmuni(), cinismpintl(), cresetintinf(), csampeintl(), csampgintl(), csampnepintl(), cseecolpossible(), csetintinf(), kelec, kphoton, and true.

Referenced by ctracking().

      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

Here is the call graph for this function:

Here is the caller graph for this function:

ctracking()

subroutine ctracking

(

)

Definition at line 82 of file ctracking.f.

References borderh, borderl, ccomppathend(), cerrormsg(), cfixmag(), cfixmodel(), checkstat(), chookeabsorb(), chookeint(), cifxobssite(), cinteraction(), cobservation(), cpush(), cputtrinfo(), csampintl(), ctruncpath(), false, kgnuc, knuc, kpion, tobeobserved, tointeract, true, and truncated.

Referenced by ctrackingall().

      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

Here is the call graph for this function:

Here is the caller graph for this function:

ctrackingall()

subroutine ctrackingall

(

)

Definition at line 2 of file ctracking.f.

References cforcessp(), cifdead(), cobas(), cpop(), cscalerprod(), ctracking(), kelec, kphoton, and true.

Referenced by ceventloop().

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

Here is the call graph for this function:

Here is the caller graph for this function:

ctruncpath()

subroutine ctruncpath

(

)

Definition at line 575 of file ctracking.f.

References cmaxmovlen(), kkaon, kmuon, kpion, tointeract, and truncated.

Referenced by ctracking().

      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

Here is the call graph for this function:

Here is the caller graph for this function: