interface.f

Go to the documentation of this file.

#include "csoftenPiK.f"
#include "cutptcl.f"
      subroutine xbgrun
      implicit none
#include  "Zmaxdef.h"
#include "Zglobalc.h"
#include "Zmanagerp.h"
#include "Ztrack.h"
#include "Ztrackp.h"
#include "Ztrackv.h"
#include "Zcode.h"
#include "Zheavyp.h"
#include "Zobs.h"
#include "Zobsp.h"
#include "Zobsv.h"
#include  "Zstackv.h"
#include "Zprivate.h"
#include "Zprivate2.h"
#include "../Zabsorb.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

!      type(track):: inci
!////////////
      type(coord):: pdir, cdir
!////////////
      type(coord):: tetafai
      
      character*128 input
      character*64 dirstr
      real sr, dr, tempr
      integer i, j, k, m,  icon
      integer ansites
      save ansites
      integer iij, code
      integer i1, i2, ic
      integer ir,  ifai, l, ridx, faiidx
      real*8  E0, cosz
      real*8  fai0, fai, sint
      real*8  delta  
      integer reducedTime
      integer NN
      integer klena
      integer w2hl(maxnoofsites)
      real*8 r, Eloss, rinmu, cosang
      real*8 dedt, dedtF, rho, dist, disto, BinFai
      real*8 aa
      real*8 wx, wy, wz, temp
      real   za
      real  de, Ek, f, molu
      real  dt, tmin 
      real*8  cvh2den
      data binfai/30./
      integer ldep
!     integer ndummy
      character*9 ptcln(4)
      data ptcln/"Photons", "Electrons","Muons", "hadron"/
      character*9 ptcl2(3)
      data ptcl2/"Electrons", "Muons","All"/
      real power(4)
      integer nstr
      data power/1.,1.,1., 1./
      real  power2(3)
      data power2/1.,1.,1./
      character*128 title
      character*96 evid(nsites)
      save evid
      real*8 cog, cog2, sumne,  obstimes, Savederg(5)
      real*8 firstcdepth, dd
      logical dosort
      real*4  wt, stime
      real*8 sumEbydEdx, sumEbyDeath,sumEbyDeathNeu,sumEbyDeathNut
      real*8 sumEbyDeathE, sumEbyDeathG, sumEbyDeathMuPiK, 
     *      sumebydeathp, sumebydeatho
      real*8 sumEcrash, sumEspace
      real*8 sumAll, sumdEinAir, sumMissing, sumUncertain
      integer vn/2/ ! version number for the fnoB output
      save 
!/////////////
      real*8 pabs, rcore, sina, cs, sn,  cf, mom(3), Ek8, u
!/////////////

!     ***********************
      include "interface1.h"
!     *********************


      do i = 1, nsites
         w2hl(i) = 0
      enddo

      do i = 1, nsites
!             histdep(i) is the layer number
         if(histdep(i) .eq. 0)  exit
         ansites = i
         w2hl( histdep(i) ) = i
      enddo

      r=rmin
      dr = 10.**bin 

      do i = 1, nrbin
!            center of the bin:   
         rbin(i) = r
         r = r* dr
      enddo
#if defined (MACOSX)
#else
!            specify bin or ascii output
      call kwhistso( binw )
#endif

      return      
