*********************************************************************** * SUBROUTINE UGINIT(DOBATCH,INPUTFILE) * * modified on: 08/04/2005 by R. Holzmann * modified on: 26/09/2007 by R. Holzmann * last modified on: 30/08/2012 by A. Mangiarotti * Because the RPC gap is very thin, * the energy loss fluctuations are * calculated correctly only with * ISTRA=1 (PAI model). *********************************************************************** implicit none integer dobatch character*(*) inputfile #include "geantdef.h" #include "geant321/gcflag.inc" #include "geant321/gclist.inc" #include "geant321/gckine.inc" #include "geant321/gcbank.inc" #include "geant321/gcphys.inc" #include "geant321/gctmed.inc" #include "geant321/gccuts.inc" #include "fields.inc" #include "user.inc" #include "kinetups.inc" #include "geatarget.inc" #include "mdcdconst.inc" #include "showerdconst.inc" #ifdef GEANT_REF #include "gurefp.inc" #include "guref.inc" #include "gufil.inc" #include "refkey.inc" #include "refeq.inc" #endif common /kineopt/ iopt, ival integer iopt, ival character*200 inifile, inifilein, mapfile character*1 ans integer ihad(4) / 8, 9, 13, 14 / integer ikey(8) /0,0,0,0,0,0,0,0/ integer ne, i, j, lin, iostat, istat parameter (ne=100) real ekina(ne) real zsum real bratio(6), bratio2(6) data bratio /100.0, 0.0, 0.0, 0.0, 0.0, 0.0/ data bratio2 /89.2, 10.8, 0.0, 0.0, 0.0, 0.0/ integer dmode1(6), dmode2(6), dmode3(6), dmode4(6), dmode5(6) data dmode1 /30201, 0, 0, 0, 0, 0/ ! see GEANT manual p. 87 data dmode2 /30201, 0, 0, 0, 0, 0/ data dmode3 /30207, 302, 0, 0, 0, 0/ data dmode4 /302, 0, 0, 0, 0, 0/ data dmode5 /302, 0, 0, 0, 0, 0/ integer nroot integer iquest common /quest/ iquest(100) integer nnmed integer GETMEDNUM external GETMEDNUM evtversion = 0 call mzlogl(ixstor,-2) ! suppress ZEBRA logging call ginit ! initialize GEANT * definition of user keys call hgeakey call richkey call mdckey call tofkey call showkey call rpckey call wallkey #ifdef GEANT_REF ireft = 1 * user keys for reference matrix call refkey #endif if(iseco(1).eq.0) iseco(2) = 0 ! make sure seco flags are consistent iopt = -1 ! default is no KINE printing itckov = 1 ! activate Cerenkov photons ilabs = 1 irayl = 1 cutgam = 0.0005 ! 500 keV cutele = 0.0005 cuthad = 0.001 ! 1 MeV cutneu = 0.001 do i=1,9 ! default vertex is (0,0,0) mm tgtpar(i) = 0.0 end do *------------------------------------------------------ * Read data cards * get file name for data cards, default: geaini.dat *----------------------------------------------------- lin = 1 call ffset('LINP',lin) ! Set lin to 1 for read inifile = inputfile open(unit=1,file=inifile,status='OLD') call gffgo ! read data cards from file close(unit=1) call raninit ! initialize random number seed *------------------------------------------------------------------ * Read file names (geometry, hits, output tuples, input events) * from initialisation file *------------------------------------------------------------------ call hgeantinit(inifile, mapfile) ! read all file names and open ntuples call gzinit ! initialize ZEBRA call ugeom ! create user geometry call gdinit ! initialize GEANT drawing package call geaevopen ! open event files for gukine if (imap.eq.1 .and. mapfile(1:1).ne.' ') & call ufldin(mapfile) ! read in field map if (imap.eq.2 .and. mapfile(1:1).ne.' ') & call ufldint(mapfile) ! read TOSCA map call gpart ! store standard GEANT particles call gpions ! store HI`s properties in JPART call hspecials ! store a few more particles (omega, rho, phi, Delta) if(iswit(4).eq.1) then ! allow only Dalitz decay for pi0, eta and omega call gsdk(7,bratio,dmode1) ! pi0 call gsdk(17,bratio,dmode2) ! eta call gsdk(52,bratio2,dmode3) ! omega call gsdk(41,bratio,dmode4) ! rho0 call gsdk(55,bratio,dmode5) ! phi endif *----------------------------------------------------------------- * Because the RPC gap is very thin, to get the energy loss * fluctuations right the ISTRA flag must be set to 1 (see * GEANT manual PHYS334). *----------------------------------------------------------------- NNMED = GETMEDNUM('GAS_RPC$') if(nnmed.eq.0) then write(6,100) ' uginit.F: rpc gap medium not found' endif CALL GSTPAR(NNMED,'STRA',1.) NNMED = GETMEDNUM('AIR_ACTIVE$') if(nnmed.eq.0) then write(6,100) ' uginit.F: ebox medium not found' endif CALL GSTPAR(NNMED,'STRA',1.) call gphysi ! prepare cross section tables *----------------------------------------------------------------- * Plot cross sections for hadrons in materials *----------------------------------------------------------------- if (lplot(1).gt.0 .and. lplot(1).le.100) then do i = 1, ne ekina(i) = 2.5/ne*(i-0.5) enddo do i = 1, 20 if (lplot(i).le.0 .or. lplot(i).gt.100) return do j = 1, 4 call gplmat(lplot(i),ihad(j),'ALL ',ne,ekina,-1) enddo enddo end if *----------------------------------------------------------------- * Open ROOT or rzt files for storing detector hit data *----------------------------------------------------------------- iquest(10) = 65000 call hgeahropen(ikey,nroot,itup) iroot = nroot if (ikey(1).eq.1) irich = 1 if (ikey(2).eq.1) imdc = 1 if (ikey(3).eq.1) itof = 1 if (ikey(4).eq.1) ishow = 1 if (ikey(5).eq.1) irpc = 1 if (ikey(6).eq.1) ikin = 1 if (ikey(7).eq.1) iuser = 1 if (ikey(8).eq.1) iwall = 1 call uhbook ! create histograms and ntuples if(nhsta.gt.0) call gbhsta ! book on demand "standard" GEANT histograms if(jver(1).eq.2) then ! use segmented target from geometry db zsum = 0.0 ! compute normalized integral target thickness do i=1,nSegments ! for vertex sampling in rantarg.F zsum = zsum + size(3,i) ! for different materials one should use: ! thickness(i)*density(i)/A(i)*cross-sect(i) targIntZ(i) = zsum enddo do i=1,nSegments targIntZ(i) = targIntZ(i)/zsum c write(6,*) "position = ", position(3,i) enddo c tgtpar(1) = size(1,1) ! assume all targets of same lateral size c tgtpar(2) = size(2,1) ! (in mm) end if *----------------------------------------------------------------- * prepare single-electron acceptance matrix : * acc(p,theta) = Naccepted/Nemitted *----------------------------------------------------------------- if(imacc.eq.1) then ! book histograms for acceptance calculation call hbook2(20,'e"A# emitted$',50,0.,1000.,50,10.,110.,0.) call hbook2(21,'e"A# accepted$',50,0.,1000.,50,10.,110.,0.) end if *----------------------------------------------------------------- * initialization of hit routines *----------------------------------------------------------------- if (irich.ge.1) call richinit if (imdc.ge.1 .or. imacc.eq.1) call mdcinit if (ishow.ge.1) call showinit if (itof.ge.1) call tofinit if (irpc.ge.1) call rpcinit if (iwall.ge.1) call wallinit return 100 format(A) end c c c subroutine setitup(iflag) #include "user.inc" integer iflag itup = iflag return end c c c c subroutine hspecials ! store more particles in Geant (data from PDG98) ! and/or reset default branchings #include "branchings.inc" real*4 ubuf(1) real bratio1(6), bratio2(6), bratio3(6), bratio45(6) ! in % real bratio6(6), bratio7(6), bratio0(6) real bratioD1(6), bratioD2(6) data bratio0 /98.798, 1.198, 0., 0., 0., 0./ ! pi0 data bratio1 /39.30, 32.20, 23.00, 4.75, 0.60, 7.1e-2/ ! eta data bratio2 /88.8, 8.5, 2.63, 8.3e-2, 5.9e-2, 7.15e-3/ ! omega data bratio3 /98.9, 0.99, 7.9e-2, 3.8e-2, 4.60e-3, 4.48e-3/ ! rho0 data bratio45 /99.95, 0.045, 0., 0., 0., 0./ ! rho+/rho- data bratio6 /49.1, 34.1, 15.5, 1.26, 2.99e-2, 1.3e-2/ ! phi data bratio7 /43.8, 30.2, 20.7, 3.01, 2.11, 0.154/ ! etaprime data bratioD1 /100., 0., 0., 0., 0., 0./ ! D++/D- data bratioD2 /66.3, 33.15, 0.55, 3.7e-3, 0., 0./ ! D+/D0 integer dmode1(6), dmode2(6), dmode3(6), dmode4(6), dmode5(6) integer dmode6(6), dmode7(6), dmode0(6) integer dmodeDpp(6), dmodeDp(6), dmodeD0(6), dmodeDm(6) data dmode0 /101, 30201, 0, 0, 0, 0/ ! pi0 data dmode1 /101, 70707, 80907, 80901, 30201, 10107/ ! eta data dmode2 /70809, 701, 809, 1701, 30207, 203/ ! omega data dmode3 /809, 80901, 701, 1701, 506, 203/ ! rho0 data dmode4 /807, 801, 0, 0, 0, 0/ ! rho+ data dmode5 /907, 901, 0, 0, 0, 0/ ! rho- data dmode6 /1112, 1016, 80907, 1701, 203, 170203/ ! phi data dmode7 /80917, 4101, 70717, 5201, 101, 70707/ ! etaprime data dmodeDpp /1408, 0 ,0, 0, 0, 0/ ! D++ data dmodeDp /1407, 1308, 1401, 140203, 0, 0/ ! D+ data dmodeD0 /1307, 1409 ,1301, 130203, 0, 0/ ! D0 data dmodeDm /1309, 0 ,0, 0, 0, 0/ ! D- logical isSet if (isSet(brpi0,2)) then do i=1,2 bratio0(i) = brpi0(i) ! take branching ratios from geaini.dat end do end if call normb(bratio0,2) call gsdk(7,bratio0,dmode0) ! get pi0 right if (isSet(breta,6)) then do i=1,6 bratio1(i) = breta(i) ! take branching ratios from geaini.dat end do end if call normb(bratio1,6) call gsdk(17,bratio1,dmode1) ! get eta right call gspart(52,"omega",3,0.7819,0.,7.808e-23 ,ubuf,0) ! omega meson if (isSet(bromeg,6)) then do i=1,6 bratio2(i) = bromeg(i) ! take branching ratios from geaini.dat end do end if call normb(bratio2,6) call gsdk(52,bratio2,dmode2) call gspart(41,"rho0",3,0.7699,0.,4.37e-24 ,ubuf,0) ! rho meson call gspart(42,"rho+",4,0.7699,1.,4.37e-24 ,ubuf,0) call gspart(43,"rho-",4,0.7699,-1.,4.37e-24 ,ubuf,0) if (isSet(brrho,6)) then do i=1,6 bratio3(i) = brrho(i) ! take branching ratios from geaini.dat end do end if call normb(bratio3,6) call normb(bratio45,2) call gsdk(41,bratio3,dmode3) call gsdk(42,bratio45,dmode4) call gsdk(43,bratio45,dmode5) call gspart(55,"phi",3,1.0194,0.,1.486e-22 ,ubuf,0) ! phi if (isSet(brphi,6)) then do i=1,6 bratio6(i) = brphi(i) ! take branching ratios from geaini.dat end do end if call normb(bratio6,6) call gsdk(55,bratio6,dmode6) call gspart(53,"etaprime",3,0.9577,0.,3.275e-21 ,ubuf,0) ! etaprime if (isSet(bretap,6)) then do i=1,6 bratio7(i) = bretap(i) ! take branching ratios from geaini.dat end do end if call normb(bratio7,6) call gsdk(53,bratio7,dmode7) call gspart(34,"Delta0",3,1.232,0.,5.48e-24 ,ubuf,0) ! Delta0 call gspart(35,"Delta++",4,1.232,2.,5.48e-24 ,ubuf,0) ! Delta++ call gspart(36,"Delta+",4,1.232,1.,5.48e-24 ,ubuf,0) ! Delta+ call gspart(37,"Delta-",4,1.232,-1.,5.48e-24 ,ubuf,0) ! Delta- call normb(bratioD2,4) call gsdk(34,bratioD2,dmodeD0) call gsdk(35,bratioD1,dmodeDpp) call gsdk(36,bratioD2,dmodeDp) call gsdk(37,bratioD1,dmodeDm) return end c c c subroutine normb(bratio,n) ! normalize sum of branching ratios to 1 real*4 bratio(1), sum integer*4 n, i sum = 0. do i=1,n sum = sum + bratio(i) end do sum = sum/100. do i=1,n bratio(i) = bratio(i)/sum end do return end logical function isSet(bratio,n) real*4 bratio(1), sum integer*4 n, i isset = .false. sum = 0. do i=1,n sum = sum + abs(bratio(i)) end do if (sum.gt.0.) isset = .true. return end