interface.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	xbgrun

Function/Subroutine Documentation

xbgrun()

subroutine xbgrun

(

)

Definition at line 4 of file interface.f.

References b, bin, c, cgetfname(), charge, cmintime2websec(), copenfw2(), cos, cqfirstid(), cqincident(), csetmuonpol(), csortstack(), d, d0, e, false, histdep(), kelec, kmuon, kphoton, kseblk(), kwhist(), kwhistai(), kwhistc(), kwhistd(), kwhistdir(), kwhistev(), kwhisti(), kwhistid(), kwhistp(), kwhists(), kwhistso(), mass, ne, nfai, ng, nmu, nrbin, p, subcode, t, tkarspec, tkrtspec, true, wgt, x, xyz, and y.

      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'

Here is the call graph for this function: