cbeginRun.f

Go to the documentation of this file.

!        Initialize simulation.  
#if  defined (KEKB) || defined (KEKA)
#define DOMPI
#endif
!
      subroutine cbeginrun
      implicit none
#include  "Zmanagerp.h"
#include  "Zmanager.h"
#include  "Zelemagp.h"
#include  "Zevhnp.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include  "Zincidentp.h" 
#include  "Zprimary.h"
#include  "Zprimaryv.h"
#include  "Zcondc.h"
#include  "Zobs.h"
#include  "Zobsp.h"
#if defined (DOMPI)
#include  "mpif.h"
#include  "Zmpi.h"
      integer intdata
#endif
      character*16 temp
      integer jold, icon
      real*8 s1
!!!!!!!!!!!!
! after reading input namleist parameter, program flow
!     comes here.  the work dependent on LatObsSite and
!     LongitObsSite  must not be done before ciniAtmos,
!     if the NRL atmosphere model is to be used by specifing
!     AtmosFile="...." in which case, (lat, long) is given
!     in that file. So (lat,long) pair is forced to replace
!     LatitOfSite and LongitObsSite values.
      
      if(cont) then
!           restore status at the end of previous run
         call crestorestatus
         cont =.true.
      elseif(job .eq. 'flesh'  .or. job .eq. 'newflesh') then
         temp = job
         call cgetskelfile 
         job = temp    !  keet as now
      elseif(job .eq. 'skeleton' .or. job .eq. 'newskel' ) then
!            to be safe, parameters are written in a file first and read
!            again ; this will avoid the possible difference of internal
!            parameter values read from original param and SkeletonParam
!
#ifndef DOMPI
         temp = job
         call copennlfw(tempdev, skeletonfile, icon)
         if(icon .ne. 0) then
            call cerrormsg(skeletonfile, 1)
            call cerrormsg('File shown above cannot be opened',0)
         endif
!           at flesh time, you don't need to rewrite Job
         if(job .eq. 'newskel') then
            job = 'newflesh'
         else
            job = 'flesh'
         endif
         call cwriteparam(tempdev, 1)
         close(tempdev)
