cputCerenkov.f

Go to the documentation of this file.

!         ready made output for cerenkov 
!         threshold dependence on the air density  is
!         taken into account.
!
       subroutine cputcerenkov
       implicit none
#include "Zmanagerp.h"
#include "Ztrack.h"
#include "Ztrackp.h"
#include "Ztrackv.h"
#include "Zcode.h"
#include "Zmass.h"
#include "Zprimary.h"
#include "Zprimaryv.h"
#include "Zheavyp.h"
#include "Zincidentv.h"
!
      integer,parameter:: maxcerenka=6
      real*8 crnth(maxcerenka)
!
!  
       real*8  thg,  thpi, thk, thmu, thnuc
       parameter(thg=40.-1., thpi=thg*maspic, thk=thg*maskc,
     1  thmu=thg*masmu, thnuc=thg*masp)
!                threshold energy in GeV at sea level
!
        type(coord)::f
        type(coord)::t
        integer chrg, code
!
        real*8 h1, den1, e1, the
!        real*8 h2, den2, e2
        real*8 cvh2den
        logical first/.true./
        real*8 denAtSea
        save first, denatsea
        integer ka, itb, it,  utrace
        data crnth(kelec)/20.e-6/, crnth(kpion)/thpi/,
     1  crnth(kkaon)/thk/, crnth(kmuon)/thmu/,
     2  crnth(knuc)/thnuc/
        
        if(first) then
           denatsea = cvh2den(0.d0)
           first = .false.
        endif

        utrace = trace
        if(utrace .gt. 160) utrace = utrace - 100
!          get air density.
        h1 = trackbefmove%pos%height
        den1 = cvh2den(h1)
!        h2 = MovedTrack.pos.height
!        den2 = cvh2den(h2)
        ka = trackbefmove%p%code
        e1 = trackbefmove%p%fm%p(4)
!        e2 = MovedTrack.p.fm.p(4)
!
        if( ka .gt. knuc) then
           the = thg * trackbefmove%p%mass
        else
           the = crnth(ka)
        endif
!
        if( e1 .gt. the*denatsea/den1 ) then
!               get transformed coord.
!  utrace:      61-70:  --> 1ry
!               71-80:  --> 1ry but z is depth above 
!               81-90:  --> det 
!               91-100 -->  det but z is depth above
!
           call ccoordfortr(utrace-60, f, t)
           itb=trackbefmove%t*100  ! centimeter/beta
           it=(movedtrack%t - trackbefmove%t)*100  !  length/beta (cm)
           chrg = trackbefmove%p%charge

           if(trace .gt. 160) then
              call chookceren(ka, chrg, e1, itb, it, f, t)
           else
              if(mod(utrace, 2) .eq. 0 ) then
                 write(tracedev) ka, chrg, sngl(e1), itb, it, 
     *           sngl(f%r(1)), sngl(f%r(2)), sngl(f%r(3)),
     *           sngl(t%r(1)), sngl(t%r(2)), sngl(t%r(3))
              else
                 write(tracedev, *) ka, chrg, sngl(e1), itb, it, 
     *           sngl(f%r(1)), sngl(f%r(2)), sngl(f%r(3)),
     *           sngl(t%r(1)), sngl(t%r(2)), sngl(t%r(3))
              endif
           endif 
       endif
       return
!        *************  start of 1 shower for cerenkov trace(called from
!                      ciniTracking
!
       entry cputcerenkovs
       utrace = trace
       if(utrace .gt. 160) utrace = utrace -100

       if(trace .gt. 160) then
          call chookcerens(eventno,  prim, angleatobscopy)
       else
          code = prim%particle%code
          if(mod(utrace, 2) .eq. 0) then
             write(tracedev) eventno,
     *       code, 
     *       prim%particle%fm%p(4), 
     *       angleatobscopy%r(1), angleatobscopy%r(2),
     *       angleatobscopy%r(3)
          else
             write(tracedev, *) eventno,
     *      code, 
     *      prim%particle%fm%p(4), 
     *      angleatobscopy%r(1), angleatobscopy%r(2),
     *      angleatobscopy%r(3)
          endif
       endif
       return
!        ************** end of  1 shower
         entry  cputcerenkove
!        ***************
         utrace = trace
         if(utrace .gt. 160) utrace = utrace - 100
         ka = 0
         chrg = 0
         e1 = 0.
         itb =0
         it =0
         f%r(1) = 0.
         f%r(2) = 0.
         f%r(3) =0
         t%r(1) = 0.
         t%r(2) = 0.
         t%r(3) = 0.
         if(trace .gt. 160) then
            call chookcerene(ka, chrg, e1, itb, it, f, t)
         else
!
            if(mod(utrace, 2)  .eq. 0) then
               write(tracedev) ka, chrg, sngl(e1), itb, it,
     *         sngl(f%r(1)), sngl(f%r(2)), sngl(f%r(3)),
     *         sngl(t%r(1)), sngl(t%r(2)), sngl(t%r(3))
            else 
               write(tracedev, *) ka, chrg, sngl(e1), itb, it,
     *         sngl(f%r(1)), sngl(f%r(2)), sngl(f%r(3)),
     *         sngl(t%r(1)), sngl(t%r(2)), sngl(t%r(3))
            endif
         endif
         end