chookSkel.f

Go to the documentation of this file.

#include "cmain.f"
#include "chookHybAS.f"
#include "ctemplCeren.f"
#if defined (KEKB) || defined (KEKA)
#define DOMPI
#endif
!
!  *************************************** 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 "Zmaxdef.h"
#include "Zmanagerp.h"
#include "Zprivate.h"
#include "Ztrack.h"
! #include "mpif.h"
! #include "Zmpi.h"
#if defined (DOMPI)
#include "mpif.h"
#include "Zmpi.h"
      real*8 etime1
      etime1 = mpi_wtime()
      write(0,*) 'rank=',mpirank, ' time start=',etime1
#else
      integer::mpirank=-1
#endif
      character*10  append
      character*100 msg
      logical  ex, opn
      integer klena
      if(mpirank .le. 0) then
!            namelist output
         call cwriteparam(errorout, 0)
!            primary information
         call cprintprim(errorout)
!            observation level information
         call cprintobs(errorout)
      endif



      call cquhooki(1, mdev)      ! get skeleton memo dev #
      call cquhooki(2, wdev)      ! get working disk dev #


      call cquhooki(3, ngmin)     ! get Nh min
      call cquhooki(4, nhmin)     ! get Ng min
      call cquhooki(5, where)     ! where to check 

      call cquhookr(1, sumegmin)      ! sum E min
      call cquhookr(2, sumehmin)

      

      call cquhookc(1, msg)       ! get file name for sekelton memo
      call cgetfname(msg, mskel)  ! add host name etc if needed
      call cquhookc(2, msg)       ! get file name for working
      call cgetfname(msg, wskel)  ! add host name etc if neededn
      call cquhookc(3, append)    ! append data, if Mskel already exists

      write(msg, *) 'Skeleton is  judged at obs.pos=', Where
      call cerrormsg(msg, 1)
      write(msg, *) ' Ngmin=',ngmin,
     *     ' SumEGmin=',sumegmin/1000.,'TeV'
      call cerrormsg(msg, 1)
      write(msg, *) ' Nhmin=',nhmin,
     *    ' SumEHmin=',sumehmin/1000.,'TeV'
      call cerrormsg(msg, 1)

!     
      inquire(file=mskel(1:klena(mskel)), opened=opn, exist=ex)
      if(opn) then
         call cerrormsg(mskel, 1)
         call cerrormsg(' already opened: starange', 0)
      elseif(ex .and.
     *      append(1:klena(append)) .eq. 'append' ) then
         open(mdev, file=mskel, form='unformatted',status='old')
         call cerrormsg('skeleton node info. will be appended', 1)
!           skip to the end of file
         do while( .true. )
            read(mdev, end=100)
         enddo
      else
         if(ex .and.
     *      append(1:klena(append)) .ne. 'append' ) then
            call cerrormsg(
     *           'Old skeleton node info. file exists', 1)
            call cerrormsg(
     *           'but node info. will NOT be appended', 1)
         endif
         open(mdev, file=mskel(1:klena(mskel)), form='unformatted',
     *       status='unknown')
      endif

 100  continue

      open(wdev, file=wskel(1:klena(wskel)), form='unformatted',
     *      status='unknown' )

      accepted = 0   ! counter;  accepted as skeleton 
!
      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"
#include "Ztrack.h"
#include "Ztrackv.h"
#include "Ztrackp.h"
#include "Zobs.h"
#include "Zobsp.h"
#include "Zobsv.h"
#include "Zcode.h"
      


      type(track)::incident
      type(coord)::angle


      integer  EventNo
      integer i, j
      integer seed(2)
      real*8 svEasWait, svEthin(4), kepn

      
      np = 0
      call cqincident( incident, angle)
      kepn = incident.p.fm.p(4)
      if( incident.p.code .eq. kgnuc ) then
         kepn = kepn/ incident.p.subcode
      endif
      ethresh = kepn * waitratio

      sveaswait = easwait       ! for safety save
      svethin(1) = ethin(1)           ! //
      svethin(2) = ethin(2)           ! //
      svethin(3) = ethin(3)           ! //
      svethin(4) = ethin(4)           ! //
      call csetemin(generate2, keminobs2(1), cutneg, cuteg)
      easwait = sveaswait       ! restore
      ethin(1) = svethin(1)
      ethin(2) = svethin(2)
      ethin(3) = svethin(3)
      ethin(4) = svethin(4)

      rewind  wdev

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

      eventno = eventno + 1     
      do i = 1, noofsites
         write(*, 1000)
     *    sngl(obssites(i).pos.depth),
     *    eventno, 
     *    incident.p.code,
     *    incident.p.subcode,
     *    incident.p.charge, 
     *    incident.p.fm.e,
     *    -angle.r(1),
     *    -angle.r(2),
     *    -angle.r(3)
 1000       format(f10.3,i9,3i4,e15.5,3(1x,f12.8))
      enddo


