chook.f File Reference

Go to the source code of this file.

Modules
module	modassize

Functions/Subroutines
subroutine	chookbgrun
integer function	csighandler (sig, code, context)
subroutine	chookbgevent
subroutine	chookobs (aTrack, id)
subroutine	chookenevent
subroutine	chookenrun
subroutine	chooktrace
subroutine	chookeint (never)
subroutine	chookgint (never)
subroutine	chooknepint (never)

Variables
integer, save	nlayers =0
integer, save	nev =0
integer, parameter	nbinr =2*4
integer, save	nthlayer =0
integer, parameter	nkind =6
integer, parameter	nbine =20*nkind
real(4), save	bine =1./(nbinE/nkind)
real(4), save	binr =1./(nbinR/4)
type(histogram2), dimension(nkind)	erhist
type(histogram2), dimension(nkind)	terhist
real(4), dimension(nkind), save	eminhist
integer, parameter	nemin =2
integer, parameter	nkindtran =7
real(8), dimension(nkindtran, nemin), save	emintran
real(8), dimension(:,:,:,:), allocatable, save	num
real(8), dimension(:,:,:), allocatable, save	pnnum
integer, parameter	fno =21

Function/Subroutine Documentation

chookbgevent()

subroutine chookbgevent

(

)

Definition at line 141 of file chook.f.

References cqinirn(), modassize::erhist, modassize::nkind, modassize::num, and modassize::pnnum.

      use modassize
      implicit none
!          ///////////
#include "Ztrack.h"
#include "Zobs.h"
#include "Zobsv.h"


      integer i

      integer seed(2)
!      write(*, *) ' bigin event generation'
      call cqinirn(seed)
!     write(*,*) ' seed=', seed
       !  clear the histogram area of this event
      do i=1, nkind
         call kwhistc2(erhist(i))
      enddo

!          clear the transtion counter of this event
      num(:,:,:,:) = 0.
      pnnum(:,:,:) = 0.

Here is the call graph for this function:

chookbgrun()

subroutine chookbgrun

(

)

Definition at line 52 of file chook.f.

References b, modassize::bine, modassize::binr, cprintobs(), cprintprim(), csighandler(), cwriteparam(), depth, modassize::eminhist, modassize::erhist, modassize::fno, kwhistso(), modassize::nbine, modassize::nbinr, modassize::nemin, modassize::nkind, modassize::nkindtran, modassize::nlayers, modassize::nthlayer, modassize::num, modassize::pnnum, and modassize::terhist.

      use modassize
      use modhistogram1
      implicit none
#include "Zmanagerp.h"
#include "Ztrack.h"
#include "Ztrackv.h"
#include "Zobs.h"
#include "Zobsp.h"
#include "Zobsv.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.

      integer:: i
      
#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)
      if(noofsites /= noofassites) then
          write(0,*)
     *    '# of Depth and ASdepth in this appli. must be the same'
          stop
      endif
      nlayers = noofsites
      if( any( asobssites(1:nlayers).pos.depth /=
     *         obssites(1:nlayers).pos.depth ) ) then
         write(0,*)
     *  'Depth and ASDepth must be the same in this appli.'
         stop
      endif
      if( nthlayer > nlayers ) then
         write(0,*)' nthlayer > # of given depths=',nlayers
         stop
      elseif( nthlayer == 0 ) then
         nthlayer = nlayers
      endif
         
         
      open(fno, file='try.hist', form='formatted')
      call kwhistso(1)   ! ascii output 
      do i =1,  nkind
         !  instance of each ptcl class: kind
         call kwhisti2(erhist(i),  
     *     eminhist(i), bine, nbine, b'00011',
     *        1.0,  binr, nbinr, b'00001' )
         call kwhisti2(terhist(i),  
     *     eminhist(i), bine, nbine, b'00011',
     *        1.0,  binr, nbinr, b'00001' )
      enddo
      !  clear the histogram area for total events
      do i=1, nkind
         call kwhistc2(terhist(i))
      enddo
      ! counter for transition of nkindTran 
      allocate( num(nlayers, nemin, nkindtran, 2) )

!     2 is for p,n (pbar,nbar included); specially take
      !   p and n tran. only down goingn
      allocate(pnnum(nlayers, nemin, 2) )

Here is the call graph for this function:

chookeint()

subroutine chookeint

(

integer

never

)

Definition at line 492 of file chook.f.

      use modassize
      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 electr.on 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
!       never = 2: same as 0
!             = 3: discard all the child of this (but not other ptcls)
!             = 4: discard the event
!
!        IntInfArray(ProcessNo).process will have one of
!       'brems', 'mscat', 'bscat', 'anihi' or 'mbrem'
!

chookenevent()

subroutine chookenevent

(

)

Definition at line 277 of file chook.f.

References depth, modassize::erhist, height, modassize::nemin, modassize::nev, modassize::nkind, modassize::nlayers, modassize::num, modassize::pnnum, and modassize::terhist.

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




      integer:: i, j

      nev= nev + 1
      do i = 1, nkind 
         if( nev == 1 ) then
            terhist(i) = erhist(i)
         else
            call kwhista2(terhist(i), erhist(i), terhist(i)) ! accumulation
         endif
      enddo

