*********************************************************************** * SUBROUTINE UGINIT(DOBATCH,INPUTFILE) * * modified on: 08/04/2005 by R. Holzmann * modified on: 26/09/2007 by R. Holzmann * 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). * modified on 13/02/2014 by I. Koenig *********************************************************************** 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, paramfile character*1 ans integer ihad(4) / 8, 9, 13, 14 / integer ikey(10) /0,0,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 emckey call rpckey call wallkey call startkey #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, paramfile) ! 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' else CALL GSTPAR(NNMED,'STRA',1.) endif NNMED = GETMEDNUM('AIR_ACTIVE$') if(nnmed.eq.0) then write(6,100) ' uginit.F: ebox medium not found' else CALL GSTPAR(NNMED,'STRA',1.) endif 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 if (ikey(9).eq.1) iemc = 1 if (ikey(10).eq.1) istart = 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 (istart.ge.1) call startinit 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 if (iemc.ge.1) call emcinit(paramfile) 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