!     ===================================================
      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).
!
!     ===================================================
      common /testcos/sumg, ng(20), nth, eventno
      real*8 sumg
      integer ng, nth, EventNo

!     ===================================================
      if( id .eq. 2 .and. atrack.p.code .ne. kneumu .and.
     *   atrack.p.code .ne. kneue) then
         np = np + 1
         if( np .gt. npmax) then
            call cerrormsg(
     *   '# of particles >NpMax in observation', 0)
         endif
         o(np).where = atrack.where
         o(np).code = atrack.p.code
         o(np).subcode = atrack.p.subcode
         o(np).charge = atrack.p.charge
         o(np).atime = atrack.t
         o(np).erg = atrack.p.fm.p(4)
         o(np).mass = atrack.p.mass
         o(np).x = atrack.pos.xyz.r(1)
         o(np).y = atrack.pos.xyz.r(2)
         o(np).wx =atrack.vec.w.r(1)
         o(np).wy =atrack.vec.w.r(2)
         o(np).wz =atrack.vec.w.r(3)
         o(np).zenith = atrack.vec.coszenith
!     ===================================================
       if( o(np).code .le. 6 .and. o(np).code .ne. 3 ) then
         write(*, 959) 
     *  o(np).where,   
     *  o(np).code,   
     *  o(np).charge,  
     *  ( o(np).erg - o(np).mass ),   
     *  ( o(np).x ), ( o(np).y ) ,  
     *  ( o(np).wx ),  
     *  ( o(np).wy ),  
     *  ( o(np).wz ),  
     *  ( o(np).zenith ) 
 959    format(3i3, g14.3,2f16.6,4(1x,f12.8))
      endif

!     ===================================================
      endif
!     ===================================================

!     ===================================================
      end

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

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

      integer i
      integer seed(2)
      real  sumeg, sumeh
      logical memorize

      integer ng, nh
      
      ng = 0
      nh = 0

      sumeg = 0
      sumeh =0
!         count sum Eg etc
      do i = 1, np
         if(o(i).where .eq. where) then
            if(o(i).code .le. kelec) then
               ng = ng + 1
               sumeg = sumeg + o(i).erg
            elseif( ( o(i).code .ge. kpion .and.
     *                o(i).code .le. knuc )
     *             .or.
     *               o(i).code .eq. kgnuc) then
               nh = nh + 1
               sumeh = sumeh + o(i).erg
            endif
         endif
      enddo

!     ===================================================
      memorize =(ng .ge. ngmin .and. sumeg .ge. sumegmin) .or.
     *  ( nh .ge. nhmin .and. sumeh .ge. sumehmin)

!/////
      write(0,*) ' memo=', memorize,
     *   ' ng=',ng, ' seg=',sumeg, ' nh=',nh,' seh=',sumeh
!//////


!     ===================================================
      if(memorize) then
         accepted = accepted + 1
!c         call cwriteSeed    
!             !  SeedFile  (seed will be written by cmemorize)
!          flag for end of 1 event on working disk
         write(wdev)  -1, p
         call cmemorize(wdev, mdev)     !  reocord this event
      endif
      end


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


!     =========================================
      call cerrormsg('++++++++++++', 1)
      write(msg, '(i8, a)') accepted,
     *       ' events are memorized as skeleton'
      call cerrormsg(msg, 1)
      call cerrormsg('their  seeds are also memorized', 1)
      call cerrormsg('-----------',1)
      close(mdev)
      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