!----------transition: Num; Emin matter
      do i = 1, nlayers
         do j = nemin-1, 1, -1
            num(i, j, :, :) = num(i, j, :, :) + num(i,j+1,:,:)
            pnnum(i, j, :) =pnnum(i,j,1) + pnnum(i, j+1, :)
         enddo
      enddo
      
      if(observeas) then
!     electron size in B approx and size weighted age
          ! counter for transition of nkindTran 
          !   Num(Nlayers, nEmin, nkindTran, 2) 
         write(*,'(a)')
     *   "# id L  v dep     H    age   Nehyb        Ne1        Ne2"
         do i = 1, nlayers
            write(*, '(a, i2, f8.1, f8.0, f5.2, 1p, 3g12.5)')
     *            "Ne ", i, obssites(i).pos.depth/10.,
     *      obssites(i).pos.height, asobssites(i).age,
     *      asobssites(i).esize, num(i, 1, 2, 1), num(i, 2, 2, 1)
         enddo
      endif
!-----------n,p
      write(*,'(a)')
     *     "# id  L   v dep    H      Np1         Nn1 "//
     *     "      Np2         Nn2"
      do i = 1, nlayers
         write(*,'(a, i2, f8.1, f8.0,  1p, 4g12.5)')
     *        "Np/n ", i, obssites(i).pos.depth/10.,
     *        obssites(i).pos.height, pnnum(i, 1, 1), pnnum(i,1, 2),
     *        pnnum(i, 2, 2), pnnum(i, 2, 2)
      enddo
!-------------  g tran
      write(*,'(a)' )
     *     "# id L  v dep     H      Ngd1        Ngu1 "//
     *     "       Ngd2       Ngu2"
      do i = 1, nlayers
         write(*,'(a, i2, f8.1, f8.0,  1p, 4g12.5)')
     *        " Ng ", i, obssites(i).pos.depth/10.,
     *        obssites(i).pos.height,
     *        num(i, 1, 1, 1), num(i, 1, 1, 2),
     *        num(i, 2, 1, 1), num(i, 2, 1, 2)              
      enddo

!-------------  e tran
      write(*,'(a)' )
     *     "# id L  v dep     H      Ned1        Neu1 "//
     *     "       Ned2        Neu2"
      do i = 1, nlayers
         write(*,'(a, i2, f8.1, f8.0,  1p, 4g12.5)')
     *        " Ne ", i, obssites(i).pos.depth/10.,
     *        obssites(i).pos.height,
     *        num(i, 1, 2, 1), num(i, 1, 2, 2),
     *        num(i, 2, 2, 1), num(i, 2, 2, 2)              
      enddo

      !-------------  mu tran
      write(*,'(a)' )
     *     "# id L  v dep     H      Nmd1        Nmu1 "//
     *     "       Nmd2        Nmu2"
      do i = 1, nlayers
         write(*,'(a, i2, f8.1, f8.0,  1p, 4g12.5)')
     *        "Nmu ", i, obssites(i).pos.depth/10.,
     *        obssites(i).pos.height,
     *        num(i, 1, 3, 1), num(i, 1, 3, 2),
     *        num(i, 2, 3, 1), num(i, 2, 3, 2)              
      enddo

!-------------  had tran
      write(*,'(a)' )
     *     "# id L  v dep     H      Nhd1        Nhu1 "//
     *     "       Nhd2        Nhu2"
      do i = 1, nlayers
         write(*,'(a, i2, f8.1, f8.0,  1p, 4g12.5)')
     *        " Nh ", i, obssites(i).pos.depth/10.,
     *        obssites(i).pos.height,
     *        num(i, 1, 4, 1), num(i, 1, 4, 2),
     *        num(i, 2, 4, 1), num(i, 2, 4, 2)              
      enddo

      !-------------  nue tran
      write(*,'(a)' )
     *     "# id L  v dep     H      Nned1       Nneu1 "//
     *     "       Nned2        Nneu2"
      do i = 1, nlayers
         write(*,'(a, i2, f8.1, f8.0,  1p, 4g12.5)')
     *        "Nnue ", i, obssites(i).pos.depth/10.,
     *        obssites(i).pos.height,
     *        num(i, 1, 5, 1), num(i, 1, 5, 2),
     *        num(i, 2, 5, 1), num(i, 2, 5, 2)              
      enddo


!-------------  numu tran
      write(*,'(a)' )
     *     "# id L  v dep     H      Nnmd1       Nnmu1 "//
     *     "       Nnmd2        Nnmu2"
      do i = 1, nlayers
         write(*,'(a, i2, f8.1, f8.0,  1p, 4g12.5)')
     *        "Nnum ", i, obssites(i).pos.depth/10.,
     *        obssites(i).pos.height,
     *        num(i, 1, 6, 1), num(i, 1, 6, 2),
     *        num(i, 2, 6, 1), num(i, 2, 6, 2)              
      enddo


      !-------------  numu tran
      write(*,'(a)' )
     *     "# id L  v dep     H      Nod1       Nou1 "//
     *     "       Nod2        Nou2"

      do i = 1, nlayers
         write(*,'(a, i2, f8.1, f8.0,  1p, 4g12.5)')
     *        "Noth ", i, obssites(i).pos.depth/10.,
     *        obssites(i).pos.height,
     *        num(i, 1, 7, 1), num(i, 1, 7, 2),
     *        num(i, 2, 7, 1), num(i, 2, 7, 2)              
      enddo

