ccompPathEnd.f

Go to the documentation of this file.

!       ccompPathEnd
!       compute path end information. including multiple scattering and mag.
!       deflection
!
      subroutine ccomppathend
      implicit none

#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
!
!
!         get endpoint coord. assuming no scat, no mag. effect. at first
!         MovedTrack has complete information at path end.
       movedtrack = trackbefmove   !   copy TrackBefMove into MovedTrack. 
       if(intinfarray(processno)%length .gt. 0.0 ) then
          call cmovestreight(intinfarray(processno)%length, 
     *             trackbefmove%vec%w)
!         If a chargde ptcl, compute energy loss by dE/dx  and  reset energy
!         in MovedTrack. If it is too large, path is truncated. Also compute
!         deflection.

          if(trackbefmove%p%charge .ne. 0) then
!            if(.not. (mod(HowGeomag, 10) .eq. 1 .and.     this existed from uv4.92 to 5.13
!     *                Zfirst .eq. 0.)) then
             call cputenergyloss
!               endif
!                if "if" is omitted, segmentation violation
!                may take place occasionally
             if(movestat .ne. dead) then
                call cputdeflection
             endif
          endif
       else
          energyloss = 0.
       endif
       end
!     --------------------------------------------------------
       subroutine cmovestreight(leng, dir)
!        move a track by leng m
!         after this is called,  MovedTrack has inf for moved pos.
!         MovedTrack.pos.where may not be correct. It must be
!         corrected by seeing if the track passes acrosss a
!         observation layer.
       implicit none

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

       real*8  leng, cnewcos
       type(coord)::dir        ! input.
!
!            new coord in xyz system.

       movedtrack%pos%xyz%r(1) = trackbefmove%pos%xyz%r(1)
!     *                      + TrackBefMove.vec.w.r(1) * leng
     *                      + dir%r(1) * leng
       movedtrack%pos%xyz%r(2) = trackbefmove%pos%xyz%r(2)
!     *                      + TrackBefMove.vec.w.r(2) * leng
     *      + dir%r(2) * leng
       movedtrack%pos%xyz%r(3) = trackbefmove%pos%xyz%r(3)
!     *                      + TrackBefMove.vec.w.r(3) * leng
     *                      + dir%r(3) * leng
       call csetpos(movedtrack%pos) ! set height, radial length, thickness
!            get new cos
       movedtrack%vec%coszenith = cnewcos(movedtrack%pos%radiallen,
     *  trackbefmove%vec%coszenith, leng)
       if(timestructure ) then
          call cgetbeta(movedtrack%p, beta)  
!             note righthand is not MovedTrack.t
          if(beta .gt. 0.) then
             movedtrack%t = trackbefmove%t + leng/beta
          endif
       endif
       end
!      -----------------------------------
       subroutine cputenergyloss
!         observation layer.
       implicit none

#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
 
!     
       real*8  rho, cvh2den, dedt, dedtF,  thick, leng, cupsilon
       save dedt, dedtf
       real*8  csyncTELoss
       integer jcut
       real*8  dedx, dedxF  ! output.  dedt, dedtF is put


       real*8  dedxmu, dedxmuF

!
!          first consider the synchrotron loss.
!
       if( trackbefmove%p%fm%p(4) .gt. magbrememin) then
         if( magbrem .eq. 1 .and. trackbefmove%p%code .eq. kelec) then
        upsilon = cupsilon(trackbefmove%p, mag)
        if(upsilon .gt. upsilonmin) then
!                 compute energy loss due to sychrotron rad.
                 energyloss = csyncteloss(upsilon) * 
     *           intinfarray(processno)%length
                 if(reverse .eq. 0) then
                    movedtrack%p%fm%p(4) = 
     *                movedtrack%p%fm%p(4) -  energyloss
                 elseif(reverse .eq. 2) then
                    movedtrack%p%fm%p(4) = 
     *                movedtrack%p%fm%p(4) +  energyloss
                 endif
!                 don't  worry about death
        endif
         endif
       endif 
!            
!           dE/dX
!
       rho = cvh2den(trackbefmove%pos%height)
       call cdedxinair(trackbefmove%p, rho, dedt, dedtf ) ! dedt; GeV/(kg/m2)

       if(trackbefmove%p%code .eq. kmuon ) then
