ciniTracking.f

Go to the documentation of this file.

      subroutine cinitracking( incident )
!          this is called after primary is fixed.
      implicit none

#include "Zmanagerp.h"
#include "Ztrack.h"
#include "Ztrackp.h"
#include "Ztrackv.h"
#include "Zobs.h"
#include "Zobsp.h"
#include "Zobsv.h"
#include "Zatmos.h"
#include "Zincidentv.h"
#include "Zearth.h"

      external cxyz2prim, cxyz2det
      integer klena, leng, icon
      real*8  zAngleSave, r1, r2, cosx, lengx
      real*8 clenbetween2h, cnewcos
      character*80 tracefile
      type(track)::incident
      type(coord)::xyz
      integer i
      data zanglesave/-1.d30/
      save zanglesave
      

      zfirst%pos%depth = 0.
      zsave = -1.d30
#if LABELING > 0
      labelcounter = 0
#endif
! 
      if(job .eq. 'flesh') then
         call csetemin(generate2, keminobs2(1), kemin2, kemincas2)
         call csetemin(generate, keminobs(1), kemin, kemincas)
      else
         call csetemin(generate, keminobs(1), kemin, kemincas)
         kemin2 = kemin
         kemincas2 = kemincas
      endif


      call cprimxyz   ! compute primary x,y,z axis in 'xyz' system
      if((trace .gt. 0 .and. trace .lt. 60) .or. 
     *    (trace .gt. 100 .and. trace .lt. 160) ) then
!            default trace output is requested. Open disk file
         write(tracefile, *) tracedir(1:klena(tracedir))//'/trace',
     *         eventno
         call kseblk(tracefile, ' ', leng)
         call copenfw(tracedev, tracefile(1:leng), icon)
         if(icon .ne. 0) then
            call cerrormsg('tracefile couldnot be opened',0)
         endif
      elseif(trace .gt. 60 .and. trace .lt. 100 .or.
     *       trace .gt. 160 .and. trace .lt. 200) then
         call cputcerenkovs     ! cerenkov output header for each event
      endif
!
      call ciniasobs
      if(abs(obsplane) .eq. horizontal  .or. 
     *        obsplane .eq. spherical) then
         call csetobsz(cxyz2det)
      elseif(abs(obsplane) .eq. perpendicular ) then
         call csetobsz(cxyz2prim)
!          inicident.where can be fixed here. not  in cmkIncidnt
         call cxyz2prim(obssites(noofsites)%pos%xyz, 
     *                  incident%pos%xyz, xyz)         
         do i = 1, noofsites
            if(xyz%r(3) .gt. obssites(i)%zpl ) then
               incident%where = i
               goto 222
            endif
         enddo
         incident%where = noofsites + 1
         incidentcopy%where = incident%where
 222     continue
      elseif(obsplane .eq.  notused ) then
!           Nothing to do. observation is at height =
!           BorderHeighL. (for neutrino)
      else
         call cerrormsg('ObsPlane value is wrong', 0)
      endif
!         check some parameters
      call  cexamparam
!         if one dimensional and large angle, use table method
!         for lenth <--> thickness conversion.  make table for that
      usetbl = onedim .eq. 3 .or.
     *  ( onedim .eq. 2  .and. 
     *    abs( angleatobscopy%r(3) ) .lt. 0.5 )
      usetbl = usetbl .and.  .not. upgoing
      if(usetbl .and. zanglesave .ne. angleatobscopy%r(3) ) then
         hbase = heightlist(noofsites) - 0.3d0  ! make little bit smaller
         htop =100.d3   ! old value was 30d3
         r1 = htop + eradius
         r2 = hbase + eradius
         lengx = clenbetween2h(r2, r1, -angleatobscopy%r(3))
         cosx = cnewcos(r2, -angleatobscopy%r(3), lengx)
         call cl2ttbl(htop, hbase, cosx, -angleatobscopy%r(3),
     *     lenstep, 
     *     lentbl, heighttbl, costbl,  thicktbl, maxl2t, numstep)
         zanglesave = angleatobscopy%r(3)
      endif