!     *********************************** 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.
!     *
      entry xbgevent

      call cqincident(inci, angle)
      e0 = inci.p.fm.p(4)
      if(inci.p.code .eq. kmuon) then
         call csetmuonpol(1.0d0)
      endif
      cosz = -angle.r(3)
      fai0 = atan2(-angle.r(2), -angle.r(1))*todeg
      sint = sqrt(1.0-cosz**2)

      if(inci.p.code .eq. 9) then
         nn= inci.p.subcode
      elseif(inci.p.code .eq. 1) then
         nn=0
      else
         nn=1
      endif

      write(0,'("i ",  i6,  i4, g13.4,3f11.7,f7.1)')
     *   eventno, inci.p.code,
     *   inci.p.fm.e,  -angle.r(1),  -angle.r(2), -angle.r(3)
       write(0,'(a, 1p, 6g15.5)')
     *    '### ', detxaxis.r(1:3), detzaxis.r(1:3)
      do i = 1, noofsites
         sumeloss(i) = 0.
         do j = 1, 4
            ng(i) = 0.
            ne(i) = 0.
            nmu(i) = 0.
            nhad(i) = 0.
         enddo
         decent(i) = 0.
         do ifai = 1, nfai
            do ir= 1, nrbin
               derfai(ir, ifai, i) = 0.
!               do j = 1, 4
!                  nrfaiAll(ir, ifai, j, i) = 0.
!               enddo
            enddo
         enddo
      enddo


!          estimate time minimum and time bin for eeach web sector
      
      do i = 1, ansites
         ldep = histdep(i)
         call cmintime2websec(obssites(ldep).pos.xyz,
     *        ldep, i,  webmin )
      enddo

#if defined (MACOSX)
#else
!     histogram: instanciate
!           t spectrum at each web sector
      if(tkrtspec) then
         do i = 1, ansites
            do j = 1, 4
!                at  center
               call kwhisti( tspec0(j, i),
     *          -5., 0.05, 200, b'00000')
               call kwhistai( tspec0(j, i),
     *         "Arrival time dist. of "//ptcln(j)//" at center",
     *         "t", "ptcls", .false., 0., "time", "ns")
!                   clear
               call kwhistc(tspec0(j, i))

               do ir=1, nrbin
                  do ifai=1, nfai
                     if(reducedtime .eq. 1) then
                        tmin = webmin(ir, 7, i)
                     else
                        tmin = webmin(ir, ifai,i)
                     endif
                     dt = 0.01*10.0**(bin*(ir-1))*100. ! approx core distnace m
                     dt = dt**0.75*1.e9/3.0e8/100. ! if sqrt 1m-->0.03 ns 10 m-->0.15 ns
                                       ! 100m 1ns 1km 5ns   4km 10ns
                                       ! dt**0.65  makes larger bin at large distance (<=x2)
                     if(j .eq. 4) dt=dt*10.0*ir/35.0 ! for delayed hadrons
                     dt= max(dt, 0.02)

                     call kwhisti( tspec(j, ir, ifai, i),
     *                    tmin,  dt, 2000,   b'00000')

                     call kwhistai( tspec(j, ir, ifai,  i),
     *                "Arrival time of "//ptcln(j)//" at (r,fai)",
     *                "rt", "ptcls", .false., 0., "time", "ns")
!                     clear 
                     call kwhistc(tspec(j, ir, ifai, i))
                  enddo
               enddo
            enddo
         enddo
      endif

!            lateral in each fai bin
      if(tkarspec) then

         do i = 1, ansites
            do j = 1, 4         ! g,e,mu,h
               do ifai = 1, nfai
                  call kwhisti(rspec(j, ifai, i),
     *                 rmin, bin, nrbin,  b'00011' )
                  call kwhistai(rspec(j, ifai, i), 
     *            "Lateral Dist. of "//ptcln(j)//" at  diff. azimuth",
     *            "ar", "ptcls", .true.,  power(j),   "r", "m.u")
!                     clear
                  call  kwhistc( rspec(j, ifai, i) )
               enddo
            enddo
         enddo
      endif
#endif

      obstimes = 0.

      return
!     ***************
      entry xobs(atrack, id)
!
!     For id =2, you need not output the z value, because it is always
!     0 (within the computational accuracy).
!

      obstimes = obstimes + 1.d0
      if(mod(obstimes, 100000.d0) .eq. 0. ) then 
         dosort=.false.
         do i = 1, min(4,stack_pos)
            if(stack(i).p.fm.p(4) .ne. savederg(i)) then
               savederg(i)=stack(i).p.fm.p(4) 
               dosort=.true.
            endif
         enddo
         if(dosort) then
            call csortstack
         endif
         write(0, *) ' obstimes=', obstimes, ' ptclE=',atrack.p.fm.p(4)
         do i = 1, min(4,stack_pos)
            write(0,*)' stack tops=', stack(i).p.fm.p(4)
         enddo
      endif
!     ***************
      code = atrack.p.code
      ldep =  atrack.where
!     ************
      if(id .eq. 2 .and. code .le. 6 ) then  ! neglect rare ptcls
         wz = atrack.vec.w.r(3) ! downgoing < 0
         if(wz .gt. 0) return
         wz = -wz
         r = sqrt( atrack.pos.xyz.x**2 +
     *                 atrack.pos.xyz.y**2 )
         molu = obssites(ldep).mu
         rinmu =r/molu
         sr = rinmu   ! single precision
         ridx = (log10( rinmu/rmin )/bin +0.5) +1 

         ek = atrack.p.fm.p(4) -atrack.p.mass
         wt = atrack.wgt  ! wt is single
         if(code .eq. kphoton) then
            ng(ldep) = ng(ldep) + atrack.wgt
         elseif(code .eq. kelec) then
            ne(ldep) = ne(ldep) + atrack.wgt
         elseif(code .eq. kmuon) then
            nmu(ldep) = nmu(ldep) + atrack.wgt
         elseif(code .le. 6) then
            nhad(ldep) = nhad(ldep) +atrack.wgt
         endif
!            ---------- compute energy loss 
         if(atrack.p.charge .ne. 0  ) then
!             -----------------
!                      /|    |
!                     / |   1g/cm2 
!                    /A |    |
!            -------------------
!                  / ptcl direction  
!         get energy loss when aTrack goes some distance
!         of which vertical gramage is 1g/cm2.
!         Gramage the particle travel is 
!         1/cos where cos is the cos of angle (i.e, A if Fig)
!          in the detctor system.
!         1g/cm^2 = 10-3kg/10-4 m^2 =10 kg/m^2.
!         To travel  1 g/cm^2, the ptcl must
!         run dist kg/m^2
            if(abs(wz) .gt. 1.d-2) then
               dist =10./wz    ! in kg/m2/(g/cm2)
            else
!                 for safety
               dist =1000.
            endif

            call cqelossrate(dedt,dedtf) !  loss rate GeV/(kg/m^2)
!                     dedtF is the full eloss ; dedt is the restricted
!                     loss.  We may better use full here.
!                       energy in 1 g/cm2 of vertical direction
            eloss =min( real(dedtF*dist), Ek)   !  GeV/(g/cm2)
            eloss = eloss*atrack.wgt !  GeV/(g/cm2)
            sumeloss(ldep)=sumeloss(ldep) + eloss
         else
            eloss=0.
         endif

         if(code .ge. 4) code=4
         if( atrack.p.charge .ne. 0   .or. 
     *       w2hl(ldep) .gt. 0 ) then 
!                  fai
!              fai is  in    -15 to 345  (for dfai=30.)                                     
            aa=atan2(atrack.pos.xyz.y, atrack.pos.xyz.x)*
     *        todeg -fai0
            fai = aa/todeg
            aa= mod(aa + 360.d0,   360.d0)
            if(aa .gt. (360.d0-dfai/2.0d0)) aa= aa-360.d0
            faiidx=(aa+dfai/2.0d0) /dfai + 1
            if(ridx .ge. 1 .and. ridx .le. nrbin) then
               derfai(ridx, faiidx, ldep) = derfai(ridx, faiidx, ldep)  
     *             +  eloss         
            elseif(ridx .le. 0) then
               decent(ldep) = decent(ldep) + eloss
            endif
