chookFlesh.f

Go to the documentation of this file.

#include "cmain.f"
#include "chookHybAS.f"
#include "ctemplCeren.f"
!
!  *************************************** hook for Beginning of a Run
!  * At this moment, all (system-level) initialization for this run
!  * has been ended.  After this routine is executed, the system goes into the
!  * event creation loop.
!  *
      subroutine chookbgrun
      implicit none
#include "Zmanagerp.h"
#include "Zprivate.h"

      real*8  temp
      character*100 msg
      integer klena
      integer i
!     ==================================================
      
      integer seed(2)
!     ==================================================

      eventno = 0
      realbegin = .true.
      topofnode = .true.

!            namelist output
      call cwriteparam(errorout, 0)
!            primary information
      call cprintprim(errorout)
!            observation level information
      call cprintobs(errorout)

      call cquhooki(1, mdev)      ! get skeleton memo dev #
      call cquhooki(6, howflesh)  ! how to flesh. non zero-->orthodox
      call cquhookc(1, msg)       ! get file name for sekelton memo
      call cgetfname(msg, mskel)  ! add host name etc if needed

!      do i = 1, NoOfSites
!         write(*, *) ' depth ',
!     *    sngl(ObsSites(i).pos.depth/10.)
!      enddo

      open(mdev, file=mskel(1:klena(mskel)), form='unformatted',
     *  status='old' )

      end

!     *********************************** hook for Beginning of  1 event
!     *  All system-level initialization for 1 event generation has been
!     *  eneded at this moment.
!     *  After this is executed, event generation starts.
!     *
      subroutine chookbgevent
      implicit none
#include "Zprivate.h"


      integer nomore
      if( realbegin ) then
         call cbegin1ev( nomore )
         if( nomore .eq. 1) then
            call cerrormsg('all events are fleshed', 1)
            stop  !!!!!!!!!!!!  
         endif
         topofnode = .true.
      endif
      call c1by1


      end
      subroutine cbegin1ev(nomore)
      implicit none
#include "Zprivate.h"
#include "Ztrack.h"
#include "Ztrackv.h"
#include "Ztrackp.h"
#include "Zobs.h"
#include "Zobsp.h"
#include "Zobsv.h"
#include "Zcode.h"
#include "Zmanager.h"
#include "Zmanagerp.h"
      
      integer nomore       !  output. 0 still there  are showers
                           !          1 no more skeleton showers to be fleshed
!          event number, primary      

      type(track)::incident
      type(track)::zsave
      type(coord)::angle
      
      integer i
      integer seed(2)
      integer cumnum, num, jeof, fin
      read( mdev, end=1000 ) cumnum, num, seedsave, zfirst

      eventsintherun = eventsintherun + 1
      eventno = eventno + 1
!          get random seed at skelton making; this can work
!          if seed file is supplied
!      call creadSeed(SeedSave, EventNo, jeof)
!      if( jeof .ne. 0 ) goto 1000

!                 reset the seed.
      call rnd1r(seedsave)

!         next incident; confirmed to be the same one as preserved one
      call cmkincident(incident, fin)

      if(fin .ne. 0 ) goto 1000
      zsave = zfirst    ! save;  this is reset in next 
      call cinitracking( incident )   
!          set first interaction pos
      zfirst = zsave
      call cresettimer(zfirst)

      realbegin = .false.

!          do your own init for a one event here
!      ==========================================================
         call cqincident( incident, angle)
      write(*,999) 
     *    eventno, 
     *    incident.p.code,
     *    incident.p.subcode,
     *    incident.p.charge, 
     *    incident.p.fm.e,
     *    -angle.r(1),
     *    -angle.r(2),
     *    -angle.r(3)
 999      format(i9,3i4,e15.5,3(1x,f12.8))


!      ==========================================================
!

      call cgethes(mdev)  ! get high energy ptlcs
      call cobshes        ! imitate their observation
      nomore = 0
      return

 1000 continue
      nomore = 1
      end

!     ************************************ hook for observation
!     *  One particel information is brought here by the system.
!     *  All information of the particle is in aTrack
!     *
      subroutine chookobs(aTrack, id)
!
!     Note that every real variable is in double  precision so
!     that you may output it in sigle precision to save the memory.
!     In some cases it is essential to put it in sigle (say,
!     for gnuplot).
! 
      implicit none
#include "Zcode.h"
#include "Ztrack.h"
#include "Zprivate.h"
      integer id  ! input.  2 ==> reached at an observation level
!                           1 ==> aTrack is going out from
!                                 outer boundery.
!                           2 ==> reached at an observation level
!                           3 ==> reached at inner boundery.
      type(track)::aTrack
!
!     For id =2, you need not output the z value, because it is always
!     0 (within the computational accuracy).
!
      if(id .eq. 2 .and. atrack.p.code .ne. kneumu .and.
     *   atrack.p.code .ne. kneue) then

!     ===================================================

       if( atrack.p.code .le. 6 .and. atrack.p.code .ne. 3 ) then
!         write(*, 959) 
!     *  aTrack.where,   
!     *  aTrack.p.code,   
!     *  aTrack.p.charge,  
!     *  sngl( aTrack.p.fm.p(4) - aTrack.p.mass ),  
!     *  sngl( aTrack.pos.xyz.r(1) ), 
!     *  sngl( aTrack.pos.xyz.r(2) ) ,  
!     *  sngl( aTrack.vec.w.r(1) ),   
!     *  sngl( aTrack.vec.w.r(2) ),  
!     *  sngl( aTrack.vec.w.r(3) ),  
!     *  sngl( aTrack.vec.coszenith )  
! 959      format(3i3,f12.3,2f16.6,4(1x,f12.8))
      endif

!     ===================================================

!         write(*,'(4i5, g15.4,g15.3)') 
!     *   aTrack.where, aTrack.p.code, aTrack.p.subcode,
!     *   aTrack.p.charge, sngl( aTrack.t ),
!     *   sngl( aTrack.p.fm.p(4) - aTrack.p.mass)
!     *   sngl( aTrack.pos.xyz.r(1) ), sngl( aTrack.pos.xyz.r(2) ),
!     *   sngl( aTrack.vec.w.r(1) ), sngl(aTrack.vec.w.r(2) ),
!     *   sngl(aTrack.vec.w.r(3) ), 
!     *   sngl(aTrack.vec.coszenith)

      endif
      end

!    *********************************** hook for end of 1 event
!    * At this moment, 1 event generation has been ended.
!    *
      subroutine chookenevent

      implicit none
#include "Zprivate.h"
#include "Ztrack.h"
#include "Ztrackv.h"
#include "Zobs.h"
#include "Zobsp.h"
#include "Zobsv.h"

      integer i
      
      if(realend) then
         call cfintracking
!           real end of 1 event; if you need to do some here is
!           the place
!         ========================================================

         if(observeas) then
!                   electron size in B approx.
            do i = 1, noofassites
               write(*, *) asobssites(i).age, asobssites(i).esize
            enddo
         endif



!        ========================================================
!  
      else
!           there is still low energy skeleton ptcls
!           nothing to do here
      endif

      end


!     ********************************* hook for end of a run
!     *  all events have been created or time lacks
!     *
      subroutine chookenrun
      implicit none
#include "Zprivate.h"
!     =========================================================

!     =========================================================
      call  cprintstatus   ! if don't like,  comment out
      end
!     ********************************* hook for trace
!     *  This is called only when trace > 100
!     *  User should manage the trace information here.
!     *  If you use this, you may need some output for trace
!     *  at the beginning of 1 event generatio and at the end of  1 event
!     *  generation so that you can identfy each event.
!     *
!     *
      subroutine chooktrace
            implicit none

#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include  "Zobs.h"
#include  "Zobsv.h"

       real*4 h1,  h2
!
!    Every time a particle is moved in the atmosphere, this routine is called,
!    if trace > 100
!         For a one track segment,
!     TrackBefMove  has  track information at the beginning of the segment.
!     MoveTrack    has   track information at the end of the segment.
!   
!     You can know the  information a track contains in the 
!     chookObs routine. (Note however, no conversion of coordinate
!     has been done.  The values are in the Earth xyz system.)
!     Besides quantities explained there, you can use, for a  given 'track'
!
!     atrack.pos.xyz.x, atrack.pos.xyz.y, atrack.pos.xyz.z    (x,y.z)
!     atrack.pos.radiallen   (distance from the center of the earth)
!     atrack.pos.depth       (vertical depth)
!     atrack.pos.height      (vertical heigth from sea level)  
!

      h1 = trackbefmove.pos.height- obssites(noofsites).pos.height
      h2 = movedtrack.pos.height - obssites(noofsites).pos.height

      end
!     ********************* this is the hook called when
!       an electron made an interaction.
!
      subroutine chookeint(never)
      implicit none
      integer never             ! input & output
      never = 1
      end