!           dE/dx due to muon pair, brem, nuc.i
          call cmudedx(muni, mubr, mupr, trackbefmove%p%fm%p(4),
     *         dedxmu, dedxmuf)   ! dedxmu in GeV/(g/cm2)
          dedxmu = dedxmu /10.  !  GeV/(kg/m2)
          dedxmuf = dedxmuf /10.  !  GeV/(kg/m2)
          dedt = dedt + dedxmu
          dedtf = dedtf + dedxmuf
       endif

       energyloss =  dedt * intinfarray(processno)%thickness
       if(reverse .eq. 0) then
           movedtrack%p%fm%p(4) = movedtrack%p%fm%p(4) -  energyloss          
!           see if <=mass
           if(movedtrack%p%fm%p(4) .lt. movedtrack%p%mass) then
              energyloss =( trackbefmove%p%fm%p(4) - movedtrack%p%mass )
              if(dedt .gt. 0.) then
                 thick =energyloss / dedt
              else
                 thick = 0.
              endif
              movedtrack%p%fm%p(4) = movedtrack%p%mass

!              call cthick2len(TrackBefMove.pos.height, 
!     *          TrackBefMove.vec.coszenith, thick, leng, thick, jcut)
!                   Mar.18
              call cthick2len(trackbefmove, thick, leng, thick, jcut)
!               reset position information in MovedTrack
              call cmovestreight(leng, trackbefmove%vec%w)
              movestat = truncated
           endif
!                reset 3 momenta px, py, pz 
!                 (but assume direction is unchanged)
           call ce2p(movedtrack)
       elseif(reverse .eq. 2) then
          movedtrack%p%fm%p(4) = movedtrack%p%fm%p(4) +  energyloss
!                reset 3 momenta px, py, pz 
!                 (but assume direction is unchanged)
          call ce2p(movedtrack)
!          MoveStat = Truncated
       else
!          MoveStat = Truncated
       endif 
       return
!      **************
       entry cqelossrate(dedx, dedxf)
!      **************
!         inquire the dedt; which may be used when particle cross 
!         an obs. level and recompute the energy
       dedx = dedt
       dedxf = dedtf
       end
!      --------------------------------------------------------------
       subroutine cputdeflection
       implicit none

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

!
       real*8 dt, dispx, dispy
       logical nodef
       if(reverse .eq. 0) then
          nodef = onedim .ne. 0  .or.
     *    (mod(howgeomag, 10) .eq. 1 .and. zfirst%pos%depth .eq. 0.)
       else
          nodef = .false.
       endif
       if(.not. nodef) then
!
!         compute   magnetic deflection first. independently of scattering.
!         system is xyz;  if Efield exists, together with it
          call cmagdef
       endif
       if(.not. nodef .and. reverse .eq. 0) then
!                this is for multiple scattering
          call celecdef(dispx, dispy) ! effect alrady put in
                    ! MovedTrack. dispx, y are dummy
       endif
       if(.not. nodef) then
          call csetpos(movedtrack%pos)
          call cgetzenith(movedtrack, movedtrack%vec%coszenith)
!          reset momentum to be compatible with direction cos.
          call cresetmom(movedtrack)
!
          if(timestructure .and. reverse .eq. 0) then  ! only for mul. scat
!                compute excess path length to be added to streight path
             call cexcesslen(dispx, dispy, dt)
             if(beta .gt. 0.) then
                movedtrack%t = movedtrack%t + dt/beta
             endif
          endif
       endif
       end
!      *************************
#include "ZsubstRec.h"
       subroutine cmagdef
       use modefield
       implicit none
#include  "Zobs.h"
#include  "Zobsp.h"
#include  "Ztrack.h"
#include  "Ztrackp.h"
#include  "Ztrackv.h"
#include  "Zelemagp.h"

!

!              by  Geomag  (dr and ddirec)
       type(coord)::dispmr, dispmd
       type(magfield)::tempmag
       real*8  leng
       real*8 temp1, temp2