! 
!  do following   for specified histo layers (typically only 1 layer)
!
#if defined  (MACOSX)
#else
            if( w2hl(ldep)  .gt. 0 ) then
               i = w2hl(ldep)
               if(tkarspec) then
                  call kwhist( rspec(code, faiidx,  i), 
     *            sr, wt)
               endif

               if( tkrtspec ) then
                  stime =  atrack.t 
                  if(reducedtime .eq. 1) then
                     delta =  r*(cos(fai) + 1.)*sint*1.d9/c ! ns
                     stime = stime + delta
                  endif
                  ir = ridx
                  if(ir .lt. 1) then
                     call kwhist( tspec0(code, i), 
     *                  stime, wt)
                  elseif(ir .le. nrbin) then
                     call kwhist( tspec(code, ir, faiidx,  i), 
     *                    stime, wt)
                  endif
               endif
            endif
#endif
         endif
      endif
      return
!     **************
      entry xenevent
!     **************
!        replace  @ # % in basefilename by hostname, etc
!        and put it in basefilename2
!        
      
      write(0,*) 'ev#=',eventno,
     *   ' generation phase finished. now writing data'

      call cgetfname(basefilename, basefilename2)
      call cqfirstid(firstcdepth)
      firstcdepth = firstcdepth* 0.1     ! in g/cm2  First col depth.



      if(observeas) then
         cog = 0.
         sumne = 0.
         do i = 1, noofassites
            if(i .gt. 1 .and. i  .lt. noofassites ) then
               dd =(asdepthlist(i+1) - asdepthlist(i-1))/2.0
            elseif(i .eq. 1) then
               dd =(asdepthlist(2) - asdepthlist(1))
            else
               dd =(asdepthlist(noofassites) -
     *              asdepthlist(noofassites-1))
            endif
            cog = cog + asobssites(i).esize*dd*asdepthlist(i)
            sumne= sumne +asobssites(i).esize*dd
         enddo
!          0.1 is for g/cm2
         cog = cog*0.1/sumne

         cog2 = 0.
         sumne = 0.
         do i = 1, noofassites
            if( asobssites(i).age .gt.
     *          (2.0-asobssites(noofassites).age))  then
               if(i .gt. 1 .and. i  .lt. noofassites ) then
                  dd =( asdepthlist(i+1) - asdepthlist(i-1))/2.0
               elseif(i .eq. 1) then
                  dd =(asdepthlist(2) - asdepthlist(1))
               else
                  dd =(asdepthlist(noofassites) -
     *              asdepthlist(noofassites-1))
               endif
               dd = dd
               cog2 = cog2 + asobssites(i).esize*asdepthlist(i)*dd
               sumne= sumne +asobssites(i).esize*dd
            endif
         enddo
         if(sumne .gt. 0.) then
            cog2 = cog2*0.1/sumne
         else
!              to deep penetration
            cog2 = asdepthlist(noofassites)*0.1
         endif

         filename = basefilename2(1:klena(basefilename2))//".hyb"
         call copenfw2(fnob, filename, 1, icon)

         write(fnob,
     *   '("h ", i4,  3i3, 1pE11.3, 0p 3f11.7, 1pE11.3, 0p,
     *     2f7.0,i2,a )')
     *      eventno,  inci.p.code,
     *      inci.p.subcode, inci.p.charge,
     *      inci.p.fm.e, -angle.r(1), -angle.r(2), -angle.r(3),
     *      firstcdepth, cog, cog2, vn, ' /'

         sumebydedx = 0.
         sumebydeathg =0.
         sumebydeathe =0.
         sumebydeathmupik =0.
         sumebydeathneu = 0.
         sumebydeathnut = 0.
         sumebydeathp = 0.
         sumebydeatho = 0.
         sumebydeath = 0.
         sumuncertain = 0.
         sumecrash = 0.
         sumespace = 0.

         do i = 1, noofassites 
            if(eabsorb(1) .ne. 0) then
               write(fnob, '("t ", i3, 2f7.1,  2f6.3,
     *         1p14g12.3 )')
     *           i, 
     *          asdepthlist(i)*0.1,  asobssites(i).mu,
     *          asobssites(i).age,   asdepthlist(i)*0.1/cog2, 
     *          ng(i), ne(i), nmu(i), nhad(i),
     *          asobssites(i).esize, sumeloss(i), 
     *          debydedx(i), debydeath(i), 
