444 lines
14 KiB
Fortran
444 lines
14 KiB
Fortran
#include "geant321/pilot.h"
|
|
c
|
|
c-----------------------------------------------------------------------------
|
|
c
|
|
subroutine guout
|
|
c
|
|
c this routine writes the data at the end of one event to NTUPLE.
|
|
c
|
|
c T. Prokscha, 15-Sep-2000, PSI
|
|
c
|
|
c Unix Version, part of the lemsr_geant program
|
|
c
|
|
c TP, 05-Apr-2001, PSI Unix, NT, Linux version
|
|
c
|
|
c-----------------------------------------------------------------------------
|
|
c
|
|
implicit none
|
|
c
|
|
#define CERNLIB_TYPE
|
|
#include "geant321/gcsets.inc"
|
|
#include "geant321/gcflag.inc"
|
|
#include "geant321/gclist.inc"
|
|
c
|
|
c ntuple variables and other common variablea
|
|
c
|
|
#include "cwn.inc"
|
|
#include "common.inc"
|
|
c
|
|
c-------------------------------------------------------------------------
|
|
c
|
|
integer*4 i, j, NVDIM, NHMAX, NHDIM
|
|
c
|
|
parameter (NVDIM = 1)
|
|
parameter (NHMAX = 100)
|
|
parameter (NHDIM = 7)
|
|
integer*4 numvs(NVDIM)
|
|
integer*4 itra(NHMAX)
|
|
integer*4 nbv(NVDIM,NHMAX)
|
|
real*4 hits(NHDIM,NHMAX)
|
|
integer*4 nhits, maxind
|
|
c
|
|
real*4 detot, descipos, descigam, desciele ! energy loss total,
|
|
c ! and in Sc. due
|
|
c ! to e+, gammas, e-
|
|
c
|
|
real*4 descil, descit, descir, descib ! energy loss in one detec.
|
|
c
|
|
integer*4 bit
|
|
integer*4 volno ! volume number; this is to decide which
|
|
c ! scintillator was hit:
|
|
c ! volno = 1 (LEFT) ----> write 1 to NTuple variable
|
|
c ! = 2 (TOP) ----> write 2
|
|
c ! = 3 (RIGHT) ----> write 4
|
|
c ! = 4 (BOTTOM)----> write 8
|
|
c ---> sum variable entry if more than one detector
|
|
c was hit.
|
|
c ! for outer scintillators, start LEFT ---> 16
|
|
c and so on.
|
|
c
|
|
integer*4 partcode ! store in NTuple which particles had contributed
|
|
c ! partcode = 1 (Gamma) ---> write 1 to NT variable
|
|
c ! partcode = 2 (e+) ---> write 2 to NT
|
|
c ! partcode = 3 (e-) ---> write 4 to NT
|
|
c ---> sum NT variable values if two or three
|
|
c different particles were generated.
|
|
c for outer scintillators start with
|
|
c GAMMA ---> 8
|
|
c
|
|
real*4 ubuf(1)
|
|
integer*4 nubuf /1/
|
|
integer*4 nvert, part
|
|
real*4 vert(3), pvert(3) ! e+ vertex values
|
|
c
|
|
c-----------------------------------------------------------------------------
|
|
c
|
|
if ( lsets(1) .eq. 5 .or. lsets(1) .eq. 500) then
|
|
c ! if initial particle was a muon,
|
|
c ! then, look for e+ momentum which
|
|
c ! is created as 2. vertex
|
|
c ! "2" is therefore the 1. argument
|
|
c ! of gfkine
|
|
c
|
|
call gfkine(2, vert, pvert, part, nvert, ubuf, nubuf)
|
|
c
|
|
nt_pxpos = 1000. * pvert(1) ! MeV/c
|
|
nt_pypos = 1000. * pvert(2) ! MeV/c
|
|
nt_pzpos = 1000. * pvert(3) ! MeV/c
|
|
c
|
|
endif
|
|
c
|
|
c now look for detector hits to get information on energy loss,
|
|
c particle and which detector has a hit, and fill NTUPLE variuables
|
|
c
|
|
c------------------------------
|
|
c
|
|
c the MCP2 detector
|
|
c
|
|
do i = 1, NHDIM ! reset first HITS(,) array
|
|
do j = 1, NHMAX
|
|
hits(i, j) = 0.
|
|
enddo
|
|
enddo
|
|
detot = 0.
|
|
c
|
|
call gfhits('MDET','DMCP',NVDIM,NHDIM,NHMAX,0,numvs,itra,nbv,
|
|
1 hits, nhits)
|
|
c
|
|
if ( nhits .gt. 0 ) then ! ~MCP2 Nhits>0
|
|
if ( nhits .le. NHMAX ) then
|
|
maxind = nhits
|
|
else
|
|
maxind = NHMAX
|
|
endif
|
|
do j = 1, maxind
|
|
c
|
|
detot = detot + hits(6,j) ! sum the energy loss per step
|
|
c
|
|
c write(6,'(1x,7E10.3)') hits(1,j),hits(2,j),hits(3,j),hits(4,j),
|
|
c + hits(5,j),hits(6,j),hits(7,j)
|
|
enddo
|
|
c
|
|
c now fill MCP relevant data to NTUPLE variables
|
|
c
|
|
nt_tmcp = int(hits(4,1)*1.e9) ! decay time at MCP in ns
|
|
nt_demcp = 1000. * detot ! total energy loss in MCP in MeV
|
|
c
|
|
endif ! ~ endif MCP2 hits >0
|
|
c
|
|
nt_mcpnhits = nhits
|
|
c
|
|
c-------------------------------------
|
|
c
|
|
c the inner scintillators
|
|
c
|
|
do i = 1, NHDIM ! reset first HITS(,) array
|
|
do j = 1, NHMAX
|
|
hits(i, j) = 0.
|
|
enddo
|
|
enddo
|
|
c
|
|
nhits = 0
|
|
detot = 0.
|
|
descipos = 0.
|
|
descigam = 0.
|
|
desciele = 0.
|
|
volno = 0
|
|
partcode = 0
|
|
descil = 0.
|
|
descit = 0.
|
|
descir = 0.
|
|
descib = 0.
|
|
c
|
|
call gfhits('ISSC','SCIS',NVDIM,NHDIM,NHMAX,0,numvs,itra,nbv,hits,
|
|
1 nhits)
|
|
c
|
|
c write(6,*) 'Number of Hits in inner Sc.: ',nhits
|
|
c
|
|
if ( nhits .gt. 0 ) then ! ~ hit in inner Sc.
|
|
c
|
|
if ( nhits .le. NHMAX ) then
|
|
maxind = nhits
|
|
else
|
|
maxind = NHMAX
|
|
endif
|
|
do j = 1, maxind
|
|
detot = detot + hits(6,j) ! sum the energy loss per step
|
|
if ( hits(7,j) .eq. 1. ) then ! particle is gamma
|
|
descigam = descigam + hits(6,j)
|
|
endif
|
|
if ( hits(7,j) .eq. 2. ) then ! particle is positron
|
|
descipos = descipos + hits(6,j)
|
|
endif
|
|
if ( hits(7,j) .eq. 3. ) then ! particle is electron
|
|
desciele = desciele + hits(6,j)
|
|
endif
|
|
c
|
|
c now look for volume
|
|
c
|
|
bit = nbv(1,j) - 1
|
|
volno = ibset(volno, bit)
|
|
c
|
|
c sum up the energy deposited in one detector
|
|
c
|
|
if ( nbv(1, j) .eq. 1 ) then ! hit is in left detector
|
|
descil = descil + hits(6,j)
|
|
endif
|
|
if ( nbv(1, j) .eq. 2 ) then ! hit is in top detector
|
|
descit = descit + hits(6,j)
|
|
endif
|
|
if ( nbv(1, j) .eq. 3 ) then ! hit is in right detector
|
|
descir = descir + hits(6,j)
|
|
endif
|
|
if ( nbv(1, j) .eq. 4 ) then ! hit is in bottom detector
|
|
descib = descib + hits(6,j)
|
|
endif
|
|
c
|
|
c look for particles
|
|
c
|
|
bit = int(hits(7,j)) - 1
|
|
partcode = ibset(partcode, bit)
|
|
c
|
|
c write(6,'(1x,7E10.3,2I5)') hits(1,j),hits(2,j),hits(3,j),hits(4,j),
|
|
c + hits(5,j),hits(6,j),hits(7,j),nbv(1,j),volno
|
|
c
|
|
enddo
|
|
c
|
|
nt_tsci = hits(4,1)*1.e9 ! decay time at Sci. in ns
|
|
nt_xsci = hits(1,1)
|
|
nt_ysci = hits(2,1)
|
|
nt_zsci = hits(3,1)
|
|
nt_desci = detot * 1000.
|
|
nt_descipos = descipos * 1000.
|
|
nt_descigam = descigam * 1000.
|
|
nt_desciele = desciele * 1000.
|
|
nt_descil = descil * 1000.
|
|
nt_descit = descit * 1000.
|
|
nt_descir = descir * 1000.
|
|
nt_descib = descib * 1000.
|
|
nt_volno = volno
|
|
nt_partcode = partcode
|
|
c
|
|
endif ! ~ hit in inner Sc.
|
|
c
|
|
nt_scinhits = nhits
|
|
c
|
|
c end of inner sc.
|
|
c
|
|
c-------------------------------
|
|
c
|
|
c the outer scintillators
|
|
c
|
|
do i = 1, NHDIM ! reset first HITS(,) array
|
|
do j = 1, NHMAX
|
|
hits(i, j) = 0.
|
|
enddo
|
|
enddo
|
|
c
|
|
nhits = 0
|
|
detot = 0.
|
|
descipos = 0.
|
|
descigam = 0.
|
|
desciele = 0.
|
|
descil = 0.
|
|
descit = 0.
|
|
descir = 0.
|
|
descib = 0.
|
|
c
|
|
c
|
|
call gfhits('OSSC','SCOS',NVDIM,NHDIM,NHMAX,0,numvs,itra,nbv,hits,
|
|
1 nhits)
|
|
c
|
|
c write(6,*) 'Number of Hits in outer Sc.: ',nhits
|
|
c
|
|
if ( nhits .gt. 0 ) then ! ~ hit in outer Sc.
|
|
c
|
|
if ( nhits .le. NHMAX ) then
|
|
maxind = nhits
|
|
else
|
|
maxind = NHMAX
|
|
endif
|
|
do j = 1, maxind
|
|
detot = detot + hits(6,j) ! sum the energy loss per step
|
|
if ( hits(7,j) .eq. 1. ) then ! particle is gamma
|
|
descigam = descigam + hits(6,j)
|
|
endif
|
|
if ( hits(7,j) .eq. 2. ) then ! particle is positron
|
|
descipos = descipos + hits(6,j)
|
|
endif
|
|
if ( hits(7,j) .eq. 3. ) then ! particle is electron
|
|
desciele = desciele + hits(6,j)
|
|
endif
|
|
c
|
|
c now look for volume
|
|
c
|
|
bit = nbv(1,j) - 1 + 4
|
|
volno = ibset(volno, bit)
|
|
c
|
|
c sum up the energy deposited in one detector
|
|
c
|
|
if ( nbv(1, j) .eq. 1 ) then ! hit is in left detector
|
|
descil = descil + hits(6,j)
|
|
endif
|
|
if ( nbv(1, j) .eq. 2 ) then ! hit is in top detector
|
|
descit = descit + hits(6,j)
|
|
endif
|
|
if ( nbv(1, j) .eq. 3 ) then ! hit is in right detector
|
|
descir = descir + hits(6,j)
|
|
endif
|
|
if ( nbv(1, j) .eq. 4 ) then ! hit is in bottom detector
|
|
descib = descib + hits(6,j)
|
|
endif
|
|
c
|
|
c look for particles
|
|
c
|
|
bit = int(hits(7,j)) - 1 + 4
|
|
partcode = ibset(partcode, bit)
|
|
c
|
|
c write(6,'(1x,7E10.3,3I5)') hits(1,j),hits(2,j),hits(3,j),hits(4,j),
|
|
c + hits(5,j),hits(6,j),hits(7,j),nbv(1,j),volno,idtype
|
|
enddo
|
|
c
|
|
c
|
|
c nt_tsco = hits(4,1)*1.e9 ! decay time at Sci. in ns
|
|
nt_xsco = hits(1,1)
|
|
nt_ysco = hits(2,1)
|
|
nt_zsco = hits(3,1)
|
|
nt_desco = detot * 1000.
|
|
nt_descopos = descipos * 1000.
|
|
nt_descogam = descigam * 1000.
|
|
nt_descoele = desciele * 1000.
|
|
nt_descol = descil * 1000.
|
|
nt_descot = descit * 1000.
|
|
nt_descor = descir * 1000.
|
|
nt_descob = descib * 1000.
|
|
nt_volno = volno
|
|
nt_partcode = partcode
|
|
c
|
|
endif ! ~ hit in inner Sc.
|
|
c
|
|
nt_sconhits = nhits
|
|
c
|
|
c end out. scintillators
|
|
c
|
|
c---------------------------------
|
|
c
|
|
c fill NTUPLE
|
|
c
|
|
call hfnt(111)
|
|
c
|
|
c reset NTUPLE variables
|
|
c
|
|
nt_x0 = 0.
|
|
nt_y0 = 0.
|
|
nt_z0 = 0.
|
|
nt_px0 = 0.
|
|
nt_py0 = 0.
|
|
nt_pz0 = 0.
|
|
nt_p0 = 0.
|
|
nt_ipart = 0
|
|
nt_xe = 0.
|
|
nt_ye = 0.
|
|
nt_ze = 0.
|
|
nt_eposend = 0.
|
|
nt_xm = 0.
|
|
nt_ym = 0.
|
|
nt_zm = 0.
|
|
nt_emend = 0.
|
|
nt_tmcp = 0
|
|
nt_tsci = 0
|
|
nt_xsci = 0.
|
|
nt_ysci = 0.
|
|
nt_zsci = 0.
|
|
nt_desci = 0.
|
|
nt_descipos = 0.
|
|
nt_descigam = 0.
|
|
nt_desciele = 0.
|
|
c
|
|
nt_xsco = 0.
|
|
nt_ysco = 0.
|
|
nt_zsco = 0.
|
|
nt_desco = 0.
|
|
nt_descopos = 0.
|
|
nt_descogam = 0.
|
|
nt_descoele = 0.
|
|
c
|
|
nt_volno = 0
|
|
nt_partcode = 0
|
|
nt_mcpnhits = 0
|
|
nt_scinhits = 0
|
|
nt_sconhits = 0
|
|
c
|
|
nt_pxpos = 0.
|
|
nt_pypos = 0.
|
|
nt_pzpos = 0.
|
|
nt_demcp = 0.
|
|
c
|
|
nt_descil = 0.
|
|
nt_descit = 0.
|
|
nt_descir = 0.
|
|
nt_descib = 0.
|
|
c
|
|
nt_descol = 0.
|
|
nt_descot = 0.
|
|
nt_descor = 0.
|
|
nt_descob = 0.
|
|
c
|
|
c-----------------------------------
|
|
c
|
|
c some information on progress of simulation
|
|
c
|
|
if ( float(idevt)/float(nevent) .ge. 0.1 .and. percent(1) .eq. 0)
|
|
1 then
|
|
write(6,*) ' fraction of processed events: 10%'
|
|
percent(1) = 1
|
|
endif
|
|
if ( float(idevt)/float(nevent) .ge. 0.2 .and. percent(2) .eq. 0)
|
|
1 then
|
|
write(6,*) ' fraction of processed events: 20%'
|
|
percent(2) = 1
|
|
endif
|
|
if ( float(idevt)/float(nevent) .ge. 0.3 .and. percent(3) .eq. 0)
|
|
1 then
|
|
write(6,*) ' fraction of processed events: 30%'
|
|
percent(3) = 1
|
|
endif
|
|
if ( float(idevt)/float(nevent) .ge. 0.4 .and. percent(4) .eq. 0)
|
|
1 then
|
|
write(6,*) ' fraction of processed events: 40%'
|
|
percent(4) = 1
|
|
endif
|
|
if ( float(idevt)/float(nevent) .ge. 0.5 .and. percent(5) .eq. 0)
|
|
1 then
|
|
write(6,*) ' fraction of processed events: 50%'
|
|
percent(5) = 1
|
|
endif
|
|
if ( float(idevt)/float(nevent) .ge. 0.6 .and. percent(6) .eq. 0)
|
|
1 then
|
|
write(6,*) ' fraction of processed events: 60%'
|
|
percent(6) = 1
|
|
endif
|
|
if ( float(idevt)/float(nevent) .ge. 0.7 .and. percent(7) .eq. 0)
|
|
1 then
|
|
write(6,*) ' fraction of processed events: 70%'
|
|
percent(7) = 1
|
|
endif
|
|
if ( float(idevt)/float(nevent) .ge. 0.8 .and. percent(8) .eq. 0)
|
|
1 then
|
|
write(6,*) ' fraction of processed events: 80%'
|
|
percent(8) = 1
|
|
endif
|
|
if ( float(idevt)/float(nevent) .ge. 0.9 .and. percent(9) .eq. 0)
|
|
1 then
|
|
write(6,*) ' fraction of processed events: 90%'
|
|
percent(9) = 1
|
|
endif
|
|
if ( idevt .eq. nevent) then
|
|
write(6,*) ' fraction of processed events: 100%'
|
|
endif
|
|
c
|
|
c-----------------------------
|
|
c
|
|
return
|
|
end
|