!cc    *   ,norm
       integer i

       integer icon
       type(coord):: middle
       logical usemiddle, high, usecmiddle
       real*8 cheight
       data cheight/30.d3/  
       real*8 newmom(3), dE, E1, E2, p1sq, p2sq, newE, norm
       save cheight

       leng = intinfarray(processno)%length
       if( howefield >= 1 ) then
          call cdefbymagande(trackbefmove,  leng,  dispmr, dispmd,
     *     newmom)
          p1sq =dot_product(trackbefmove%p%fm%p(1:3),
     *           trackbefmove%p%fm%p(1:3))
          e1 = sqrt( p1sq + trackbefmove%p%mass**2)
          p2sq =dot_product(newmom(1:3),newmom(1:3))
          e2=sqrt(p2sq +  trackbefmove%p%mass**2)
          de=e2-e1  !  may be + or -  depending on Ef and charge
          ! dE/dx loss has been put in MovedTrack. adjust it again
          
          newe = movedtrack%p%fm%p(4) + de
          if( newe <  trackbefmove%p%mass ) then
             newe =  trackbefmove%p%mass
             newmom(:) = 0
          else
             ! keeping the direction,adjust new momentum to be
             ! consistent with newE
             if( p2sq <= 0. ) then
                write(0,*) ' p2sq = ',p2sq
                write(0,*)
     *          ' TrackBefMove%p=',trackbefmove%p%fm%p(1:4)
                write(0,*) ' dE, newE=',de, newe
                write(0,*) ' leng =', leng
                write(0,*) 'dispmr=',dispmr%r(:)
                write(0,*) 'dispmd=',dispmd%r(:)
                write(0,*)
     *            'code, chg=',trackbefmove%p%code,
     *            trackbefmove%p%charge

                newe = trackbefmove%p%mass
                newmom(:) = 0
                dispmd%r(:) = (/0.,0.,1./)
             else
                norm = sqrt( (newe**2 - trackbefmove%p%mass**2)/p2sq)
                newmom(:) = newmom(:) *norm
             endif
          endif
          movedtrack%p%fm%p(4) = newe
!          note. dispmr is r(new)-r(old) vector and different
!          from other routines below. 
!           dispmd is new dir and set at 100
          movedtrack%pos%xyz%r(:) =  trackbefmove%pos%xyz%r(:) +
     *         dispmr%r(:) 
          movedtrack%p%fm%p(1:3) = newmom(:)
          goto 100
       endif
!
!            UseRungeKutta         Height       
!                0                  any         Mag and cmagneticDef

!                1                 >cheight     middle Mag and cmagneticDef
!                1                 <            Mag and cmagneticDef

!                2                 >            middle Mag and cmagneticDef 
!                                               or cbDefByRK2        
!                2                 <            Mag and cmagneticDef 

!                3                 >            middle Mag and cmagneticDef 
!                                               or cbDefByRK        
!                3                 <            Mag and cmagneticDef 

!                4                 >            use Mag at curved middle point
!                                               and estimate end point by 
!                                               cmagneticDef or cbDefByRK2
!                4                 <            Mag and cmagneticDef   

!                5                 >            use Mag at curved middle point
!                                               and estimate end point by 
!                                               cmagneticDef or cbDefByRK
!                5                 <            Mag and cmagneticDef   

!                6                 >            cbDefByRK2
!                6                 <            Mag and cmagneticDef 

!                7                 >            cbDefByRK
!                7                 <            Mag and cmagneticDef 
!
!                8      at any height           cbDefUser ; interface is
!                                               same as cbDefByRK
!
       if(userungekutta .eq. 8 ) then
          call cbdefuser(trackbefmove,  leng,  dispmr, dispmd,
     *    movedtrack%p%fm%p(1:3))
#if defined SUBSTREC             
             movedtrack%pos%xyz = dispmr ! this is not a dispalcement vector
#else 
             call csubstcoord(dispmr, movedtrack%pos%xyz)
#endif
          goto 100
       endif
       high = trackbefmove%pos%height .gt. cheight
       usemiddle = ( userungekutta .ge. 1 .and.
     *               userungekutta .le. 3 .and.
     *               high )
       usecmiddle = ( userungekutta .ge. 4 .and.
     *               userungekutta .le. 5 .and.
     *               high )


       if(usemiddle) then
          do i = 1, 3
             middle%r(i) = trackbefmove%pos%xyz%r(i)+ 
     *          leng/2 * trackbefmove%vec%w%r(i)
          enddo
          middle%sys='xyz'

          call cgeomag(yearofgeomag,  middle,
     *                tempmag, icon)
          call ctransmagto('xyz', middle,
     *        tempmag, tempmag)
       elseif(usecmiddle) then
!            get curved middle point
          call cmagneticdef(trackbefmove, mag, leng/2.0d0,
     *     dispmr, dispmd)  
          do i = 1, 3
             middle%r(i) = trackbefmove%pos%xyz%r(i) + dispmr%r(i)
          enddo
          middle%sys ='xyz'
          call cgeomag(yearofgeomag,  middle,
     *                tempmag, icon)
          call ctransmagto('xyz', middle,
     *        tempmag, tempmag)
       endif

       if(usemiddle .and. userungekutta .eq. 1) then
          mag = tempmag
       elseif( usemiddle .and.  userungekutta .eq. 2 ) then
          temp1 = tempmag%x**2 + tempmag%y**2 + tempmag%z**2
          temp2 = mag%x**2 + mag%y**2 + mag%z**2
