SUBROUTINE EGSINI C--------------------------------------------------------------------- C E(LECTRON) G(AMMA) S(HOWER) INI(TIALIZATION) C C INITIALIZES EGS4 PACKAGE C THIS SUBROUTINE IS CALLED FROM START C C DESIGN : D. HECK IK3 FZK KARLSRUHE C--------------------------------------------------------------------- *KEEP,ATMOS. COMMON /ATMOS/ AATM,BATM,CATM,DATM DOUBLE PRECISION AATM(5),BATM(5),CATM(5),DATM(5) *KEND. COMMON/BOUNDS/ECUT(6),PCUT(6),VACDST *KEEP,ELABCT. COMMON /ELABCT/ ELCUT DOUBLE PRECISION ELCUT(4) *KEND. COMMON/GEOM/ZALTIT,BOUND(6),NEWOBS,OBSLVL(10) COMMON /MEDIA/ NMED, RLC,RLDU,RLDUI,RHO,MSGE,MGE,MSEKE,MEKE,MLEKE, *MCMFP,MRANGE,IRAYLM,HBARO(6),HBAROI(6) CHARACTER MEDIA*24 COMMON/MEDIAC/MEDIA DOUBLE PRECISION PRRMMU COMMON/MUON/PRRMMU,RMMU,RMMUT2 COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(6),RHOR(6),IRAYLR(6) *KEEP,OBSPAR. COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP, * THETPR,PHIPR,NOBSLV DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10), * THETAP,THETPR(2),PHIP,PHIPR(2) INTEGER NOBSLV *KEEP,PAM. COMMON /PAM/ PAMA,SIGNUM DOUBLE PRECISION PAMA(6000),SIGNUM(6000) *KEEP,PARPAR. COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C, * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14), * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH INTEGER ITYPE,LEVL *KEND. DOUBLE PRECISION PI0MSQ COMMON/PION/PI0MSQ,PITHR,PICMAS,PI0MAS,AMASK0,AMASKC,AMASPR,AMASNT *KEEP,REJECT. COMMON /REJECT/ AVNREJ, * ALTMIN,ANEXP,THICKA,THICKD,CUTLN,EONCUT, * FNPRIM DOUBLE PRECISION AVNREJ(10) REAL ALTMIN(10),ANEXP(10),THICKA(10),THICKD(10), * CUTLN,EONCUT LOGICAL FNPRIM *KEEP,RUNPAR. COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB, * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN, * MONIOU,MDEBUG,NUCNUC, * CETAPE, * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL, * N1STTR,MDBASE, * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR, * FIX1I,FMUADD,FNKG,FPRINT,FDBASE * ,GHEISH,GHESIG COMMON /RUNPAC/ DSN,HOST,USER DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB REAL STEPFC INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC, * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL, * N1STTR,MDBASE INTEGER CETAPE CHARACTER*79 DSN CHARACTER*20 HOST,USER LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR, * FIX1I,FMUADD,FNKG,FPRINT,FDBASE * ,GHEISH,GHESIG *KEND. COMMON/THRESH/RMT2,RMSQ,ESCD2,AP,API,AE,UP,UE,TE,THMOLL DOUBLE PRECISION PZERO,PRM,PRMT2,RMI,VC COMMON/USEFUL/PZERO,PRM,PRMT2,RMI,VC,RM,MEDIUM,MEDOLD,IBLOBE,ICALL COMMON/ACLOCK/NCLOCK,JCLOCK CHARACTER MEDARR*24 DATA MEDARR/'AIR-NTP '/ C--------------------------------------------------------------------- IF((DEBUG))WRITE(MDEBUG,* )'EGSINI:' C*** INITIALISATION BEFORE THE FIRST CALL OF EGS4 ICALL = 1 IF (( DEBUG )) THEN KMPO = MDEBUG ELSE KMPO = MONIOU END IF WRITE(KMPO,10) 10 FORMAT (' START AIR SHOWER SUBROUTINE EGS4 VERSION 64S (OCT 94)'/ * ' MU+MU- PAIRS ARE FORMED IN ANALOGY WITH E+E- PAIRS. THIS IS ', * 'CERTAINLY NOT'/' CORRECT BECAUSE OF DIFFERENT FORM FACTORS. T', * 'HEREFORE TREAT MUON PAIRS WITH CAUTION.') C*** SET PARTICLE MASSES AND PHYSICAL CONSTANTS RM = PAMA(2)*1.D3 RMT2 = 2.*RM RMSQ = RM**2 PRRMMU = PAMA(5)*1.D3 RMMU = PRRMMU RMMUT2 = 2.*RMMU PICMAS = PAMA(8)*1.D3 PI0MAS = PAMA(7)*1.D3 PI0MSQ = PI0MAS**2 AMASKC = PAMA(11)*1.D3 AMASK0 = PAMA(10)*1.D3 AMASPR = PAMA(14)*1.D3 AMASNT = PAMA(13)*1.D3 C*** INVERSE OF VELOCITY OF LIGHT VC = 1.D0/C(25) C*** PION-PRODUCTION THRESHOLD PITHR = 152.0 C*** NMED AND DUNIT DEFAULT TO 1,I.E. ONE MEDIUM AND WE WORK IN CM MEDIUM=1 DO 21 I=1,24 MEDIA(I:I)=MEDARR(I:I) 21 CONTINUE 22 CONTINUE C*** BOUNDARY 1= TOP OF ATMOSPHERE (SEE SUBROUTINE HOWFAR) BOUND(1)=11300000. BOUND(2)= 4000000. BOUND(3)= 1000000. BOUND(4)= 400000. BOUND(5)= 0. BOUND(6)= -0. MED(1)=0 MED(6)=0 C*** VACUUM IN REGIONS 1 AND 6, AIR IN REGION 2 TO 5 DO 31 IRL=2,5 MED(IRL)= 1 C *** VALUE OF BARO-HEIGHT FROM HILLAS, PRIV. COMMUN. (AUG 1988) HBARO(IRL)=CATM(6-IRL) HBAROI(IRL)=1./HBARO(IRL) RHOR(IRL) =BATM(6-IRL)*HBAROI(IRL) C *** NEEDED FOR REGION 2 TO 5 SINCE NO TRANSPORT ELSEWHERE C *** ECUT IS TOTAL ENERGY C *** TERMINATE ELECTRON HISTORIES AT ECUT (GEV TO MEV CONVERTED) ECUT(IRL)=1000.D0*ELCUT(3)+RM C *** TERMINATE PHOTON HISTORIES AT PCUT (GEV TO MEV CONVERTED) PCUT(IRL)=1000.D0*ELCUT(4) 31 CONTINUE 32 CONTINUE OPEN(UNIT=KMPI,FILE='EGSDAT2',STATUS='OLD') C*** PICK UP CROSS SECTION DATA FOR AIR-NTP FROM UNIT KMPI=12 CALL HATCH CLOSE(UNIT=KMPI) C*** INVERTED PHOTON THRESHOLD API=1./AP WRITE(KMPO,40) (AE-RM)*.001,AP*.001,ECUT(2)*.001,PCUT(2)*.001 40 FORMAT (' ELECTRONS CAN BE CREATED AND ANY ELECTRON FOLLOWED DO', * 'WN TO'/T40,F8.4,' GEV KINETIC ENERGY'/' GAMMAS CAN BE CREATED', * ' AND ANY GAMMA FOLLOWED DOWN TO'/T40,F8.4,' GEV ENERGY'/' ELE', * 'CTRON HISTORIES ARE TERMINATED AT',F10.4,' GEV'/' GAMMA HISTO', * 'RIES ARE TERMINATED AT ',F10.4,' GEV'/) cc IF((DEBUG.OR.(JCLOCK.GT.1)))WRITE(KMPO,50) cc50 FORMAT (7X,' PART|TOT.ENERGY|ANGLE Z|ANGLE X|ALTITUDE|', cc * ' TIME | POS. X | POS. Y |GENER|',/,8X,'ICLE|', cc * ' (GEV) |COSTHET| (RAD) | (CM) | (MSEC) | (CM) |', cc * ' (CM) |ATION|') C*** CALCULATE THE LAYER THICKNESS BELOW EACH DETECTOR DO 61 IDET=1,NOBSLV C *** NECESSARY BECAUSE OF DOUBLE PRECIS. OBSLVL(IDET)=OBSLEV(IDET) DO 71 JREG=2,5 IF (OBSLVL(IDET).GE.BOUND(JREG)) THEN KREG=JREG GO TO 80 END IF 71 CONTINUE 72 CONTINUE WRITE(KMPO,90) IDET,OBSLVL(IDET)*0.01 90 FORMAT (' EGS4 :', ' DETECTOR ',I2,' AT ',E10.3,' M IS OUT OF A *TMOSPHERE') STOP 80 THICKD(IDET)=EXP(-OBSLVL(IDET)*HBAROI(KREG)) THICKA(IDET)=RHOR(KREG)*HBARO(KREG)*(1.-THICKD(IDET)) C *** MIN ALTITUDE FOR REJECT IS OBSERVATION LEVEL+3*36.6 G/CM**2 ALTMIN(IDET)=-HBARO(KREG)*LOG(MAX(1.E-37,(1.-(THICKA(IDET)+109.8 * )* HBAROI(KREG)/RHOR(KREG)))) ALTMIN(IDET)=MIN(ALTMIN(IDET),BOUND(1)) 61 CONTINUE 62 CONTINUE RETURN END