!          fix B at the starting point
      if( howgeomag .le. 2 .or. howgeomag .eq. 31 ) then
         call cgeomag(yearofgeomag, incident%pos%xyz, mag, icon)
         call ctransmagto('xyz', incident%pos%xyz,  mag, mag)
      else
         mag = magfieldxyz
      endif

      end
!     --------- compute primary system, x,y,z.  
!               The deepest detector is the reference.
      subroutine cprimxyz
      implicit none
#include "Ztrack.h"
#include "Zincidentv.h"
#include "Zobs.h"
#include "Zobsv.h"
!           all are unit vector in 'xyz' system.
!
!        X = Z x V  ( Z  is primary direction; V is vertical axis )
!          =  V x(-Z)
!          =  V x DcAtObsXyz = DetZaxis.r(1) DcAtObsXyz
!        Y = Z x X  = X x (-Z)
!          = X x DcAtObsXyz
!              
!      compute Xprimary, Yprimary, Zprimary in 'xyz' system.

      real*8  temp
!      
      zprimary%r(:) = - dcatobsxyz%r(:)
      call cvecprod(zprimary, detzaxis, xprimary)
!         see if Zprimary // DetZaxis; if so reset Xprimary
      temp= sum( xprimary%r(:)**2 )
      if(temp .lt. 1.e-12) then
         xprimary = detxaxis
      else
         temp =sqrt(temp)
!                ifort v14 dies here if (:) is used  
         xprimary%r(1) = xprimary%r(1)/temp         
         xprimary%r(2) = xprimary%r(2)/temp         
         xprimary%r(3) = xprimary%r(3)/temp         
      endif
      call cvecprod(zprimary, xprimary, yprimary)
      txyz2prim(1,:) = xprimary%r(:)
      txyz2prim(2,:) = yprimary%r(:)
      txyz2prim(3,:) = zprimary%r(:)
      tprim2xyz(:,:) = transpose(txyz2prim(:,:))
       
      end
!     ******************************
      subroutine csetobsz(converter)
!     *****************************
!      compute Z of observation plane 
!      
      implicit none
#include "Zcoord.h"
#include "Zpos.h"
#include "Zmagfield.h"
#include "Zobs.h"
#include "Zobsp.h"
#include "Zobsv.h"

      external converter  ! cxyz2prim or cxyz2det
!                           depending on |ObsPlane| =2, 1
!                           or  cxyz2det for ObsPlane=3
      integer i
      type(coord)::temp
!        compute z from the base detector.
      do i = 1, noofsites
         call converter(obssites(noofsites)%pos%xyz, 
     *     obssites(i)%pos%xyz,
     *     temp)
         obssites(i)%zpl = temp%r(3)
      enddo
!            this is for A.S
      do i = 1, noofassites
         call converter(asobssites(noofassites)%pos%xyz, 
     *     asobssites(i)%pos%xyz,
     *     temp)
         asobssites(i)%zpl = temp%r(3)
      enddo
      end
!      **************************
      subroutine cqfirstid(depth)
!         inquire the first interaction depth
      implicit none

#include "Ztrack.h"
#include "Ztrackv.h"
      real*8 depth   ! output. to get the first I%D
      
      depth = zfirst%pos%depth   ! vertical depth
      end
!      **************************
      subroutine cqfirstipi(ptrack)
!         inquire the complete first interaction point info. 
!         as a track of the primary.
      implicit none

#include "Ztrack.h"
#include "Ztrackv.h"
      type(track)::ptrack   ! output.    if ptrack%pos%depth=0.
                              !  no interaction has been occurred
      
      ptrack = zfirst
      end
      subroutine cqfirstcolmedia(A, Z, xs)
!        retrns first col. element and Xsecion on it
      implicit none