!           if  dB/B~ dB^2/B^2/2 > 0.4 %, use RungeKutta.
          if( abs(temp1-temp2)/temp1/2.0 .gt. 4.0d-3) then
             call cbdefbyrk2(trackbefmove,  leng, dispmr, dispmd,
     *         movedtrack%p%fm%p(1:3))
#if defined SUBSTREC             
             movedtrack%pos%xyz = dispmr ! this is not a dispalcement vector
#else 
             call csubstcoord(dispmr, movedtrack%pos%xyz)
#endif
             goto 100
          else
             mag = tempmag
          endif
       elseif(usemiddle .and. userungekutta .eq. 3) then
          temp1 = tempmag%x**2 + tempmag%y**2 + tempmag%z**2
          temp2 = mag%x**2 + mag%y**2 + mag%z**2
!           if  dB/B~ dB^2/B^2/2 > 0.45 %, use RungeKutta.
!                  next 4.5d-3 is very sensitive to cpu time
!                 if it was 4.0d-3, cpu time becomes twice.
          if( abs(temp1-temp2)/temp1/2.0 .gt. 4.5d-3) then
             call cbdefbyrk(trackbefmove,  leng, dispmr, dispmd,
     *         movedtrack%p%fm%p(1:3))
#if defined SUBSTREC             
             movedtrack%pos%xyz = dispmr ! this is not a dispalcement vector
#else 
             call csubstcoord(dispmr, movedtrack%pos%xyz)
#endif
             goto 100
          else
             mag = tempmag
          endif
       elseif(usecmiddle .and. userungekutta .eq. 4) then
          temp1 = tempmag%x**2 + tempmag%y**2 + tempmag%z**2
          temp2 = mag%x**2 + mag%y**2 + mag%z**2
!           if  dB/B~ dB^2/B^2/2 > 0.1 %, use RungeKutta.
          if( abs(temp1-temp2)/temp1/2.0 .gt. 4.0d-3) then
             call cbdefbyrk2(trackbefmove,  leng, dispmr, dispmd,
     *         movedtrack%p%fm%p(1:3))
#if defined SUBSTREC             
             movedtrack%pos%xyz = dispmr ! this is not a dispalcement vector
#else 
             call csubstcoord(dispmr, movedtrack%pos%xyz)
#endif
             goto 100
          else
             mag = tempmag
          endif
       elseif(usecmiddle .and. userungekutta .eq. 5) then
          temp1 = tempmag%x**2 + tempmag%y**2 + tempmag%z**2
          temp2 = mag%x**2 + mag%y**2 + mag%z**2
!           if  dB/B~ dB^2/B^2/2 > 0.45 %, use RungeKutta.
          if( abs(temp1-temp2)/temp1/2.0 .gt. 4.5d-3) then
             call cbdefbyrk(trackbefmove,  leng, dispmr, dispmd,
     *         movedtrack%p%fm%p(1:3))
#if defined SUBSTREC             
             movedtrack%pos%xyz = dispmr ! this is not a dispalcement vector
#else 
             call csubstcoord(dispmr, movedtrack%pos%xyz)
#endif
             goto 100
          else
             mag = tempmag
          endif
       endif
!!!!!!   default comes here
       if(userungekutta .le. 5 .or. .not. high ) then
          call cmagneticdef(trackbefmove, mag, leng,
     *     dispmr, dispmd)

          do i = 1,  3
             movedtrack%pos%xyz%r(i) = trackbefmove%pos%xyz%r(i) +
     *         dispmr%r(i)
          enddo
       elseif(userungekutta .eq. 6) then
          call cbdefbyrk2(trackbefmove,  leng, dispmr, dispmd,
     *         movedtrack%p%fm%p(1:3))
#if defined SUBSTREC
          movedtrack%pos%xyz = dispmr  ! this is not a dispalcement vector
#else 
          call csubstcoord(dispmr, movedtrack%pos%xyz)
#endif
                                      !  but the vector itself.
       else
          call cbdefbyrk(trackbefmove,  leng, dispmr, dispmd,
     *         movedtrack%p%fm%p(1:3))
#if defined SUBSTREC
          movedtrack%pos%xyz = dispmr  ! this is not a dispalcement vector
                                      !  but the vector itself.
#else 
          call csubstcoord(dispmr, movedtrack%pos%xyz)
#endif
       endif
 100   continue