#include  "Zprivate.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Zcode.h"
!  #include  "Ztrackp.h"
      integer never             ! input & output
      integer nlow
! 
!    If Job = 'newskel', input  "never" may be < 0,  and MovedTrack
!    may not be an electron.
!      never = -1 :  MovedTrack K.E becomes < KEminObs due to
!                    energy loss during traveling.
!      never = -2 :  The same as above, but the particle crosses an
!                    observation depth, so the calling point to this
!                    routine is different from the never = -1 case.
!      never = -3 :  K.E >= KEminiObs.  The ptcl is observed at the
!                    current deepest Obs. level. but at the flesh time
!                    the deepest level will be more deep so that
!                    this must be memorized.
!
!         For above cases, the product is only MovedTrack and 
!         no particle is in PWork.
!  Otherwise,
!   MovedTrack is the electron that made interaction
!   Pwork contains produced particles (normally gamma, but mayb  e).
!   Nproduced has the number of particles in Pwork (normally 1)
!   IntInfArray(ProcessNo) contains the type of interaction
!        IntInfArray(ProcessNo).process will have one of
!       'brems', 'mscat', 'bscat'  'anihi' or 'mbrem'
!     
      if(never .lt. 0) then
         nproduced = 1
!           IBM user must modify next
         pwork(1) = movedtrack.p
      endif

!         high energy parent at node might be used
!        for hybrid AS generation if it is in some
!        energy region.
      if( movedtrack.p.code .eq. kelec ) then
         if( movedtrack.asflag .eq. 0 ) then
            if( movedtrack.p.fm.p(4) .lt. ethresh ) then
               movedtrack.asflag = -1
            endif
         endif
      endif

      call cmemonode(wdev, never)
      if(movedtrack.asflag .eq. -1) then
!          decendent should not be used to  generae A.S
         movedtrack.asflag  = -2
      endif

 10   continue
      never = 0
      end

!     ********************* this is the hook called when
!       a gamma ray made an interaction.
!
      subroutine chookgint(never)
      implicit none
#include  "Zprivate.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 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
!         IntInfArray(ProcessNo).process will have one of
!        'pair', 'comp', 'photoe' 'photop' 'mpair'
      call cmemonode(wdev, 1)
      never = 0
      end

!     ********************* this is the hook called when
!       non e-g particle made an interaction.
!
      subroutine chooknepint(never)
      implicit none
#include  "Zprivate.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
!

      call cmemonode(wdev, 1)
      never = 0
      end

      subroutine cmemonode(dev, flag)
      implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include  "Zprivate.h"
!
      integer dev  ! input. fortran logical dev. number for data write
      integer flag ! input. If this is -3, high energy particles must
                   !       also be memorized.
!
!
!       memorize nodal information at interaction.
!     
!        \
!         \  projectile    = MovedTrack
!         /|\  produced particles.: Pwork(i), i=1~Nproduced
!        / | \
!
!   From projectile, the track information is extracted
!   From produced particles, only those of K.E< KEminObs is
!   extracted and information needed for further tracking is 
!   obtaned and memorized. The position information is common
!   to all the children.
!
!   Track   information;   pos =
!                                xyz = r(1~3), sys
!                                radiallen, depth, height, colheight
!                            t
!                          vec  =
!                                 w = r(1~3), sys
!                                 coszenith
!                          wgt
!                         where
!                        asflag 
!
!      Among these, we don't memorize 
!         sys which is always 'xyz'
!       radiallen: can be computed from height
!         vec; children knows  their direction
!         wgt: normally not needed, it should be 1 for skeleton making
!              So thinsamping must not be used when making skeleton.
!      asflag: should be always F, (assume for skeleton making, hybrid
!              air shower is not requested)       
!
!   write  track info. of projectile
!
      integer i, nlow
      real*8 ke

!
!
      nlow = 0 
      do i = 1, nproduced
         ke = pwork(i).fm.p(4) - pwork(i).mass 
         if( ( pwork(i).code .le. kelec .and. 
     *       ke .ge.  cuteg ) .or. 
     *       ( pwork(i).code .gt. kelec .and.          
     *       ke .ge.  cutneg ) ) then
