chookASbyH.f

Go to the documentation of this file.

#include "../chookHybAS.f"
#include "../ctemplCeren.f"
#include "ZcosmosBD.h"
!
!     chookASbyH:  This is a special user hook to compute
!                  AS generated by heavy primaries with
!                  generate='qas'.
!  **Important**  
!      You almost need not modify this file, although
!      you have to read this.  Those parts that you 
!      have to modify are gathered in chookASbyH2.f
!      
!
!       If the primary is a proton, we need about only 10 min even at
!       10^21 eV to produce AS by the hybrid method.  The time needed
!       at lower energies is shorter (but not propotinal to the energy).
!         However, for heavy primaries such as Fe, 
!       we need few to several hours to simulate 1 A.S even at 10^18 eV,
!       so we must coin out some new method.
!       The idea here is to use chookNEPInt where we get control 
!       when a hadron makes an interaction.  We may see interacting
!       nucleons there, and replace each of their further development by 
!       a precomputed proton shower (let's call it componet shower)
!       and discard further development of the shower by that nucleon.
! 
!       To use this routine for a heavy primary of mass A with total energy E,
!       you need to create AS by protons with energy E/A beforehand.
!       The first iteraction point of the protons should be fixed somewhere,
!       or should be known.
!       (To fix the first collision point when making c.s,
!        put Freec = f, and  say, HeightOfInj =20000). 
!       (As far as the one  dimensional  dvelopment of electrons is 
!       concerned, the first interaction point  can be almost arbitrary; 
!       The LPM effect will not be the problem).
!
!         The number of component showers should be 100~300.  For each of
!       component showers,  you must have the first interaction point, 
!       air shower sizes at different depths  (better at equisteped depths)
!       In some case you need  ages, muon size, etc, depending on your demands.
!       An arbitrary shower among them should be randomly selected as 
!       a component shower.  For this end, I recommend you to use
!       direct access file so that  you may read an arbitrary c.s from 
!       a disk file  when a c.s is requested.  Alternatively, you may
!       define a large array to store c.s data so that quick random
!       selection is possible.
!
!    Thus,  (in chookASbyH2.f)
!     a)  In chookBgRun, 
!           you may open the component shower file.
!           get some overall information
!     b)  In chookBgEvent,
!          you may clear the AS size counters.  This is auomatically
!         performed if you use standard counters.
!     c)  In chookNEPInt, 
!         c-1) you should inquire the interacting nucleons by
!              calling  cqIntePtcl,
!         c-2) get a component shower randomly
!         c-3) increment the size counters
!              You need some interpolation for this end.
!         c-4) Finally you must give "never = 2" to discard further
!              development of showers  by the interacting nucleons.
!     d)  In chookEnEvent,
!          you should output AS size.
!
!     In this template, we assume that log10( electron size) is stored in
!     a direct access file.  If size is 0, log10(size) may be -1.
!
!    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!      All actual interface routines are in chookASbyH2.f
!      In this file, there are lines to call routines in chookASbyH2.f
!      see %%%% parts below.
!
!  *************************************** 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"

!
!         If you feel writing the parameters on stderr is
!         a bother, comment out the next or
!         use other device than ErrorOut.
!         Also you may comment out all output routines below.
#ifdef sun4
      external csighandler
      integer  ieeer, ieee_handler

      ieeer = ieee_handler('set', 'invalid', csighandler)
#endif
!
!            namelist output
      call cwriteparam(errorout, 0)
!            primary information
      call cprintprim(errorout)
!            observation level information
      call cprintobs(errorout)
!     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!      open component disk file, and get some over all
!     informaiton such as number of component showers
     
      call chookbgrunas 
      
      end

#ifdef sun4
      integer function csighandler(sig, code, context)
      implicit none
#include "Zmanagerp.h"
      integer sig, code, context(5)
      write(errorout, *)  ' f.p exception content=' , context(4)
!      call abort()
      end
#endif

!     *********************************** 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
!
      integer seed(2)
!      write(*, *) ' bigin event generation'
      call cqinirn(seed)
!      write(*,*) ' seed=', seed


      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"
      integer id  ! input.  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) then