#include "Ztrack.h"
#include "Ztrackv.h"
      integer,intent(out):: A
      integer,intent(out):: Z
      real(8),intent(out):: xs
      a = firstcola
      z = firstcolz
      xs = firstcolxs
      end
! --------------------------------
      subroutine ciniasobs
      implicit none

#include "Zcoord.h"
#include "Zpos.h"
#include "Zmagfield.h"
#include "Zobs.h"
#include "Zobsp.h"
#include "Zobsv.h"

      integer i

      do i = 1, noofassites
         asobssites(i)%esize = 0.
         asobssites(i)%age = 0.
      enddo
      end
!     ***********************
      subroutine  csetemin(gen, eminob,  emin, emCas)
!     ***********************
      implicit none
#include  "Zmanagerp.h"
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackp.h"
#include  "Ztrackv.h"
#include  "Zincidentv.h"
#include  "Zheavyc.h"
#include  "Zheavyp.h"
!
      character*(*)  gen ! input. a copy of Generate. 
      real*8  eminob   !   input. minim observational energy (kientic)
      real*8  emin     !  output. Minimum energy to be followed for non e-g
                       !  ptcls.
      real*8  emCas    !  output. Minimum energy to be followed for e-g
!
!
      real*8  ergpn    ! energy / nucleon
      logical cas      ! 
      logical obas     ! 

!
      cas = index(gen, 'em') .gt. 0 
      obas = index(gen, 'as') .gt. 0 .or.
     *   index(gen, 'lat') .gt. 0


!          energy / nucleon for heavy
      if(incidentcopy%p%code .eq. kgnuc) then
         ergpn =
     *        incidentcopy%p%fm%p(4)/incidentcopy%p%subcode
      elseif(incidentcopy%p%code .ge. kalfa .and.
     *     incidentcopy%p%code .le.  khvymax) then
!           will  not come here
         ergpn =
     *        incidentcopy%p%fm%p(4)/code2massn(incidentcopy%p%code)
      else
         ergpn = incidentcopy%p%fm%p(4)
      endif
      if(obas ) then
               ! wait until electon  energy becomes < EasWait for AS
               ! generation. WaitRatio < 1 should be used when 1ry is 
               ! g or e.   (say, 0.01)
         easwait =ergpn * waitratio
      endif


      emcas = 1.d30

      if(obas) then
         emcas = easwait
      endif
      if(cas) then
         emcas =min(eminob, emcas)
      endif
      emin = eminob
      if(obas) then
         emin = min(emin, max(ratiotoe0* ergpn, 1.d0) )
      endif 

      if(thinsampling) then
!           thin sampling
         if(ethinratio(1) .lt. 0.) then
            ethin(1) = abs(ethinratio(1))
         else
            ethin(1) = ethinratio(1)*ergpn
         endif
         if(ethinratio(3) .lt. 0.) then
            ethin(3) = abs(ethinratio(3))
         elseif(ethinratio(3) .eq. 0.) then
            ethin(3) = ethin(1)/10.
         else
            ethin(3) = ethinratio(3)*ergpn
         endif


!         if(EthinRatio(2) .lt. 0.) then
!            Ethin(2) = abs(EthinRatio(2))
!         else
!            Ethin(2) = EthinRatio(2)*ergpn
!         endif
!             maximum weight above which no thinning
          ethin(2) =  ethinratio(2)
          if( ethinratio(4) .eq. 0.) then
             ethin(4) =  ethin(2)/10.
          else
             ethin(4) = ethinratio(4)
          endif
!///////////////////
!          write(0, *) ' ****** Ethin1,2=', Ethin
!          write(0, *) ' ergpn=', ergpn, ' EthinRatio=',EthinRatio
!////////////////

      else
           ethin(1) = 0.
           ethin(2) = 0.
           ethin(3) = 0.
           ethin(4) = 0.
      endif
      end
      subroutine cexamparam
      implicit none
#include "Ztrackp.h"
#include "Zelemagp.h"

      end