!                   next ones are from 7.51
     *          debydeathg(i),  debydeathe(i), debydeathmupik(i), 
     *          debydeathp(i),  debydeathnut(i), debydeatho(i)

               if(i .le. eabsorb(2) ) then
!                    to see E consv. we should not count
!                    level > Eabsorb(2).
                  sumebydedx = sumebydedx + debydedx(i)
                  sumebydeath = sumebydeath + debydeath(i)
                  sumebydeathneu = sumebydeathneu +debydeathneu(i)
                  sumebydeathnut = sumebydeathnut +debydeathnut(i)
                  sumebydeathg = sumebydeathg + debydeathg(i)
                  sumebydeathe = sumebydeathe + debydeathe(i)
                  sumebydeathmupik = sumebydeathmupik +
     *                   debydeathmupik(i)
                  sumebydeathp = sumebydeathp +debydeathp(i)
                  sumebydeatho = sumebydeatho +debydeatho(i)
               endif
            else
               write(fnob, '("t ", i3, 2f7.1,  2f6.3,
     *         1p6E11.3 )')
     *           i, 
     *          asdepthlist(i)*0.1,  asobssites(i).mu,
     *          asobssites(i).age,   asdepthlist(i)*0.1/cog2, 
     *          ng(i), ne(i), nmu(i), nhad(i),
     *          asobssites(i).esize, sumeloss(i)
            endif
         enddo
         if(eabsorb(1) .ne. 0) then
            do i = 1, 7
               sumecrash = sumecrash + ecrash(i)
               sumespace = sumespace + espace(i)
            enddo
            write(fnob,'("b ", 1p7E11.3)') (espace(i), i=1,7)
            write(fnob,'("b ", 1p7E11.3, i4)') (ecrash(i), i=1,7),
     *      eabsorb(2)
            write(fnob,
     *       '("c ",1p7E11.3)' ) 
     *      maxebreak, maxrelebreak, sumediff, sumabsediff,
     *      maxebreak(1)/inci.p.fm.p(4)

            summissing =  sumecrash + sumespace + sumebydeathneu
            sumuncertain = sumebydeathnut
            sumdeinair =  sumebydedx + sumebydeath
            sumall = sumdeinair + summissing + sumuncertain

            write(fnob,'("s ", 1p8E11.3)') 
     *       sumebydedx, sumebydeath, sumdeinair,
     *       sumecrash, sumespace, sumebydeathnut,
     *       sumebydeathneu, sumall

            write(fnob,'("r ", 1p4E11.3)') 
     *      sumdeinair/e0, sumuncertain/e0, summissing/e0, sumall/e0
!                normalized one
            write(fnob,'("n ", 1p4E11.3)') 
     *      sumdeinair/sumall, sumuncertain/sumall,summissing/sumall,
     *      1.0
!                additional info for more details
            write(fnob,'("a ", 1p5g12.3 )')
     *      sumebydeathg,  sumebydeathe,  sumebydeathmupik,
     *      sumebydeathp,  sumebydeatho
         endif
         write(fnob,*)
         close(fnob)
      endif

      do i = 1, ansites
         j=histdep(i)
         write(evid(i), 
     *   '(i3, i5,  f5.2, f5.2,
     *   f7.1,  i4)')
     *   histdep(i), int( asdepthlist(j)*0.1 ),  
     *   asobssites(j).age, asdepthlist(j)*0.1/cog2,
     *   asobssites(j).mu, int(cog2)
      enddo