!     ********************* this is the hook called when
!       a gamma ray made an interaction.
!
      subroutine chookgint(never)
      implicit none
      integer never   ! input & output
      never = 1
      end

!     ********************* this is the hook called when
!       non e-g particle made an interaction.
!
      subroutine chooknepint(never)
      implicit none
      integer never   ! input & output
      never = 1
      end


      subroutine cgethes(from)
      implicit none
#include "Zprivate.h"
      integer from

      integer i

      read(from)  np
      do i = 1, np
         read(from) o(i)
      enddo
      end

      subroutine cobshes
      implicit none
#include "Zprivate.h"
#include "Ztrack.h"
!
!           memorized high energy showers at the skeleton making
!     time is put into the chookObs as if they are really observed
      type(track)::aTrack

      integer i

      do i = 1, np
         atrack.where =  o(i).where 
         atrack.p.code =  o(i).code 
         atrack.p.subcode = o(i).subcode 
         atrack.p.charge = o(i).charge 
         atrack.t = o(i).atime 
         atrack.p.fm.p(4) = o(i).erg
         atrack.p.mass = o(i).mass 
         atrack.pos.xyz.r(1) = o(i).x 
         atrack.pos.xyz.r(2) = o(i).y 
         atrack.vec.w.r(1) = o(i).wx
         atrack.vec.w.r(2) = o(i).wy
         atrack.vec.w.r(3) = o(i).wz
         atrack.vec.coszenith = o(i).zenith 
         call chookobs(atrack, 2)
      enddo
      end


!        process low energy partilces in the skeleton 1 by 1

      subroutine c1by1
      implicit none
#include "Zprivate.h"
#include "Ztrack.h"
#include "Ztrackv.h"

      character*100 msg

      call cinitstack  ! empty the stack

      if( topofnode ) then
         read(mdev)  nooflowe,  p
         if( p.asflag .eq. -1 .and. observeas 
     *      .and. howflesh .ne. 0 )  then
            call embedas
         endif
         nlowcounter = 0
         if( nooflowe .eq. -1 ) then
            realend = .true.
            realbegin = .true.
            return    ! ************
         endif
      endif

      realbegin = .false.
      realend = .false.


      if( nlowcounter .eq. nooflowe ) then
         topofnode =.true.
         return   ! **************
      endif

      topofnode = .false.
!         still not the  end of 1 event

      read(mdev)  c

      nlowcounter = nlowcounter + 1
      call cmove_c_stack   ! move c into stack

      end
!
      subroutine embedas
      implicit none
#include "Zprivate.h"
#include "Ztrack.h"
#include "Zearth.h"

      type(track)::el
      
      el.pos.depth = p.depth
      el.vec.coszenith = p.coszenith
      el.pos.radiallen = p.height + eradius
      el.pos.height = p.height
      el.p.fm.p(4) = p.erg
      el.wgt = 1.0
      call cobas(el)
      end

      
      

      subroutine cmove_c_stack
      implicit none

#include "Zprivate.h"
#include "Ztrack.h"
#include "Zearth.h"

      type(track)::aTrack
!
!          a child of the current parent is moved to stack
!      as a track info.
!
      atrack.pos.xyz.r(1) = p.posx
      atrack.pos.xyz.r(2) = p.posy 
      atrack.pos.xyz.r(3) = p.posz
      atrack.pos.depth = p.depth 
      atrack.pos.height = p.height
      atrack.pos.colheight = p.colheight 
      atrack.t = p.atime 

      atrack.where = p.where 

      atrack.p.code = c.code 
      atrack.p.subcode = c.subcode
      atrack.p.charge  = c.charge
      atrack.p.fm.p(1) = c.fm(1) 
      atrack.p.fm.p(2) = c.fm(2) 
      atrack.p.fm.p(3) = c.fm(3) 
      atrack.p.fm.p(4) = c.fm(4) 
      atrack.p.mass = c.mass

!        --------------- next must be compute here

      atrack.pos.radiallen = eradius +atrack.pos.height
      atrack.pos.xyz.sys = 'xyz'
      atrack.vec.w.sys = 'xyz'
      atrack.wgt = 1.0
      if(p.asflag .ne. 0 .and. howflesh .ne. 0 ) then 
!           asflag=   -1 or -2 may happen
         atrack.asflag = 1
      else
!           parent is made to be unconceivable, so
!           don't put flag for AS generation
         atrack.asflag =  0
      endif


         

      call cresetdirec( atrack )    ! set vec.w and coszenith 

      call cpush(atrack)
      end