!           macIFC/MACOSX cannot put correct namelist char data 
!           (' is missing) so cannot read SkeletonFile
!           avoid read it. the user must do the equiv.  by hand 
!            call copenNLf(TempDev, SkeletonFile, icon)
!            call creadParam(TempDev)
!            close(TempDev)
!            restore Job
         job = temp
#endif
      endif
!        reset parameters
      if(borderheighth .eq. 0.0) then
         borderheighth = heightofinj + 1.0d0    ! from v6.10
      endif

      eventsintherun = 0      ! moved from ciniTracking0. 2001.09/30

      if(desteventno(2) .eq. 0) then
         desteventno(2) = desteventno(1)
      elseif(desteventno(1) .lt. 0) then
         desteventno(2) = -abs(desteventno(2))
      endif

!         to lower case letters.
      temp = generate
      call c2lowercase(temp, generate)
      temp = job
      call c2lowercase(temp, job)

      call rndsw(jold, 1)       ! specify random number generator 1.

!  Now (v7.651) next random number init is placed before cixsec
!     otherwise qgsjetII-03 generates the same event every run,
!     (ciniRN was placed inside cwhatJob bef. v7.651)
      call cinirn

!           set Target and do related xsec business; moved before
!           cintModels (for qgsjet dummy ptcl generation in cQGSini)
      call cixsec
      call chowmcs  ! MCSmodel is analysed ...
      call cinimcs  ! init for MCS and set doNewMCS (t/f)

!       cintModels probably reads disk file(s) internally, 
!       we may better to avoid simultaneous access to that disk 
!       by different ranks; reading  starts from rank0, rank1...
#if defined DOMPI
      if(mpirank .eq.  0) then
         call cintmodels('cosmos') ! analysis of interaction models and init.
         if(mpisize .gt. 1) then
            call mpi_send(mpirank, 1, mpi_integer, 1, 1,
     *           mpi_comm_world, mpierr)
         endif
      else
         call mpi_recv(intdata, 1, mpi_integer, mpirank-1, 1,
     *      mpi_comm_world, mpistat, mpierr)
         call cintmodels('cosmos') ! analysis of interaction models and init.
         if( mpirank .lt. mpisize-1 ) then
            call mpi_send(mpirank, 1, mpi_integer, mpirank+1, 1,
     *       mpi_comm_world, mpierr)
         endif
      endif
#else
      call cintmodels('cosmos') ! analysis of interaction models and init.
#endif
       ! init for Atmosphere
      call ciniatmos
      
!         init for geomag
      call crdgeomag(geomagfile, yearofgeomag)
!        init for LPM effect energy sampling.
      s1 = (targetatomicn**(1./3.d0)/183.d0)**2
      call csetlpmcnst(s1, log(s1), 1.d-4, x0)
!        initialize  observation  
      call cinitobs
!        init for   primary sampling
      call cinisprim(primaryfile)

      if(cutofffile .ne. ' ') then
         call crigcut0(cutofffile) ! read cutoff talbe and init.
      endif
!       init for muon interaction routines; specific for Air.
      call crdmutab        ! set various consts for mu int.
      call csetmu(targetatomicn, targetmassn)
      fromepics = .false.  ! muon interaction routines for Air are
                           ! inside Cosmos 
      call cinisprimang    !  this is in csPrimaAgn.f in Tracking dir.
!        check job
      call cwhatjob
!          this is moved here; before v6.10 it was before cwhatJob
      call cinitracking0  !  init for  tracking for all events
!           init for knockon process; Knockon below is not used now.
      if(knockonratio .lt. 1.d0) then
!              ********* Knockon is not used now ******
         if(job .eq. 'newskel') then
            call cdedxeleci(keminobs(1)*knockonratio, knockon)
         elseif(keminobs2(1)*knockonratio .gt. 0.) then
!                 this must come after cwhatJob, since KEminObs2 must
!           be fixed.  For skeleton-flesh job, KnockOnRatio should be
!           small enough so that KEminObs2*KnockOnRatio < KEminObs at
!           flesh time  

            call cdedxeleci(keminobs2(1)*knockonratio, knockon)
         else
            call cerrormsg('KnockOnRatio<1 and others mismatch', 0)
         endif
      else
!!!         call cdedxEleci(RecoilKineMinE, Knockon)  replaced by next few lines
!!             v7.643
         if( knockonratio == 1.0d0 .or. recoilkinemine == 0.) then
            recoilkinemine= keminobs(1)
         endif
         call cdedxeleci(recoilkinemine, knockon)
      endif
!         user hook
      call chookbgrun
      
      end

      subroutine cinitracking0
      implicit none
#include  "Zmanager.h"
#include  "Zmanagerp.h"
#include  "Zevhnp.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include  "Zprimary.h"
#include  "Zprimaryv.h"
#include  "Zobs.h"
#include  "Zobsv.h"
#include  "Zincidentp.h"
#include  "Zheavyv.h"

!

      integer i
       real*8 dstep 
!  
!
!
      observeas = index(generate, 'as')  .gt. 0 .or.
     *            index(generate, 'lat') .gt. 0 
      if(index(generate, 'qas') .gt. 0) then
         skipptclgen = 1  ! quick as generation for heavies.
      else
         skipptclgen = 0
      endif
      
!
!    d = min( dZ/2, geneal min) for e+/e- other charged.
!    general min = same as so far for e+/e-, p, heavy
!                  mu, pi, K, --> dE/Ek< 1%; dE=Ek/100.
!                  min=Ek/100/2e-3 g/cm2 = Ek/10/2e-3 kg/m2
!                     =Ek/2e-2= 50Ek 
!                  Ek=1  -->50 kg/m2  = 5g/cm2 ~5000cm ~50m
!                  Ek=0.1-->5 kg/m2=0.5g/cm2   ~ 5m
!                    0.01-->0.05kg/m2=0.05g/cm2 ~ 0.5m
!                    0.001-->                     0.5m  
!        
!             ////////                   
      do i = 1, noofsites
         if( i .eq. 1) then
            dstep = obssites(1)%pos%depth 
         else
            dstep = obssites(i)%pos%depth - obssites(i-1)%pos%depth
         endif
         if(dstep  .lt. 15. )  then
            stepcontrol=2
         else
            stepcontrol = dstep/25.0 
         endif
         maxstep(i) =   dstep/stepcontrol
      enddo
!         added v7.651
      maxstep(0) = maxstep(1)
      maxstep(noofsites+1) = maxstep(noofsites)

!          compute the offset point in 'xyz' system
!        the deepest detector origin + Offset is the point
!        to which the primary is directed.
!         offset in the detector system.
      offset%r(1) = 0.
      offset%r(2) = 0.
      offset%r(3) = offsetheight 
!        convert it to xyz system.
      call cdet2xyz(obssites(noofsites)%pos%xyz, offset, offset)
!        make it offset 
      do i= 1, 3
         offset%r(i) = offset%r(i) -
     *    obssites(noofsites)%pos%xyz%r(i)
      enddo

      if(eabsorb(1) .ne. 0) then
         if(eabsorb(2) .le. 0) then
            eabsorb(2) = noofsites
         elseif( eabsorb(2) .gt. noofsites) then
            call cerrormsg("Eabsorb(2) > NoOfSites", 0)
         endif
      endif
      end
!     ************
      subroutine  cinirn
!     initialize random # generator. This was a part of
!     cwhatJob before v7.651 but separated as ciniRN
!     and is called from before cixsec. 
!     ************
      implicit none
#include  "Zmanager.h"
#include  "Zmanagerp.h"
!
!
      character*190  msg
      real*8 u
      integer::i,   now(2)

!      move from ciniTracking0  
      refreshir = initrn(1) .lt. 0 .and. 
     *        ( job .ne. 'flesh' .and. job .ne. 'newflesh')

      if(initrn(1) .gt. 0 .and. initrn(2) .gt. 0 ) then
         call rnd1r(initrn)     ! init randeom number generator
!      *****************
      elseif(.not. refreshir .and. initrn(2) .lt. 0) then
         call cmkseed(0, now)   ! make seed using timer and hostname
         call rnd1r(now)
!           dummy use of 1000 times
         do i = 1, 1000
            call rndc(u)
         enddo
      endif
!          this is almost ok but later once more saved.
      call rnd1s(seedsave)
      end
!     ******************
      subroutine  cwhatjob
!     ************
      implicit none
#include  "Zmanager.h"
#include  "Zmanagerp.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include  "Zobs.h"
#include  "Zobsp.h"
#include  "Zobsv.h"
!   
      integer klena, icon
      character*8  uid 
      character*16 temp

      integer i
      character*190  msg

      if(keminobs(2) .ne. keminobs(1)) then
         write(0,*) ' KEminObs(2) is forced to be the same as'
         write(0,*) ' KEminObs(1)=',keminobs(1)
         keminobs(2)= keminobs(1)
      endif

      if(job .eq. ' ' .or. job .eq. 'skeleton' .or. 
     *   job .eq. 'newskel'  ) then
         if(job .ne. 'newskel') then
!             save present conditions
            do i = 1, 8
               keminobs2(i) = keminobs(i)
            enddo
            generate2 = generate
            EndLevel2 = endlevel
         elseif(job .eq. 'newskel') then
            if( keminobs2(1) .ge. keminobs(1) .and.
     *          endlevel2 .le. endlevel .and.
     *          index(generate2,'as') .eq. 0 .and.
     *          index(generate2,'lat') .eq. 0 ) then
               call cerrormsg(
     *          'Doing newskel job seems nonsense', 1)
               call cerrormsg(
     *          'Check Generate2, KEminObs2(1), EndLevel2',0)
            endif
         endif
         noofsites2 = noofsites    ! probably not needed
         if(job .eq. ' ') then
            if(seedfile .ne. ' ') then
!                  open seed file for output
               write(msg, *) 'opening SeedFile=',
     *         seedfile(1:klena(seedfile))
               call cerrormsg(msg, 1)
               call copenfw(seedfiledev, seedfile, icon)
               if(icon .ne. 0) then
                  call cerrormsg(seedfile, 1)
                  call cerrormsg('File shown above cannot be opened',0)
               endif
               if(cont) then
                  call cskiptoeof(seedfiledev)
               endif
            endif               
         elseif(job .eq. 'skeleton' .or. job .eq. 'newskel' ) then
            call
     *      cerrormsg('  ********** skeleton making **********', 1)
            write(msg, *)  '      Generate=', generate
            call cerrormsg(msg, 1)
!                 save skeleton inf. in skelotonFile file.
!             The file will be modified when the distjob command
!             processes Job = 'flesh' later. You need not modify
!             skeleton file if distjob is employed.
            if(.not. cont) then
               temp = job    ! save current Job
               if(job .eq. 'newskel') then
                  job = 'newflesh'
               else
                  job = 'flesh'
               endif
               call copennlfw(tempdev, skeletonfile, icon)
               if(icon .ne. 0) then
                  call cerrormsg(skeletonfile, 1)
                  call cerrormsg(
     *              'File shown above cannot be opened',0)
               endif
               call cwriteparam(tempdev, 1)
               close(tempdev)
               job = temp
            endif

!                open SeedFile 
            if(seedfile .eq. ' ') then
!                error. you need file; not needed actually
!               write(msg, *)
!     *         ' SeedFile must not be blank for skelton making'
!               call cerrorMsg(msg, 0)
            else
               write(msg, *) 'opening SeedFile=',
     *         seedfile(1:klena(seedfile))
               call cerrormsg(msg, 1)
               call copenfw(seedfiledev, seedfile, icon)
               if(icon .ne. 0) then
                  call cerrormsg(seedfile, 1)
                  call cerrormsg('File shown above cannot be opened',0)
               endif
               if(cont) then
                  call cskiptoeof(seedfiledev)
               endif
            endif
         endif
      elseif(job   .eq. 'flesh' .or. job .eq. 'newflesh') then
!             don't worry about KEminObs2 etc.  They have been read
!             from &Hparam
          call cerrormsg('  ********** fleshing job   *********', 1)
          if(job .eq. 'flesh') then
             if(endlevel .gt. endlevel2) then
!                to deeper detph than skeleton
                write(msg, *)
     *          ' fleshing will be done to deeper depth than'//
     *               ' skeleton making time'
                call cerrormsg(msg, 1)
                write(msg, *) ' No of old levels=', endlevel2,
     *               ' No of new levels=', endlevel
                call cerrormsg(msg, 1)
             elseif(endlevel .lt. endlevel2) then
                call cerrormsg('EndLevel must be >= skelton time', 0)
             endif
             write(msg, *)  '      Old Generate=', generate2
             call cerrormsg(msg, 1)
             write(msg, *)  '      New Generate=', generate
             call cerrormsg(msg, 1)
          else
             if(endlevel .lt. endlevel2) then
!                to deeper detph than skeleton
                write(msg, *)
     *          ' fleshing will be done to deeper depth than'//
     *           ' skeleton making time'
                call cerrormsg(msg, 1)
             endif
!
!              copy  old Generate, KEminObs2 to current value
!              they are future values at newskel time
             do i = 1, 8
                keminobs(i) = keminobs2(i)
             enddo
             generate = generate2
             EndLevel = endlevel2
          endif
!                open SeedFile
          if(seedfile .eq. ' ') then
!c             write(*, *) ' SeedFile must not be blank for flesh job'
!c             call cerrorMsg(msg, 0)
!               seed will be read from Mdev.
          else
             write(msg, *) 'opening SeedFFile=', 
     *        seedfile(1:klena(seedfile))
             call cerrormsg(msg, 1)
             call copenf(seedfiledev, seedfile, icon)
             if(icon .ne. 0) then
                call cerrormsg(seedfile, 0)
                call cerrormsg('File shown above seems missing',0)
             endif
          endif
      else
           write(msg,*) ' Job=',job, ' undefined'
           call cerrormsg(msg, 0)
      endif
      if((trace .gt. 0 .and. trace .lt. 60) .or. trace .gt. 100) then
!           defalut trace.  fix the dirctor
          if(tracedir .eq. ' ') then
             call cgetloginn(uid)
             tracedir = '/tmp/'//uid(1:klena(uid))
          endif
       endif
      end
!        **************************************** read cont job info
      subroutine crestorestatus
      implicit none
#include "Zmanagerp.h"
      integer icon

      call copennlf(tempdev, contfile,icon)
      if(icon .ne. 0) then
         call cerrormsg(contfile, 1)
         call cerrormsg('File shown above seems missing',0)
      endif
      call creadparam(tempdev)
      close(tempdev)
      end
      subroutine cgetskelfile
      implicit none
#include "Zmanagerp.h"
      integer icon
      call copennlf(tempdev, skeletonfile, icon)
      if(icon .ne. 0 ) then
         call cerrormsg(skeletonfile, 1)
         call cerrormsg('File shown above seems missing',0)
      endif
!          read skelton parameters for flesing
      call creadparam(tempdev)
      close(tempdev)
      end