!       separator
      write(*, *)

chookenrun()

subroutine chookenrun

(

)

Definition at line 427 of file chook.f.

      use modassize
      implicit none

chookgint()

subroutine chookgint

(

integer

never

)

Definition at line 525 of file chook.f.

      use modassize
      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
!       never = 2: same as 0
!             = 3: discard all the child of this (but not other ptcls)
!             = 4: discard the event

!         IntInfArray(ProcessNo).process will have one of
!        'pair', 'comp', 'photoe' ,'photop' or 'mpair'

!       

chooknepint()

subroutine chooknepint

(

integer

never

)

Definition at line 559 of file chook.f.

      use modassize
      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 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
!
!        default setting
      never = 1
!       never = 2: if proton is the current particle, give 2 
!                  if you make a hybrid air shower from that and
!                  want to discard that proton.
!             = 3: discard all the child of this (but not other ptcls)
!             = 4: discard the event

!
!        IntInfArray(ProcessNo).process  will have
!             'col' or 'decay'

chookobs()

subroutine chookobs	(	type(track)	aTrack,
		integer	id
	)

Definition at line 171 of file chook.f.

References charge, code, coszenith, modassize::emintran, modassize::erhist, mass, modassize::nemin, modassize::nthlayer, modassize::num, p, modassize::pnnum, wgt, x, xyz, and y.

!
!     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).
!
      use modassize
      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

      integer:: kind, updw, iemin, kindtran
      real(4)::ke,r
!
!     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).
!c                        !  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
!c     *  sngl(aTrack.vec.w.x),  ! direc. cos.x in the current detector system.
!c     *  sngl(aTrack.vec.w.y),  ! direc. cos.y
!c     *  sngl(aTrack.vec.w.z),  ! direc. cos.z
!     *  sngl(aTrack.vec.coszenith) ! cos of zenith angle
!c         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
         kind=atrack.p.code
         
         ke = atrack.p.fm.p(4) - atrack.p.mass
         do iemin = nemin, 1, -1
            if( ke >= emintran(min(kind,9),iemin) ) goto 10
         enddo
         iemin= 0
 10      continue
         if( iemin > 0  ) then
            if( kind == 6 .and. atrack.vec.coszenith >0.) then
               ! p,n (or pbar,nbar)
               if( atrack.p.charge == 0 ) then
                  pnnum(atrack.where, iemin, 2) =
     *                 pnnum(atrack.where, iemin, 2) + atrack.wgt
               else
                  pnnum(atrack.where, iemin, 1) =
     *                 pnnum(atrack.where, iemin, 1) + atrack.wgt
               endif
            endif
            if( kind <= 3 ) then
               kindtran= kind
            elseif(kind >= 4 .and. kind <=6 .or. kind == 9 ) then
               kindtran = 4         !  pi,K, n,A
            elseif( kind == 7 ) then
               kindtran = 5  ! 7,8  neue
            elseif( kind ==8 ) then
               kindtran = 6    ! 8  neumu
            elseif(kind > 8) then
               kindtran = 7     ! others
            endif
            if(atrack.vec.coszenith > 0.) then
               updw = 1
            else
               updw = 2
            endif

            num(atrack.where, iemin, kindtran,  updw) =
     *       num(atrack.where, iemin, kindtran, updw) + atrack.wgt
         endif
!-------
         if( atrack.where == nthlayer ) then
            if(atrack.vec.coszenith > 0.) then
               if(kindtran < 7 ) then
                 r = sqrt(atrack.pos.xyz.x**2 + atrack.pos.xyz.y**2)
                 call kwhist2( erhist(kindtran),  ke, r,
     *                  sngl(atrack.wgt) )
               endif
            endif
         endif
      endif

chooktrace()

subroutine chooktrace

(

)

Definition at line 439 of file chook.f.

      use modassize
      implicit none

#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include  "Zobs.h"
#include  "Zobsv.h"
#include  "Zcode.h"
       type(coord)::f
       type(coord)::t
       logical compress/.true./
       save compress


!       real*4 h1,  h2
!
!    If trace > 100.
!    every time a particle is moved in the atmosphere, this routine is called.
!
!         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
!      This example here is to put only muons with the same format
!       the standard trace infomation.  Trace-100 is treated as
!       the standard Trace value.
!     'compress' is to indicate compressed output.
!
!      if(TrackBefMove.p.code .eq. kmuon) then
!         call ccoordForTr(Trace-100, f, t)   ! convert coordinate
!         call cwrtTrInfo(compress, f, t)    ! write trace data
!      endif

csighandler()

integer function csighandler	(	integer	sig,
		integer	code,
		integer, dimension(5)	context
	)

Definition at line 127 of file chook.f.

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