!              count low energy ptcls
            if(flag .ne. -3) then
               if( ke .lt.  keminobs(1)) then
                  nlow = nlow + 1
               endif
            else
!               all ptcl must be memorized
               nlow = nlow + 1
            endif
         endif
      enddo
!         next is harmful for B-approx air shower
!         because, parent is not memorized and cannot
!         generate AS at flesh time at proper positon
!         asflag == -1 must apear in flesh
!      if(nlow .eq. 0 )  return !  *************
!       if there is no return here, the size of skeleton data
!       becomes 3.2 times larger
      if(nlow .eq. 0  .and. movedtrack.asflag .ne. -1) return


      p.posx = movedtrack.pos.xyz.r(1)
      p.posy = movedtrack.pos.xyz.r(2)
      p.posz = movedtrack.pos.xyz.r(3)
      p.depth = movedtrack.pos.depth
      p.height = movedtrack.pos.height
      if( movedtrack.pos.colheight .gt. 1.e36) then
         p.colheight = 1.e36    ! no col. yet.
      else
         p.colheight =  movedtrack.pos.colheight
      endif
      p.atime = movedtrack.t
      p.where = movedtrack.where
      p.coszenith = movedtrack.vec.coszenith
      p.code = movedtrack.p.code
      p.erg  = movedtrack.p.fm.p(4)
      p.asflag = movedtrack.asflag


!c
!c     *           p.posx, p.posy, p.posz, p.depth, p.height, 
!c     *           p.colHeight, p.atime, p.where
!c
!       write particle info
!           p(1~4)
!           mass
!           code
!           subcode
!           charge
! 


      write(dev) nlow, p


      do i = 1, nproduced
         ke = pwork(i).fm.p(4) - pwork(i).mass 
         if( ( pwork(i).code .le. kelec .and. 
     *       ke .ge.  cuteg ) .or. 
     *       ( pwork(i).code .gt. kelec .and.          
     *       ke .ge.  cutneg ) ) then

            if(flag .eq. -3 .or.  ke .lt. keminobs(1)) then
               c.code = pwork(i).code
               c.subcode =  pwork(i).subcode 
               c.charge =  pwork(i).charge
               c.fm(1) = pwork(i).fm.p(1)
               c.fm(2) = pwork(i).fm.p(2)
               c.fm(3) = pwork(i).fm.p(3)
               c.fm(4) = pwork(i).fm.p(4)
               c.mass = pwork(i).mass

               write(dev) c
            endif
         endif
      enddo
      end

!        
      subroutine cmemorize(from, to)
      implicit none
#include "Zmaxdef.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#ifdef  DOMPI       
#include "mpif.h"
#include  "Zmpi.h"
#else
       integer::mpirank = -1
#endif      
      integer from      !  working disk file 
      integer to        !  permanent disk file where skeleton is sotered




      integer num, cumnum, irevent(2)
!      record /track/ incident
!      record /coord/ angle


      rewind from
!          need not memorize, can be generated at flesh
!      call cqIncident(incident, angle)      

      call cqeventno(num, cumnum)

      call cqinirn(irevent)   ! seed of the event

      write(to) cumnum, num, irevent, 
#if defined (KEKB) || defined (KEKA)
#include "ZavoidUnionMap.h"      
#else
     *  zfirst
#endif
      call cputhes(to)   ! put high energy showers.
      write(0,*) 'rank=',mpirank,
     * ' first Z=',zfirst.pos.depth*0.1, ' g/cm2',
     *  zfirst.pos.height, ' m ir=', irevent 
      
      call cputnodinfo(from, to)  ! put nordal info.

      end

      subroutine cputhes(to)
      implicit none
#include "Zprivate.h"
      integer to
!
!        write high energy sekeleton data into 'to'      
!
      integer i

      write(to) np

      do i = 1, np
         write(to) o(i)
      enddo

      end

      subroutine cputnodinfo(from, to)
      implicit none
#include "Zprivate.h"

      integer from, to
      
      integer i, nlow

      nlow = 1
      do while ( nlow .ge. 0 )
         read(from) nlow,  p
         write(to)  nlow,  p
         do i = 1, nlow
            read(from) c
            write(to) c
         enddo
      enddo

      end