#if defined  (MACOSX)
#else
      if(tkarspec) then
         filename = basefilename2(1:klena(basefilename2))//"-r.hist"
         call copenfw2(fnol, filename, binw, icon)
        do i = 1, ansites
            k=histdep(i)
            do j = 1, 4
               write(dirstr,'(a,"/d",i4, "/")')
     *              ptcln(j), int( depthlist(k)*0.1 )
               call kseblk(dirstr, "|", nstr)
               do l = 1, nfai
                  call kwhistdir(rspec(j, l,  i),  dirstr)
!                  call kwhists(  rspec(j, l, i), 0. )
                  call kwhists(  rspec(j, l, i), 0. )
                  call kwhistev( rspec(j, l, i), eventno)
                  call kwhistid( rspec(j, l,  i), evid(i))
                  call kwhistp( rspec(j, l, i),  fnol)
!                        *********** deallocate ********                            
                  call kwhistd( rspec(j, l, i) )
               enddo  ! code loop                                                   
            enddo ! fai loop                                                        
         enddo ! depth loop      

         close(fnol)
      endif

      if( tkrtspec ) then
         filename = basefilename2(1:klena(basefilename2))//"-t.hist"
         call copenfw2(fnot, filename, binw, icon)

         do i = 1, ansites
            do j = 1, 4
               call kwhists( tspec0(j,i), 0.)
               call kwhistev( tspec0(j,i), eventno)
               call kwhistid( tspec0(j,i), evid(i))
               k=histdep(i)
               dirstr = " "
               write( dirstr,'(a,"/d",i4, "/")')
     *              ptcln(j), int( asdepthlist(k)*0.1 )
               call kseblk( dirstr, "|", nstr)
               call kwhistdir( tspec0(j,i),  dirstr )
               call kwhistp( tspec0(j,i),  fnot )
!                 *********** deallocate ********         
               call kwhistd( tspec0(j,i) )
            enddo
         enddo

         do i = 1, ansites
            do j = 1, 4
               do ifai= 1, nfai
                  do ir= 1, nrbin
                     call kwhists( tspec(j,ir, ifai,i), 0.)
                     call kwhistev(tspec(j,ir, ifai,i), eventno)
                     call kwhistid( tspec(j,ir, ifai,i), evid(i))
                     dirstr = " "
                     write(dirstr,'(a,"/d",i4, "/F",i2,"/")')
     *                    ptcln(j), int( depthlist(k)*0.1), ifai
                     call kseblk(dirstr, "|", nstr)
                     call kwhistdir(tspec(j,ir, ifai,i),  dirstr)
                     call kwhistp( tspec(j,ir, ifai,i),  fnot)
!                        *********** deallocate ********                            
                     call kwhistd( tspec(j, ir, ifai, i) )
                  enddo
               enddo
            enddo  ! code loop                                                      
         enddo   ! depth loop    
         close(fnot)
      endif
#endif
!          output web data
      if(tkweb) then
         filename = basefilename2(1:klena(basefilename2))//".nrfai"
         call copenfw2(fnon, filename, 1, icon)
         
         write(fnon,
     *   '(i4,1pE11.3, 0p,i3, f8.4, 1pE11.3,0p, 4i4, 1p,8g11.3)')
     *   eventno, e0, nn, cosz, firstcdepth, nrbin, nfai, ansites,
     *   noofsites, keminobs   ! this is not exist in the older version          
!                                                                      
!           dE/dx lateral                                              
         do i = 1, noofsites
            do k = 1, nfai
               write(fnon, '("dE/dx",f7.1, 3i4)' )
     *            depthlist(i)*0.1, i, i, k
               write(fnon, '(1p10E11.3)')
     *             ( derfai(m,k,i), m=1,nrbin  ), decent(i)
!                                same  center value is put for all fai
            enddo
         enddo
         close(fnon)
      endif

      write(0,*) 'ev#=',eventno,' finished completely'

      end