program Optical ********************************************************************* * Geant Monte Carlo for Optical photon tracing * 5 volumes: GLA1:GLA2:GLA3:GL2:GLA1 * All volumes are almost same and made of Si02 * X/Y sizes are 1mx1m, Z thicknes: GLA1:3mm, GLA2:10um, GLA3:1mm * and they have different initial refractive index: * GLA1: 1.5 GLA2: 1.42 and GLA3: 1.62 * All volumes are inside MARS volume, which is Air * Optical characteristics for Air: Rindex = 1. is given as well as * for all 3 SiO2 volumesL GLA1,GLA2 and GLA3 * Photons are produced in Ebergy range 1.75-7.0 eV just in frint * of SiO2 volumes and propagated into the above described volumes * * 04/30/2006 Vasily Dzhordzhadze * ********************************************************************* parameter (mxhbk =1000000) common /pawc/ h(mxhbk) parameter (kwbank=10000000) COMMON /GCBANK/ NZEBRA,GVERSN,ZVERSN,IXSTOR,IXDIV,IXCONS, + FENDQ(16),LMAIN,LR1,WS(kwbank) DIMENSION IQ(2),Q(2),LQ(8000),ISW(2) EQUIVALENCE (Q(1),IQ(1),LQ(9)), (LQ(1),LMAIN),(ISW(1),WS(1)) EQUIVALENCE (JCG, JGSTAT) COMMON/GCLINK/JDIGI ,JDRAW ,JHEAD ,JHITS ,JKINE ,JMATE ,JPART + ,JROTM ,JRUNG ,JSET ,JSTAK ,JGSTAT,JTMED ,JTRACK,JVERTX + ,JVOLUM,JXYZ ,JGPAR ,JGPAR2,JSKLT C C allocate memory for ZEBRA and HBOOK CALL TIMEST(40000000.) call gzebra (kwbank) mxahbk = -mxhbk call hlimit (mxahbk) * * initialization phase call uginit * * processing phase call grun * * termination phase call uglast * stop end C =============================================================== subroutine uginit * * user initialization * *------------------------------------------------------- * COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN + ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2) COMMON/GCFLAX/BATCH, NOLOG LOGICAL BATCH, NOLOG REAL UB(2),bratio(2) INTEGER mode(2),NWB DATA NWB/2/ * * initialize GEANT variables call ginit * * read data cards * open the FFREAD cards file CALL FFSET('LINP',9) OPEN(9,ERR=991,FILE='cardsfile', STATUS='UNKNOWN') OPEN(6,ERR=991,FILE='outputfile', STATUS='UNKNOWN') CALL GFFGO CLOSE(9) * * initialize data strutctures call gzinit * * initialize graphics if (iswit(3) .ne. 0) then call iginit(0) call igsse(6,4) call igmeta(10,-111) call igsa(0) call gdinit endif * * define standard particles and materials call gpart call gmate * * define gemoetrical setup call ugeom * * compute cross-sections and energy loss tables call gphysi * * initialize user histograms call hbhist * 991 continue end C =========================================================== subroutine hbhist * histogram booking * initialize user histograms common/cvari/xph,yph,zph,phx,phy,phz,pphot,phthet,phphi,vl common/ivari/thet,vol,phzz,ipho,istp,iev real thet,vol,phzz integer ipho,istp,iev * initialize user histograms * *-------------- call hbook1(1,'x target',50,-2.5,2.5,0.) call hbook1(2,'y target',50,-2.5,2.5,0.) call hbook1(3,'z target',70,-350.0,350.0,0.) call hbook1(4,'momentum ',100,0.,0.01,0.) call hbook1(5,'Geant particle type',15,0.,15.0,0.) call hbook1(6,'photon number ',600,0.,6000.0,0.) call hbook1(7,'photon energy distribution ', * 100,0.,0.00000001,0.) call hbook1(8,'cerenkov theta ',180,0.,180.,0.) call hbook1(9,'cerenkov phi ',360,0.,360.,0.) call hropen(101,'myfile1','photon0.ntp','N',1024,istat) call hbnt(721,'iphot',' ') call hbname(721,'par',thet,'thet:r,vol:r,phzz:r,'// * 'ipho:i,istp:i,iev:i') call hropen(102,'myfile2','photon.ntp','N',1024,istat) call hbnt(722,'phot',' ') call hbname(722,'parz',xph,'xph:r,yph:r,zph:r,phx:r,phy:r,'// * 'phz:r,pphot:r,phthet:r,phphi:r,vl:r') end ************************************************************* SUBROUTINE UGEOM C. C. DEFINE USER GEOMETRY SET UP C. C.---------------------------------------------------------------------- C COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN + ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2) COMMON/GCFLAX/BATCH, NOLOG LOGICAL BATCH, NOLOG C COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP + ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD C INTEGER MXGKIN PARAMETER (MXGKIN=100) COMMON/GCKING/KCASE,NGKINE,GKIN(5,MXGKIN), + TOFD(MXGKIN),IFLGK(MXGKIN) INTEGER KCASE,NGKINE ,IFLGK REAL GKIN,TOFD C C COMMON/GCMATE/NMAT,NAMATE(5),A,Z,DENS,RADL,ABSL INTEGER NMAT,NAMATE REAL A,Z,DENS,RADL,ABSL C COMMON/GCNUM/NMATE ,NVOLUM,NROTM,NTMED,NTMULT,NTRACK,NPART + ,NSTMAX,NVERTX,NHEAD,NBIT COMMON /GCNUMX/ NALIVE,NTMSTO C COMMON/GCONST/PI,TWOPI,PIBY2,DEGRAD,RADDEG,CLIGHT,BIG,EMASS COMMON/GCONSX/EMMU,PMASS,AVO C COMMON/GCPHYS/IPAIR,SPAIR,SLPAIR,ZINTPA,STEPPA + ,ICOMP,SCOMP,SLCOMP,ZINTCO,STEPCO + ,IPHOT,SPHOT,SLPHOT,ZINTPH,STEPPH + ,IPFIS,SPFIS,SLPFIS,ZINTPF,STEPPF + ,IDRAY,SDRAY,SLDRAY,ZINTDR,STEPDR + ,IANNI,SANNI,SLANNI,ZINTAN,STEPAN + ,IBREM,SBREM,SLBREM,ZINTBR,STEPBR + ,IHADR,SHADR,SLHADR,ZINTHA,STEPHA + ,IMUNU,SMUNU,SLMUNU,ZINTMU,STEPMU + ,IDCAY,SDCAY,SLIFE ,SUMLIF,DPHYS1 + ,ILOSS,SLOSS,SOLOSS,STLOSS,DPHYS2 + ,IMULS,SMULS,SOMULS,STMULS,DPHYS3 + ,IRAYL,SRAYL,SLRAYL,ZINTRA,STEPRA COMMON/GCPHLT/ILABS,SLABS,SLLABS,ZINTLA,STEPLA + ,ISYNC + ,ISTRA * COMMON/GCSETS/IHSET,IHDET,ISET,IDET,IDTYPE,NVNAME,NUMBV(20) C COMMON/GCTIME/TIMINT,TIMEND,ITIME,IGDATE,IGTIME INTEGER ITIME,IGDATE,IGTIME REAL TIMINT,TIMEND C COMMON/GCTMED/NUMED,NATMED(5),ISVOL,IFIELD,FIELDM,TMAXFD,STEMAX + ,DEEMAX,EPSIL,STMIN,CFIELD,PREC,IUPD,ISTPAR,NUMOLD COMMON/GCTLIT/THRIND,PMIN,DP,DNDL,JMIN,ITCKOV,IMCKOV,NPCKOV c COMMON/GCTLIT/ITCKOV,IMCKOV,NPCKOV C COMMON/GCUNIT/LIN,LOUT,NUNITS,LUNITS(5) INTEGER LIN,LOUT,NUNITS,LUNITS COMMON/GCMAIL/CHMAIL CHARACTER*132 CHMAIL C PARAMETER (MAXMEC=30) COMMON/GCTRAK/VECT(7),GETOT,GEKIN,VOUT(7),NMEC,LMEC(MAXMEC) + ,NAMEC(MAXMEC),NSTEP ,MAXNST,DESTEP,DESTEL,SAFETY,SLENG + ,STEP ,SNEXT ,SFIELD,TOFG ,GEKRAT,UPWGHT,IGNEXT,INWVOL + ,ISTOP ,IGAUTO,IEKBIN, ILOSL, IMULL,INGOTO,NLDOWN,NLEVIN + ,NLVSAV,ISTORY PARAMETER (MAXME1=30) COMMON/GCTPOL/POLAR(3), NAMEC1(MAXME1) C COMMON/GCVOLU/NLEVEL,NAMES(15),NUMBER(15), +LVOLUM(15),LINDEX(15),INFROM,NLEVMX,NLDEV(15),LINMX(15), +GTRAN(3,15),GRMAT(10,15),GONLY(15),GLX(3) C * INTEGER I C Parameters include file INTEGER Air,Vacuum,Glass1,Glass2,Glass3,Aerogel COMMON/CALOPAR/Air,Vacuum,Glass1,Glass2,Glass3,Aerogel DATA Air/1/ DATA Vacuum/2/ DATA Glass1/3/ DATA Glass2/4/ DATA Glass3/5/ DATA Aerogel/6/ CHARACTER*4 WORLD PARAMETER (WORLD='CDET') REAL PAR(3) C GLASS1 REAL AGLASS1(2),ZGLASS1(2),WGLASS1(2) REAL DENGLASS1 C GLASS2 REAL AGLASS2(2),ZGLASS2(2),WGLASS2(2) REAL DENGLASS2 C GLASS3 REAL AGLASS3(2),ZGLASS3(2),WGLASS3(2) REAL DENGLASS3 C AEROGEL REAL AAERO(3),ZAERO(3),WAERO(3),WAERO1(3) C AEROGEL DATA AAERO/16.00, 1.00, 28.08/ DATA ZAERO/ 8.00, 1.00, 14.00/ DATA WAERO/ 3.00, 2.00, 1.00/ DATA WAERO1/0.7263373, 0.06113339, 0.2125293/ DATA DENAERO/0.2/ C C ---- GLASS Si02 C DATA AGLASS1/28.09 , 15.9994/ DATA ZGLASS1/14.0 , 8.0/ DATA WGLASS1/ 1.0 , 2.0/ DATA DENGLASS1/2.20/ DATA AGLASS2/28.09 , 15.9994/ DATA ZGLASS2/14.2 , 8.2/ DATA WGLASS2/ 1.0 , 2.0/ DATA DENGLASS2/2.22/ DATA AGLASS3/28.09 , 15.9994/ DATA ZGLASS3/14.3 , 8.3/ DATA WGLASS3/ 1.0 , 2.0/ DATA DENGLASS3/2.23/ COMMON/Z_THICK/Z_M,Z_O,Z_T REAL Z_M,Z_O,Z_T C DATA Z_M/0.1/ DATA Z_M/0.01/ DATA Z_O/0.1/ DATA Z_T/0.001/ C cerenkov variables integer MAXNP parameter (MAXNP=50) C Photon energy at which index of refraction is given real E_index parameter (E_index=4.0) integer ce_npckov c logical ok real ce_ppckov(MAXNP) real ce_eff_va(MAXNP) real ce_eff_gl(MAXNP) real ce_absco_va(MAXNP) real ce_absco_gl1(MAXNP) real ce_absco_gl2(MAXNP) real ce_absco_gl3(MAXNP) real ce_rindex_va(MAXNP) real ce_rindex_gl1(MAXNP) real ce_rindex_gl2(MAXNP) real ce_rindex_gl3(MAXNP) real E_low, E_hi, E real n_vac,dn_vac,abs_vac real n_glass1,dn_glass1,abs_glass1 real n_glass2,dn_glass2,abs_glass2 real n_glass3,dn_glass3,abs_glass3 data abs_vac/1e5/ data abs_glass1/1e5/ data abs_glass2/1e5/ data abs_glass3/1e5/ data n_vac/1.0008/ data n_glass1/1.5/ data n_glass2/1.42/ data n_glass3/1.59/ data dn_vac/835e-6/ data dn_glass1/835e-6/ data dn_glass2/835e-6/ data dn_glass3/835e-6/ data E_low/1.75/ data E_hi/7.0/ C C-- C-- DEFINES USER PARTICULAR MATERIALS C CALL GSMIXT(21,'SiO1',AGLASS1,ZGLASS1,DENGLASS1,-2,WGLASS1) CALL GSMIXT(22,'SiO2',AGLASS2,ZGLASS2,DENGLASS2,-2,WGLASS2) CALL GSMIXT(23,'SiO2',AGLASS3,ZGLASS3,DENGLASS3,-2,WGLASS3) CALL GSMIXT(24,'AEROGEL',AAERO,ZAERO,DENAERO,-3,WAERO) cc CALL GSMIXT(25,'AEROGEL',AAERO,ZAERO,DENAERO,3,WAERO1) C Parameters include file C-- C-- DEFINES USER TRACKING MEDIA PARAMETERS C if (iswit(6).gt.0) then FIELDM = float(iswit(6)) IFIELD = 3 else fieldm = 0. IFIELD = 0 endif TMAXFD = 10.0 DEEMAX = 0.25 STMIN = 1. stemax = 100. EPSIL = 0.005 DEEMAX = 0.20 STMIN = 0.10 stemax = 1.0 EPSIL = 0.005 DEEMAX = 0.2 stemax = 0.0001 EPSIL = 0.0001 STMIN = 0.0001 CALL GSTMED(1,'AIR$' , 15 , 0 , IFIELD, * FIELDM,TMAXFD,stemax,DEEMAX, EPSIL, STMIN, 0 , 0) CALL GSTMED(2,'VACUUM$' , 16 , 0 , IFIELD, * FIELDM,TMAXFD,stemax,DEEMAX, EPSIL, STMIN, 0 , 0) CALL GSTMED(3,'SiO1' , 21 , 0 , IFIELD, * FIELDM,TMAXFD,stemax,DEEMAX, EPSIL, STMIN, 0 , 0) CALL GSTMED(4,'SiO2' , 22 , 0 , IFIELD, * FIELDM,TMAXFD,stemax,DEEMAX, EPSIL, STMIN, 0 , 0) CALL GSTMED(5,'SiO3' , 23 , 0 , IFIELD, * FIELDM,TMAXFD,stemax,DEEMAX, EPSIL, STMIN, 0 , 0) CALL GSTMED( 6,'AER0GEL$' , 24 , 0 , IFIELD, * FIELDM,TMAXFD,stemax,DEEMAX, EPSIL, STMIN, 0 , 0 ) C-- C-- DEFINES USER'S VOLUMES C C- SURROUNDINGS C C --- Mars Volume PAR(1) = 200.00 PAR(2) = 200.00 PAR(3) = 200.00 CALL GSVOLU('HALL','BOX ',Air,PAR,3,IVOLU) IF(IVOLU.lt.0) STOP 'ERROR WITH GSVOLU' c call gsrotm(1,90.0,0.0,90.0,90.0,0.0,0.0) C --- GLASS1 Si02 PAR(1) = 50.0 ! n = 1.5 PAR(2) = 50.0 PAR(3) = Z_O/2. CALL GSVOLU('GLA1','BOX ',Glass1,PAR,3,IVOLU) IF(IVOLU.lt.0) STOP 'ERROR WITH GSVOLU' C --- GLASS2 Si02 PAR(1) = 50.0 ! n = 1.42 PAR(2) = 50.0 PAR(3) = Z_T/2. CALL GSVOLU('GLA2','BOX ',Glass2,PAR,3,IVOLU) IF(IVOLU.lt.0) STOP 'ERROR WITH GSVOLU' C --- GLASS3 Si02 PAR(1) = 50.0 ! n = 1.59 PAR(2) = 50.0 PAR(3) = Z_M/2. CALL GSVOLU('GLA3','BOX ',Glass3,PAR,3,IVOLU) IF(IVOLU.lt.0) STOP 'ERROR WITH GSVOLU' CALL GSPOS('GLA3',1,'HALL',0.0,0.0,0.0,0,'ONLY') Z_ST = -Z_M/2. - Z_T/2. CALL GSPOS('GLA2',1,'HALL',0.0,0.0,Z_ST,0,'ONLY') Z_ST = Z_M/2. + Z_T/2. CALL GSPOS('GLA2',2,'HALL',0.0,0.0,Z_ST,0,'ONLY') Z_ST = -Z_M/2. - Z_T - Z_O/2. CALL GSPOS('GLA1',1,'HALL',0.0,0.0,Z_ST,0,'ONLY') Z_ST = Z_M/2. + Z_T + Z_O/2. CALL GSPOS('GLA1',2,'HALL',0.0,0.0,Z_ST,0,'ONLY') ce_npckov = MAXNP do i = 1, ce_npckov E = E_low + (E_hi-E_low)*float(i-1)/float(ce_npckov-1) ce_ppckov(i) = E * 1e-9 ! in GeV ce_absco_va(i) = abs_vac ce_absco_gl1(i) = abs_glass1 ce_absco_gl2(i) = abs_glass2 ce_absco_gl3(i) = abs_glass3 ce_eff_va(i) = 1.0 ce_eff_gl(i) = 1.0 ce_rindex_va(i) = n_vac ! + dn_vac*(E-E_index) ce_rindex_gl1(i) = n_glass1 ! + dn_glass1*(E-E_index) ce_rindex_gl2(i) = n_glass2 ! + dn_glass2*(E-E_index) ce_rindex_gl3(i) = n_glass3 ! + dn_glass3*(E-E_index) enddo write (6,*) 'npckov: ', ce_npckov n = min(4,ce_npckov) write (6,*) (ce_ppckov(i),i=1,n) write (6,*) (ce_absco_va(i),i=1,n) write (6,*) (ce_absco_gl1(i),i=1,n) write (6,*) (ce_absco_gl2(i),i=1,n) write (6,*) (ce_absco_gl3(i),i=1,n) write (6,*) (ce_eff_va(i),i=1,n) write (6,*) (ce_eff_gl(i),i=1,n) write (6,*) (ce_rindex_va(i),i=1,n) write (6,*) (ce_rindex_gl1(i),i=1,n) write (6,*) (ce_rindex_gl2(i),i=1,n) write (6,*) (ce_rindex_gl3(i),i=1,n) call gsckov (Air,ce_npckov,ce_ppckov, 1 ce_absco_va, ce_eff_va, ce_rindex_va) call gsckov (Glass1,ce_npckov,ce_ppckov, 1 ce_absco_gl1, ce_eff_gl, ce_rindex_gl1) call gsckov (Glass2,ce_npckov,ce_ppckov, 1 ce_absco_gl2, ce_eff_gl, ce_rindex_gl2) call gsckov (Glass3,ce_npckov,ce_ppckov, 1 ce_absco_gl3, ce_eff_gl, ce_rindex_gl3) C C-- CLOSE GEOMETRY BANKS. (OBLIGATORY SYSTEM ROUTINE) CALL GGCLOS C-- C-- DEFINE SENSITIVE SETS/DETS FOR HIT DETECTION C IF (ISWIT(4) .EQ. 1) THEN CALL GPPART (0) CALL GPTMED (0) CALL GPMATE (0) CALL GPVOLU (0) CALL GPROTM (0) CALL GPRINT('KINE',0) CALL GPSETS('*','*') ENDIF C END ******************************************************************** subroutine gustep * * steps in Detector set-up * *----------------------------------------------------------- * COMMON/GCSETS/IHSET,IHDET,ISET,IDET,IDTYPE,NVNAME,NUMBV(20) C COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN + ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2) COMMON/GCFLAX/BATCH, NOLOG LOGICAL BATCH, NOLOG COMMON/GCONST/PI,TWOPI,PIBY2,DEGRAD,RADDEG,CLIGHT,BIG,EMASS C COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP + ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD C INTEGER MXGKIN PARAMETER (MXGKIN=100) COMMON/GCKING/KCASE,NGKINE,GKIN(5,MXGKIN), + TOFD(MXGKIN),IFLGK(MXGKIN) INTEGER KCASE,NGKINE ,IFLGK REAL GKIN,TOFD C INTEGER MXPHOT PARAMETER (MXPHOT=800) COMMON/GCKIN2/NGPHOT,XPHOT(11,MXPHOT) INTEGER NGPHOT REAL XPHOT C COMMON/GCKIN3/GPOS(3,MXGKIN) REAL GPOS C PARAMETER (MAXMEC=30) COMMON/GCTRAK/VECT(7),GETOT,GEKIN,VOUT(7),NMEC,LMEC(MAXMEC) + ,NAMEC(MAXMEC),NSTEP ,MAXNST,DESTEP,DESTEL,SAFETY,SLENG + ,STEP ,SNEXT ,SFIELD,TOFG ,GEKRAT,UPWGHT,IGNEXT,INWVOL + ,ISTOP ,IGAUTO,IEKBIN, ILOSL, IMULL,INGOTO,NLDOWN,NLEVIN + ,NLVSAV,ISTORY PARAMETER (MAXME1=30) COMMON/GCTPOL/POLAR(3), NAMEC1(MAXME1) C COMMON/GCNUM/NMATE ,NVOLUM,NROTM,NTMED,NTMULT,NTRACK,NPART + ,NSTMAX,NVERTX,NHEAD,NBIT COMMON /GCNUMX/ NALIVE,NTMSTO COMMON/GCTMED/NUMED,NATMED(5),ISVOL,IFIELD,FIELDM,TMAXFD,STEMAX + ,DEEMAX,EPSIL,STMIN,CFIELD,PREC,IUPD,ISTPAR,NUMOLD C COMMON/GCVOLU/NLEVEL,NAMES(15),NUMBER(15), +LVOLUM(15),LINDEX(15),INFROM,NLEVMX,NLDEV(15),LINMX(15), +GTRAN(3,15),GRMAT(10,15),GONLY(15),GLX(3) * character*4 names C common/cvari/xph,yph,zph,phx,phy,phz,pphot,phthet,phphi,vl common/ivari/thet,vol,phzz,ipho,istp,iev real thet,vol,phzz integer ipho,istp,iev integer ilabel integer iphoton common/ipht/iphoton common/thetaphi/theta,phi,p real PI,theta,phi,p,vl data PI/3.141592654/ integer NS,NW,NM,MEC,I,IIP,IMIN DIMENSION MECNAM(20) LOGICAL ROTATE if(iswit(5).gt.0.)then call gpcxyz endif if(ipart.eq.4)ISTOP =10 IF (ITRTYP.EQ.7) then if(inwvol.eq.1.and.names(nlevel).eq.'HALL')then xph = vect(1) yph = vect(2) zph = vect(3) phx = vect(4)*vect(7) phy = vect(5)*vect(7) phz = vect(6)*vect(7) pphot = vect(7) phthet= acos(vect(6))*180./PI phphi = atan(vect(5)/vect(4))*180./PI if(vect(3).le.0.1) vl=1. if(vect(3).ge.0.1) vl=2. if(vect(2).ge.50.) vl=3. if(vect(2).le.-50.) vl=4. if(vect(1).ge.50.) vl=5. if(vect(1).le.-50.) vl = 6. call hf1(7,vect(7),1.) call hf1(8,phthet,1.) call hf1(9,phphi,1.) call hfnt(722) ISTOP = 10 endif if(inwvol.eq.0.and.vect(7).eq.pphot)then istp = istp + 1 endif if(inwvol.eq.1)then istp = 0 iphoton = iphoton + 1 ipho = iphoton iev = IEVENT thet = phthet vol = vl phzz = phz call hfnt(721) endif IF (iphoton.eq.1) then NGPHOT = 1 I = 1 XPHOT(1,I) = VECT(1) XPHOT(2,I) = VECT(2) XPHOT(3,I) = VECT(3) XPHOT(4,I) = VECT(4) XPHOT(5,I) = VECT(5) XPHOT(6,I) = VECT(6) XPHOT(7,I) = VECT(7) IMIN = 1 if(abs(VECT(5)).LT.XMIN) IMIN =2 if(abs(VECT(6)).LT.XMIN) IMIN =3 IF(IMIN.eq.1)THEN XPHOT(8,I) = 0. XPHOT(9,I) = -VECT(6)/SQRT(VECT(5)**2+VECT(6)**2) XPHOT(10,I)= VECT(5)/SQRT(VECT(5)**2+VECT(6)**2) ENDIF IF(IMIN.eq.2)THEN XPHOT(8,I) = -VECT(6)/SQRT(VECT(4)**2+VECT(6)**2) XPHOT(9,I) = 0. XPHOT(10,I)= VECT(4)/SQRT(VECT(4)**2+VECT(6)**2) ENDIF IF(IMIN.eq.3)THEN XPHOT(8,I) = VECT(5)/SQRT(VECT(4)**2+VECT(5)**2) XPHOT(9,I) = -VECT(4)/SQRT(VECT(4)**2+VECT(5)**2) XPHOT(10,I)= 0. ENDIF XPHOT(11,I)= 0. ISTOP = 10 endif if(inwvol.eq.1.and.NSTEP.gt.0)then I = 1 XPHOT(1,I) = VECT(1) XPHOT(2,I) = VECT(2) XPHOT(3,I) = VECT(3) XPHOT(4,I) = VECT(4) XPHOT(5,I) = VECT(5) XPHOT(6,I) = VECT(6) XPHOT(7,I) = VECT(7) XMIN = abs(VECT(4)) IMIN = 1 if(abs(VECT(5)).LT.XMIN) IMIN =2 if(abs(VECT(6)).LT.XMIN) IMIN =3 IF(IMIN.eq.1)THEN XPHOT(8,I) = 0. XPHOT(9,I) = -VECT(6)/SQRT(VECT(5)**2+VECT(6)**2) XPHOT(10,I)= VECT(5)/SQRT(VECT(5)**2+VECT(6)**2) ENDIF IF(IMIN.eq.2)THEN XPHOT(8,I) = -VECT(6)/SQRT(VECT(4)**2+VECT(6)**2) XPHOT(9,I) = 0. XPHOT(10,I)= VECT(4)/SQRT(VECT(4)**2+VECT(6)**2) ENDIF IF(IMIN.eq.3)THEN XPHOT(8,I) = VECT(5)/SQRT(VECT(4)**2+VECT(5)**2) XPHOT(9,I) = -VECT(4)/SQRT(VECT(4)**2+VECT(5)**2) XPHOT(10,I)= 0. ENDIF XPHOT(11,I)= 0. endif endif call gsxyz if(ngkine.gt.0) then call gsking(0) endif IF(NGPHOT.GT.0) then call gskpho(0) endif if(iswit(3).gt.0)then if (ievent .ge.iswit(3)+1.and.ievent .le.iswit(3)+3) then call gdcxyz endif endif * 998 continue * * 999 end c =================================================================== subroutine guout * * final operations per event * *----------------------------------------------------------- * COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN + ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2) COMMON/GCFLAX/BATCH, NOLOG LOGICAL BATCH, NOLOG C PARAMETER (MAXMEC=30) COMMON/GCTRAK/VECT(7),GETOT,GEKIN,VOUT(7),NMEC,LMEC(MAXMEC) + ,NAMEC(MAXMEC),NSTEP ,MAXNST,DESTEP,DESTEL,SAFETY,SLENG + ,STEP ,SNEXT ,SFIELD,TOFG ,GEKRAT,UPWGHT,IGNEXT,INWVOL + ,ISTOP ,IGAUTO,IEKBIN, ILOSL, IMULL,INGOTO,NLDOWN,NLEVIN + ,NLVSAV,ISTORY PARAMETER (MAXME1=30) COMMON/GCTPOL/POLAR(3), NAMEC1(MAXME1) C COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP + ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD C C C ---- TERMINATION IF WANTED NUMBER OF TRIGGERS IS STORED ON THE TAPE C if (iswit(3) .ne. 0) then call gdhits('*','*',0,851,0.1) CALL GDXYZ(0) call gdpart(0,11,0.35) call iclrwk(0,0) CALL GDSPEC('CDET') CALL GDHEAD(1101,'GLASS $',0.4) endif * end C ===================================================================== subroutine uglast * * final operations * *----------------------------------------------------------- * C COMMON/GCLINK/JDIGI ,JDRAW ,JHEAD ,JHITS ,JKINE ,JMATE ,JPART + ,JROTM ,JRUNG ,JSET ,JSTAK ,JGSTAT,JTMED ,JTRACK,JVERTX + ,JVOLUM,JXYZ ,JGPAR ,JGPAR2,JSKLT C COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN + ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2) COMMON/GCFLAX/BATCH, NOLOG LOGICAL BATCH, NOLOG C COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP + ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD C COMMON/GCUNIT/LIN,LOUT,NUNITS,LUNITS(5) INTEGER LIN,LOUT,NUNITS,LUNITS COMMON/GCMAIL/CHMAIL CHARACTER*132 CHMAIL C COMMON/GCTIME/TIMINT,TIMEND,ITIME,IGDATE,IGTIME INTEGER ITIME,IGDATE,IGTIME REAL TIMINT,TIMEND * write(6,*) 'end after',nevent,' events' * if (iswit(3) .ne. 0) then call igend endif ******************* * NEW ****************** C CALL HRPUT(0,'ceren.hst','N') call hrout(0, icycle,' ') call hrend('myfile1') call hrend('myfile2') ccc call histdo C write(6,*) '***************************' write(6,*) '* END OF RUN *' write(6,*) '***************************' C C. ------------------------------------------------------------------ C. WRITE (CHMAIL,1000) IEVENT CALL GMAIL(0,0) C CALL GRNDMQ(NRNDM(1),NRNDM(2),0,'G') C WRITE (CHMAIL,3000) NRNDM CALL GMAIL(0,0) C C COMPUTE ONE EVENT PROCESSING TIME C IF(IEVENT.GT.0)THEN CALL TIMEL(TIMLFT) XMEAN = (TIMINT - TIMLFT)/IEVENT WRITE(CHMAIL,4000)XMEAN CALL GMAIL(0,2) ENDIF C C Print ZEBRA statistics C C CALL MZEND C C CALL GPSTAT C 1000 FORMAT('1',9X,'**** NUMBER OF EVENTS PROCESSED =',I10) 3000 FORMAT(10X,'**** RANDOM NUMBER GENERATOR AFTER' +,' LAST COMPLETE EVENT ',2I12) 4000 FORMAT(10X,'**** TIME TO PROCESS ONE EVENT IS =',F10.4, + ' SECONDS') * STOP 1 * end c ================================================================= subroutine gukine * * fix kinematics * *----------------------------------------------------------- C COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN + ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2) COMMON/GCFLAX/BATCH, NOLOG LOGICAL BATCH, NOLOG C COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP + ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD COMMON/GCTMED/NUMED,NATMED(5),ISVOL,IFIELD,FIELDM,TMAXFD,STEMAX + ,DEEMAX,EPSIL,STMIN,CFIELD,PREC,IUPD,ISTPAR,NUMOLD c COMMON/GCTLIT/ITCKOV,IMCKOV,NPCKOV COMMON/GCTLIT/THRIND,PMIN,DP,DNDL,JMIN,ITCKOV,IMCKOV,NPCKOV c INTEGER MXPHOT PARAMETER (MXPHOT=800) COMMON/GCKIN2/NGPHOT,XPHOT(11,MXPHOT) INTEGER NGPHOT REAL XPHOT COMMON/Z_THICK/Z_M,Z_O,Z_T REAL Z_M,Z_O,Z_T common/cvari/xph,yph,zph,phx,phy,phz,pphot,phthet,phphi,vl common/ivari/thet,vol,phzz,ipho,istp,iev real thet,vol,phzz integer ipho,istp,iev common/ipht/iphoton common/thetaphi/theta,phi,p real pbeam(3),theta,phi,PI,p,pt integer ippr data ippr/0/ data PI/3.1415926/ * iphoton = 0 ccc istp = 0 vol = 0. vert(1) = (rndm(1.)-0.5)*100. vert(2) = (rndm(1.)-0.5)*100 vert(3) = 0.0 25 continue c p = .005 c ipart = 3 ipart = 50 gpart = ipart p=rndm(1.)*7.0e-9 if(p.lt.1.75e-9) goto 25 223 phi = rndm(1.)*360.0 222 theta = (rndm(1.))*180.0 c write(6,47) theta,phi 47 format('theta = ',f12.6,2x,'phi = ',f12.6) theta = theta*PI/180. phi = phi*PI/180. pbeam(1)= p*sin(theta)*cos(phi) pbeam(2)= p*sin(theta)*sin(phi) pbeam(3)= p*cos(theta) pt = sqrt(pbeam(1)*pbeam(1)+pbeam(2)*pbeam(2)) call hfill(1,vert(1),0.,1.) call hfill(2,vert(2),0.,1.) call hfill(3,vert(3),0.,1.) call hfill(4,p,0.,1.) call hfill(5,gpart,0.,1.) call gsvert(vert, 0, 0, 0, 0, nvxp) call gskine(pbeam,ipart,nvxp,0,0,ntd1) if(ippr.le.4)then write(6,*) 'magnetic field : ',IFIELD,FIELDM write(6,*) 'p, id ',pbeam(1),pbeam(2),pbeam(3),ipart endif ippr=ippr + 1 991 continue return end c ================================================================= subroutine gutrev * * fix kinematics * *----------------------------------------------------------- * COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN + ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2) COMMON/GCFLAX/BATCH, NOLOG LOGICAL BATCH, NOLOG COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP + ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD C INTEGER MXGKIN PARAMETER (MXGKIN=100) COMMON/GCKING/KCASE,NGKINE,GKIN(5,MXGKIN), + TOFD(MXGKIN),IFLGK(MXGKIN) INTEGER KCASE,NGKINE ,IFLGK REAL GKIN,TOFD C INTEGER MXPHOT PARAMETER (MXPHOT=800) COMMON/GCKIN2/NGPHOT,XPHOT(11,MXPHOT) INTEGER NGPHOT REAL XPHOT C COMMON/GCKIN3/GPOS(3,MXGKIN) REAL GPO common/ipht/iphoton real xphoton ******************************************************************* CALL GTREVE write(6,*) ' iphoton = ', iphoton xphoton = iphoton call hf1(6,xphoton,1.) end c ================================================================= subroutine gutrak * * *----------------------------------------------------------- * COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN + ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2) COMMON/GCFLAX/BATCH, NOLOG LOGICAL BATCH, NOLOG COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP + ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD C INTEGER MXGKIN PARAMETER (MXGKIN=100) COMMON/GCKING/KCASE,NGKINE,GKIN(5,MXGKIN), + TOFD(MXGKIN),IFLGK(MXGKIN) INTEGER KCASE,NGKINE ,IFLGK REAL GKIN,TOFD C INTEGER MXPHOT PARAMETER (MXPHOT=800) COMMON/GCKIN2/NGPHOT,XPHOT(11,MXPHOT) INTEGER NGPHOT REAL XPHOT C COMMON /GCKIN3/ GPOS(3,MXGKIN) REAL GPOS CALL GTRACK end * * FUNCTION GUPLSH(MED0,MED1) C. C. ****************************************************************** C. * * C. * * C. * ==>Called by : GLISUR * C. * * C. ****************************************************************** C. COMMON/GCTMED/NUMED,NATMED(5),ISVOL,IFIELD,FIELDM,TMAXFD,STEMAX + ,DEEMAX,EPSIL,STMIN,CFIELD,PREC,IUPD,ISTPAR,NUMOLD COMMON/GCTLIT/THRIND,PMIN,DP,DNDL,JMIN,ITCKOV,IMCKOV,NPCKOV C. C. ------------------------------------------------------------------ C. * * *** By default this defines perfect smoothness GUPLSH = 1.0 C END