!            output typical quantities.
!        write(*, *) 
!     *  aTrack.where,   !  observation level. integer*2.  1 is highest.
!     *  aTrack.p.code,    !  ptcl code.  integer*2.
!     *  aTrack.p.charge,  !  charge,  integer*2 
!     *  sngl(aTrack.t), !  relateive arrival time in nsec (NOT sec).
!                        !  if TimeStructure is F, nonsense.
!     *  sngl(aTrack.p.fm.e)  ! total energy in GeV.
!     *  sngl(aTrack.pos.xyz.x), sngl(aTrack.pos.xyz.y),  !  x, y in m
!     *  sngl(aTrack.vec.w.x),  ! direc. cos.x in the current detector system.
!     *  sngl(aTrack.vec.w.y),  ! direc. cos.y
!     *  sngl(aTrack.vec.w.z),  ! direc. cos.z
!     *  sngl(aTrack.vec.coszenith) ! cos of zenith angle
!         if(aTrack.p.code .eq. kelec) then
!            write(*, *) aTrack.where
!         endif
!      endif
!         you may need in some case other information such as
!       aTrack.p.subcode   ! sub code of the particle integer*2
!       aTrack.p.mass      ! mass 
!       aTrack.wgt         ! weight of the particle (may not be 1. if
!                           ! ThinSampling =T)
!       aTrack.p.fm.x      ! momentum x component.  Note. Momentum is
!                            given in the  Earth xyz system.

!       aTrack.p.fm.y      !          y
!       aTrack.p.fm.z      !          z

      endif
      end

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

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



!      type(track):: inci
!      type(coord):: angle
      integer i
      real*8 fdepth

!

      if(observeas) then
         call cqfirstid(fdepth)
!         call cqIncident(inci, angle)
!         bsin = cgetBsin(inci.p, Mag)*1.e4
!                   electron size in B approx.
!         write(*, *) (ASObsSites(i).esize, i=1, NoOfASSites)
!                   size weighted age
!         write(*, *) (ASObsSites(i).age,   i=1, NoOfASSites)
         do i = 1, noofassites
             write(*, *)  sngl(asobssites(i).pos.depth),
     *           sngl(asobssites(i).esize),
     *           sngl(asobssites(i).age)
!     *          ,sngl(fdepth)
!     *          ,sngl(bsin)
         enddo
         write(*, *)
      endif

 
      end


!     ********************************* hook for end of a run
!     *  all events have been created or time lacks
!     *
      subroutine chookenrun

      implicit none
      end
!     ********************************* hook for trace
!     *  This is called only when trace > 60
!     *  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 > 60. 
!         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

#include  "Ztrack.h"
#include  "Ztrackv.h"
!  #include  "Ztrackp.h"
      
      integer never   ! input & output
      
!         don't make never = 1, if you want to get
!         information after an electron made interaction
!         if this is made non zero, this routine will never be called.
!
!   MovedTrack is the electron that made interaction
!   Pwork contains produced particles.
!   Nproduced has the number of particles in Pwork
!   IntInfArray(ProcessNo) contains the type of interaction
!
!        default setting
      never = 1
!
!        IntInfArray(ProcessNo).process will have one of
!       'brems', 'mscat', 'bscat', 'anihi' or 'mbrem'
!
      end

!     ********************* this is the hook called when
!       a gamma ray made an interaction.
!
      subroutine chookgint(never)
            implicit none

#include  "Ztrack.h"
#include  "Ztrackv.h"
!  #include  "Ztrackp.h"
      
      integer never   ! input & output
      
!         don't make never = 1, if you want to get
!         information after a gamma ray made interaction
!         if this is made non zero, this routine will never be called.
!
!   MovedTrack is the gamma that made interaction
!   Pwork contains produced particles.
!   Nproduced has the number of particles in Pwork
!   IntInfArray(ProcessNo) contains the type of interaction
!
!        default setting
      never = 1
!         IntInfArray(ProcessNo).process will have one of
!        'pair', 'comp', 'photoe' 'photop' 'mpair'
!       
      end

!     ********************* this is the hook called when
!       non e-g particle made an interaction.
!
      subroutine chooknepint(never)
            implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
!  #include  "Ztrackp.h"
      
      integer never   ! input & output
      
!         don't make never = 1, if you want to get
!         information after a non-e-g particle  made interaction.
!         if this is made non zero, this routine will never be called.
!
!   MovedTrack is the particle that made interaction
!   Pwork contains produced particles.
!   Nproduced has the number of particles in Pwork
!   IntInfArray(ProcessNo) contains the type of interaction
!
!        IntInfArray(ProcessNo).process  will have
!             'col' or 'decay'
!   %%%%%%%%%%%%%%%%%%%%%%%%%%
!         replace interacting nucleons by c.s. give never = 2
!     so that no further development by the interacting nucleons.
!
      call chooknepinta
      never = 2
      end