#if defined SUBSTREC
       movedtrack%vec%w = dispmd
#else
       call csubstcoord(dispmd, movedtrack%vec%w)
#endif
       end
!        for IBM AIX
      subroutine csubstcoord(c1, c2)
      implicit none
#include "Zcoord.h"
      type(coord)::c1, c2   
      c2 = c1
      end

!      ==============================================================
       subroutine celecdef(dx, dy)
       implicit none

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

      real*8 dx, dy             ! output. scatttering displacement
!     
      type(coord)::dircos
!              by  Multiple Scattering
      type(coord)::dsa   ! dire ccos of scattering angle
      type(coord)::w
      real*8  sint, cs, sn, tmp, avx, avy, disp, dt, dl
      real*8 r, g1, g2, gf1, gf2, beta2, tetarms
      real*8 theta
      
      call cmulscat(theta)
      if(theta .lt. 0.01d0) then
!                 cos
!     dsa%z = 1.-theta**2/2
         dsa%r(3) = 1.-theta**2/2
         sint = theta
      else
!     dsa%z = cos(theta)
         dsa%r(3) = cos(theta)
         sint = sin(theta)
      endif
!        azimuthal angle
      call kcossn(cs, sn)
!     dsa%x = sint * cs
!     dsa%y = sint * sn
!     dsa%z = cos(theta)
      dsa%r(1) = sint * cs
      dsa%r(2) = sint * sn
      dsa%r(3) = cos(theta)    
      dt = intinfarray(processno)%thickness/ x0 ! r%l
      dl = intinfarray(processno)%length !  m

!      if(Moliere .ge. 0) then
!     if( ALateCor ) then       ! >v7.645
!         v>=7.655      
      if( (moliere == 0 .and. alatecor>=1) .or. 
     *     alatecor == 2) then
!     ALateCor D=1.  if Gaussian (Moliere=0), correlation is taken
!     into account.
!     If 0, no correlation is considered.
!        2, correlation based on Gaussian formula is forced even if
!           Moliere is not 0.
         
!         sample displacement correlated to theta
!                 this is the same as P.D.B though look like
!                diff.
        tmp = dl/2.d0
        avx = tmp * dsa%r(1)
        avy = tmp * dsa%r(2)
!                dispersion
        gf1 = trackbefmove%p%fm%p(4)/movedtrack%p%mass
        gf2 = movedtrack%p%fm%p(4)/movedtrack%p%mass
        beta2 = 1.d0 - 1.d0/gf1/gf2
        if(beta2 .le. 0.) then
           disp = 0.d0
        else
           if(dt .gt. 1.d-3) then
              tetarms = es/trackbefmove%p%fm%p(4)*
     *            abs(movedtrack%p%charge)*
     *            sqrt(dt)*(1.0 + 0.038*log(dt))
           else
              tetarms = es/trackbefmove%p%fm%p(4)*
     *                  abs(movedtrack%p%charge)*
     *              sqrt(dt)
           endif
           disp=tetarms/sqrt(6.d0*beta2)*dl/2.d0
!               sample 2 independent gaussian variables
!             with mean 0 and var 1
        endif
        call kgauss2(0.d0, 1.0d0, g1, g2)
        dx = g1 * disp + avx
        dy = g2 * disp + avy
!                  displacement
        r=sqrt(dx*dx+dy*dy)     ! in m
!              direction cos of vector r in original sys.
        if(r .ne. 0.) then
!           w%x = dx/r
!           w%y = dy/r
!           w%z = 0.
           w%r(1) = dx/r
           w%r(2) = dy/r
           w%r(3) = 0.d0
           
!                 transform wx,wy,wz to original sys.
!                    TrackBefMove is better
           call ctransvectz(trackbefmove%vec%w, w, w)
!               r is already in m.
!              add scattering effect.
!              r*w is displacement by scattering
!           MovedTrack%pos%xyz%x = r*w%x + MovedTrack%pos%xyz%x
!           MovedTrack%pos%xyz%y = r*w%y + MovedTrack%pos%xyz%y
!           MovedTrack%pos%xyz%z = r*w%z + MovedTrack%pos%xyz%z
           movedtrack%pos%xyz%r(1:3) =
     *         r*w%r(1:3)+ movedtrack%pos%xyz%r(1:3)
        endif
      else
         dx = 0.
         dy = 0.
      endif
!        convert scattering angle at end of path to
!        original system . MovedTrack is better since
!        mag. def is contained there already.
      call ctransvectz(movedtrack%vec%w, dsa,
     *      movedtrack%vec%w)

      end