C======================================================================= SUBROUTINE INPRM C----------------------------------------------------------------------- C IN(PUT) PR(I)M(ARY) C C TAKES INPUT PRIMARY ENERGY FROM SPECIFIED SPECTRUM C CHECKS INPUT VARIABLES FOR CONSISTENCY AND LIMITATIONS C WRITES DATA BASE FILE C THIS SUBROUTINE IS CALLED FROM MAIN C----------------------------------------------------------------------- IMPLICIT NONE *KEEP,ATMOS. COMMON /ATMOS/ AATM,BATM,CATM,DATM DOUBLE PRECISION AATM(5),BATM(5),CATM(5),DATM(5) *KEEP,BUFFS. COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH INTEGER MAXBUF,MAXLEN PARAMETER (MAXBUF=39*7) PARAMETER (MAXLEN=12) REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF), * RUNE(MAXBUF),DATAB(MAXBUF) INTEGER LH CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE) EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE) *KEEP,CONST. COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER *KEEP,DPMFLG. COMMON /DPMFLG/ NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM INTEGER NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM *KEEP,ELABCT. COMMON /ELABCT/ ELCUT DOUBLE PRECISION ELCUT(4) *KEEP,ETHMAP. COMMON /ETHMAP/ ECTMAP,ELEFT DOUBLE PRECISION ECTMAP,ELEFT *KEEP,LONGI. COMMON /LONGI/ APLONG,HLONG,PLONG,SPLONG,THSTEP,THSTPI, * NSTEP,LLONGI,FLGFIT DOUBLE PRECISION APLONG(0:1040,9),HLONG(0:1024),PLONG(0:1040,9), * SPLONG(0:1040,9),THSTEP,THSTPI INTEGER NSTEP LOGICAL LLONGI,FLGFIT *KEEP,MAGANG. COMMON /MAGANG/ ARRANG,ARRANR,COSANG,SINANG DOUBLE PRECISION ARRANG,ARRANR,COSANG,SINANG *KEEP,MAGNET. COMMON /MAGNET/ BX,BZ,BVAL,BNORMC,BNORM,COSB,SINB,BLIMIT DOUBLE PRECISION BX,BZ,BVAL,BNORMC REAL BNORM,COSB,SINB,BLIMIT *KEEP,NKGI. COMMON /NKGI/ SEL,SELLG,STH,ZEL,ZELLG,ZSL,DIST, * DISX,DISY,DISXY,DISYX,DLAX,DLAY,DLAXY,DLAYX, * OBSATI,RADNKG,RMOL,TLEV,TLEVCM,IALT DOUBLE PRECISION SEL(10),SELLG(10),STH(10),ZEL(10),ZELLG(10), * ZSL(10),DIST(10), * DISX(-10:10),DISY(-10:10), * DISXY(-10:10,2),DISYX(-10:10,2), * DLAX (-10:10,2),DLAY (-10:10,2), * DLAXY(-10:10,2),DLAYX(-10:10,2), * OBSATI(2),RADNKG,RMOL(2),TLEV(10),TLEVCM(10) INTEGER IALT(2) *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,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 *KEEP,PARPAE. DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE), * (CURPAR(4), PHI ), (CURPAR(5), H ), * (CURPAR(6), T ), (CURPAR(7), X ), * (CURPAR(8), Y ), (CURPAR(9), CHI ), * (CURPAR(10),BETA), (CURPAR(11),GCM ), * (CURPAR(12),ECM ) *KEEP,PRIMSP. COMMON /PRIMSP/ PSLOPE,LLIMIT,ULIMIT,LL,UL,SLEX,ISPEC DOUBLE PRECISION PSLOPE,LLIMIT,ULIMIT,LL,UL,SLEX INTEGER ISPEC *KEEP,RANDPA. COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR DOUBLE PRECISION FAC,U1,U2 REAL RD(3000) INTEGER ISEED(103,10),NSEQ LOGICAL KNOR *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 *KEEP,STACKF. COMMON /STACKF/ STACK,STACKP,EXST,NSHIFT,NOUREC,ICOUNT,NTO,NFROM INTEGER MAXSTK PARAMETER (MAXSTK = 12*340*2) DOUBLE PRECISION STACK(MAXSTK) INTEGER STACKP,EXST,NSHIFT,NOUREC,ICOUNT,NTO,NFROM *KEEP,VERS. COMMON /VERS/ VERNUM,MVDATE,VERDAT DOUBLE PRECISION VERNUM INTEGER MVDATE CHARACTER*18 VERDAT *KEEP,VENUS. COMMON /VENUS/ ISH0,IVERVN,MTAR99,FVENUS,FVENSG INTEGER ISH0,IVERVN,MTAR99 LOGICAL FVENUS,FVENSG *KEEP,CEREN1. COMMON /CEREN1/ CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD, * CERSIZ,LCERFI DOUBLE PRECISION CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD REAL CERSIZ LOGICAL LCERFI *KEEP,CEREN2. COMMON /CEREN2/ PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS, * DCERX,DCERY,ACERX,ACERY, * XCMAX,YCMAX,EPSX,EPSY, * DCERXI,DCERYI,FCERX,FCERY, * XSCATT,YSCATT,CERXOS,CERYOS, * NCERX,NCERY,ICERML REAL PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS, * DCERX,DCERY,ACERX,ACERY, * XCMAX,YCMAX,EPSX,EPSY, * DCERXI,DCERYI,FCERX,FCERY, * XSCATT,YSCATT,CERXOS(20),CERYOS(20) INTEGER NCERX,NCERY,ICERML *KEND. DOUBLE PRECISION EFRAC,HEIGH,H0,OOO,THICK REAL VERVEN INTEGER I,IBL,IDPM,ILONG,ISO,J,L LOGICAL LTHIN EXTERNAL HEIGH,THICK CHARACTER*1 MARK CHARACTER*9 LSTDSN C----------------------------------------------------------------------- WRITE(MONIOU,504) 504 FORMAT(//' ',10('='),' SHOWER PARAMETERS ', 50('=') ) C WRITE ENERGY SPECTRUM TO HEADER RUNH(16) = PSLOPE RUNH(17) = LLIMIT RUNH(18) = ULIMIT EVTH(58) = PSLOPE EVTH(59) = LLIMIT EVTH(60) = ULIMIT IF ( PRMPAR(1) .GE. 6000.D0 .OR. PRMPAR(1) .LE. 0.D0 ) THEN WRITE(MONIOU,*)'INCORRECT SELECTION OF PRIMARY PARTICLE TYPE = ' * ,INT(PRMPAR(1)) WRITE(MONIOU,*)'PLEASE READ THE MANUALS' STOP ENDIF C CHECK WETHER NUCLEUS IS A SINGLE NUCLEON IF (PRMPAR(1) .EQ. 100.D0 ) PRMPAR(1) = 13.D0 IF (PRMPAR(1) .EQ. 101.D0 ) PRMPAR(1) = 14.D0 WRITE(MONIOU,*)'PRIMARY PARTICLE IDENTIFICATION IS ', * NINT(PRMPAR(1)) C CHECK RECOMMENDED ENERGY RANGE IF ( FVENUS .AND. * ULIMIT.GT.2.D7 .AND. PRMPAR(1).GE.8.D0 ) THEN WRITE(MONIOU,502) ULIMIT 502 FORMAT(' INTERACTION MODEL DOUBTFUL FOR THE SELECTED PRIMARY ', * 'ENERGY OF ',E10.3,' GEV'/' PLEASE READ THE MANUALS') STOP ENDIF c> *** modified by fs (22/09/98) ******************************* c IF ( PRMPAR(1) .GT. 101.D0 ) THEN c IF ( GHEISH ) THEN cC GHEISHA CAN TREAT ONLY DEUTERONS, TRITONS, AND ALPHA PARTICLES c IF ( PRMPAR(1) .NE. 201.D0 .AND. PRMPAR(1) .NE. 301.D0 c * .AND. PRMPAR(1) .NE. 402.D0 ) THEN c IF ( LLIMIT .LT. HILOELB * INT(PRMPAR(1)/100.D0) ) THEN c WRITE(MONIOU,503) INT(PRMPAR(1)/100.D0),LLIMIT c STOP c ENDIF c ENDIF c ELSE c IF ( LLIMIT .LT. HILOELB * INT(PRMPAR(1)/100.D0) ) THEN c WRITE(MONIOU,503) INT(PRMPAR(1)/100.D0),LLIMIT c 503 FORMAT(' NUCLEUS WITH A =',I2,' AND PRIMARY ENERGY =', c * 1PE10.3,' GEV TOO LOW FOR HIGH ENERGY INTERACTION MODEL'/ c * ' AND CANNOT BE TREATED BY LOW ENERGY INTERACTION MODEL'/ c * ' PLEASE READ THE MANUALS') c STOP c ENDIF c ENDIF c ENDIF c> *** end of modification **************************************** C DEFINE ENERGY RANGE AND ENERGY SPECTRUM OF PRIMARY IF ( LLIMIT .EQ. ULIMIT ) THEN ISPEC = 0 WRITE(MONIOU,506) LLIMIT 506 FORMAT(' PRIMARY ENERGY IS FIXED AT ',1PE10.3, * ' GEV' ) ELSE ISPEC = 1 WRITE(MONIOU,505) PSLOPE,LLIMIT,ULIMIT 505 FORMAT(' PRIMARY ENERGY IS TAKEN FROM SPECTRUM VIA MONTE CARLO'/ * 5X,' SLOPE OF PRIMARY SPECTRUM = ',1P,E10.3/ * 5X,' LOWER LIMIT CUT-OFF FOR PRIMARY SPECTRUM = ',E10.3,' GEV'/ * 5X,' UPPER LIMIT CUT-OFF FOR PRIMARY SPECTRUM = ',E10.3,' GEV'/) IF ( PSLOPE .NE. -1.D0 ) THEN LL = LLIMIT ** (PSLOPE + 1.D0) UL = ULIMIT ** (PSLOPE + 1.D0) SLEX = 1.D0 / (PSLOPE + 1.D0) ELSE LL = ULIMIT / LLIMIT ENDIF ENDIF C FIRST INTERACTION TARGET FIXED ? IF ( N1STTR .EQ. 1 ) THEN WRITE(MONIOU,508) 'NITROGEN' 508 FORMAT(' TARGET OF FIRST INTERACTION IS FIXED TO ',A8) ELSEIF ( N1STTR .EQ. 2 ) THEN WRITE(MONIOU,508) 'OXYGEN ' ELSEIF ( N1STTR .EQ. 3 ) THEN WRITE(MONIOU,508) 'ARGON ' ELSE N1STTR = 0 WRITE(MONIOU,*)'TARGET OF FIRST INTERACTION IS CHOSEN RANDOMLY' ENDIF C CHECK ANGULAR SETTINGS IF ( THETPR(1) .LT. 0.D0 ) THEN WRITE(MONIOU,*)'UNALLOWED CHOICE OF THETPR = ',SNGL(THETPR(1)), * ' DEGREES' WRITE(MONIOU,*)'PLEASE READ THE MANUALS' STOP ENDIF c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> c> IF ( THETPR(2) .GT. 70.D0 ) THEN c> WRITE(MONIOU,*)'UNALLOWED CHOICE OF THETPR = ',SNGL(THETPR(2)), c> * ' DEGREES' c> WRITE(MONIOU,*)'PLEASE READ THE MANUALS' c> STOP c> ELSEIF ( THETPR(2) .GT. 45.D0 ) THEN c> WRITE(MONIOU,*)'UNALLOWED CHOICE OF THETPR = ',SNGL(THETPR(2)), c> * ' DEGREES' c> WRITE(MONIOU,*)'#########################################' c> WRITE(MONIOU,*)'# IN DOUBTFUL CASES CONTACT THE AUTHORS #' c> WRITE(MONIOU,*)'#########################################' c> STOP c> ENDIF c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C INCIDENCE ANGLE FIXED ? IF ( THETPR(1) .EQ. THETPR(2) .AND. PHIPR(1) .EQ. PHIPR(2) ) THEN FIXINC = .TRUE. WRITE(MONIOU,517) THETPR(1),PHIPR(1) 517 FORMAT(' THETA OF INCIDENCE IS FIXED TO ',F10.2,' DEGREES'/ * ' PHI OF INCIDENCE IS FIXED TO ',F10.2,' DEGREES') ELSE FIXINC = .FALSE. WRITE(MONIOU,527) THETPR,PHIPR 527 FORMAT(' THETA OF INCIDENCE CHOSEN FROM ',F10.2,'...',F10.2, * ' DEGREES'/ * ' PHI OF INCIDENCE CHOSEN FROM ',F10.2,'...',F10.2, * ' DEGREES') ENDIF EVTH(81) = THETPR(1) EVTH(82) = THETPR(2) EVTH(83) = PHIPR(1) EVTH(84) = PHIPR(2) THETPR(1) = THETPR(1)*PI/180.D0 THETPR(2) = THETPR(2)*PI/180.D0 PHIPR(1) = PHIPR(1) *PI/180.D0 PHIPR(2) = PHIPR(2) *PI/180.D0 C----------------------------------------------------------------------- C PRMPAR, OBSLEV, NOBSLV PRMPAR(2) = 0.D0 PRMPAR(6) = 0.D0 PRMPAR(7) = 0.D0 PRMPAR(8) = 0.D0 C ORDERING OF OBSERVATION LEVELS FROM TOP TO BOTTOM IF ( NOBSLV .GT. 1 ) THEN 215 CONTINUE DO 11 I = 2,NOBSLV IF ( OBSLEV(I) .GT. OBSLEV(I-1) ) THEN OOO = OBSLEV(I) OBSLEV(I) = OBSLEV(I-1) OBSLEV(I-1) = OOO GOTO 215 ENDIF 11 CONTINUE ENDIF C CHECK WETHER OBSERVATION LEVELS ARE IN ALLOWED RANGE DO 12 I = 1,NOBSLV IF ( OBSLEV(I) .GE. HEIGH(0.D0) ) THEN WRITE(MONIOU,120)I,OBSLEV(I),HEIGH(0.D0) 120 FORMAT(' UNALLOWED CHOICE OF OBSLEV '/' OBSERVATION LEVEL ', * I2,' IS AT ',F12.3,' CM, WHICH IS ABOVE ', * F12.3,' CM'/' PLEASE READ THE MANUALS') STOP ENDIF IF ( OBSLEV(I) .LE. -1.D5 ) THEN WRITE(MONIOU,121)I,OBSLEV(I) 121 FORMAT(' UNALLOWED CHOICE OF OBSLEV '/' OBSERVATION LEVEL ', * I2,' IS AT ',F12.3,' CM, WHICH IS BELOW ', * '-1.D5 CM'/' PLEASE READ THE MANUALS') STOP ENDIF THCKOB(I) = THICK(OBSLEV(I)) 12 CONTINUE C WRITE OBSERVATION LEVELS TO HEADER (IN CM) RUNH(5) = REAL(NOBSLV) EVTH(47) = REAL(NOBSLV) DO 114 I = 1,NOBSLV RUNH(5+I) = OBSLEV(I) EVTH(47+I) = OBSLEV(I) 114 CONTINUE C FIRST INTERACTION HEIGHT FIXED ? IF ( FIX1I ) THEN IF ( FIXHEI .GE. HEIGH(0.D0) ) THEN WRITE(MONIOU,122)FIXHEI,HEIGH(0.D0) 122 FORMAT(' UNALLOWED CHOICE OF FIXHEI '/' FIRST INTERACTION ', * 'IS FIXED AT ',F12.3,' CM, WHICH IS ABOVE ', * F12.3,' CM'/' PLEASE READ THE MANUALS') STOP ENDIF IF ( FIXHEI .LE. OBSLEV(NOBSLV) ) THEN WRITE(MONIOU,123)FIXHEI,OBSLEV(NOBSLV) 123 FORMAT(' UNALLOWED CHOICE OF FIXHEI '/' FIRST INTERACTION ', * 'IS FIXED AT ',F12.3,' CM, '/' WHICH IS BELOW ', * 'LOWEST OBSERVATION LEVEL AT ',F12.3,' CM' * /' PLEASE READ THE MANUALS') STOP ENDIF WRITE(MONIOU,507) FIXHEI 507 FORMAT(' HEIGHT OF FIRST INTERACTION IS FIXED TO ',1P,E10.2, * ' CM') IF ( N1STTR .GE. 1 .AND. N1STTR .LE. 3 ) THEN IF ( PRMPAR(1) .LE. 3.D0 ) THEN WRITE(MONIOU,516) INT(PRMPAR(1)) 516 FORMAT(' TARGET OF FIRST INTERACTION CANNOT BE FIXED FOR ', * 'PRIMARY TYPE ',I5/' PLEASE READ THE MANUALS') STOP ELSEIF ( N1STTR .EQ. 1 ) THEN WRITE(MONIOU,*)'TARGET OF FIRST INTERACTION IS NITROGEN' ELSEIF ( N1STTR .EQ. 2 ) THEN WRITE(MONIOU,*)'TARGET OF FIRST INTERACTION IS OXYGEN' ELSEIF ( N1STTR .EQ. 3 ) THEN WRITE(MONIOU,*)'TARGET OF FIRST INTERACTION IS ARGON' ENDIF ELSE WRITE(MONIOU,*) * 'TARGET OF FIRST INTERACTION IS CHOSEN AT RANDOM' ENDIF ELSE FIXHEI = 0.D0 WRITE(MONIOU,*) 'HEIGHT OF FIRST INTERACTION IS CHOSEN RANDOMLY' ENDIF C STARTING ALTITUDE WITHIN ATMOSPHERE? IF ( THICK0 .LT. 0.D0 ) THEN WRITE(MONIOU,130)THICK0 130 FORMAT(' UNALLOWED STARTING ALTITUDE WITH NEGATIVE MASS OVERLAY' * ,E12.3/' PLEASE READ THE MANUALS') STOP ENDIF IF ( THICK0 .GE. THCKOB(NOBSLV) ) THEN WRITE(MONIOU,131) THICK0 131 FORMAT(' UNALLOWED STARTING ALTITUDE AT ',F12.3,' G/CM**2', * ' WHICH IS BELOW LOWEST OBSERVATION LEVEL'/ * ' PLEASE READ THE MANUALS') STOP ENDIF H0 = HEIGH(THICK0) IF ( THICK0 .EQ. 0.D0 ) THEN WRITE(MONIOU,518) H0, THICK0 WRITE(MONIOU,*)' WHICH IS AT TOP OF ATMOSPHERE' ELSE WRITE(MONIOU,518) H0, THICK0 ENDIF 518 FORMAT(' STARTING ALTITUDE AT ',1P,F13.2,' CM (=', * E7.1,' G/CM**2)') WRITE(MONIOU,203) (OBSLEV(I),THCKOB(I),I=1,NOBSLV) 203 FORMAT(/' OBSERVATION LEVELS IN CM AND IN G/CM**2 ', * 1P /(5X, 2E20.8 /)) C LONGITUDINAL SHOWER DEVELOPMENT IF ( LLONGI ) THEN THSTEP = NINT(THSTEP) THSTEP = MAX(THSTEP,1.D0) THSTEP = MIN(THSTEP,1040.D0) THSTPI = 1.D0/THSTEP NSTEP = INT(THCKOB(NOBSLV)*THSTPI) IF ( NSTEP .GE. 1040 ) THEN NSTEP = 1040 THSTEP = THCKOB(NOBSLV)/NSTEP WRITE(MONIOU,*)'LONGITUDINAL SHOWER SAMPLING MODIFIED' ENDIF WRITE(MONIOU,925) NSTEP+1,THSTEP 925 FORMAT(/' LONGITUDINAL SHOWER DEVELOPMENT:'/ * ' SHOWER IS SAMPLED IN ',I4, * ' STEPS OF ',F6.1,' G/CM**2') C GET HEIGHT VALUES IN CM FOR USE IN EGS DO 478 J = 0,NSTEP HLONG(J) = HEIGH(J*THSTEP) 478 CONTINUE IF ( FLGFIT ) THEN WRITE(MONIOU,*) * ' FIT TO CHARGED PARTICLE LONG. DISTRIBUTION ENABLED' ELSE WRITE(MONIOU,*) * ' FIT TO CHARGED PARTICLE LONG. DISTRIBUTION DISABLED' ENDIF WRITE(MONIOU,*) ENDIF C----------------------------------------------------------------------- C CHECK INPUT OF ENERGY CUTS IF ( GHEISH .AND. ELCUT(1) .LT. 0.05D0 ) THEN WRITE(MONIOU,*)'ELCUT(1) SELECTED INCORRECT TO ',ELCUT(1),' GEV' WRITE(MONIOU,*)'PLEASE READ THE MANUALS ' STOP ELSEIF ( .NOT.GHEISH .AND. ELCUT(1) .LT. 0.3D0 ) THEN WRITE(MONIOU,*)'ELCUT(1) SELECTED INCORRECT TO ',ELCUT(1),' GEV' WRITE(MONIOU,*)'PLEASE READ THE MANUALS ' STOP ENDIF IF ( ELCUT(2) .LT. 0.05D0 ) THEN WRITE(MONIOU,*)'ELCUT(2) SELECTED INCORRECT TO ',ELCUT(2),' GEV' WRITE(MONIOU,*)'PLEASE READ THE MANUALS ' STOP ENDIF IF ( ELCUT(3) .LT. 0.003D0 ) THEN WRITE(MONIOU,*)'ELCUT(3) SELECTED INCORRECT TO ',ELCUT(3),' GEV' WRITE(MONIOU,*)'PLEASE READ THE MANUALS ' STOP ENDIF IF ( ELCUT(4) .LT. 0.003D0 ) THEN WRITE(MONIOU,*)'ELCUT(4) SELECTED INCORRECT TO ',ELCUT(4),' GEV' WRITE(MONIOU,*)'PLEASE READ THE MANUALS ' STOP ENDIF WRITE(MONIOU,703) ECTMAP,ELCUT 703 FORMAT (' PARTICLES WITH LORENTZ FACTOR LARGER THAN',1P,E15.4, * ' ARE PRINTED OUT'/' SHOWER PARTICLES ENERGY CUT :'/ * ' FOR HADRONS : ',E15.4,' GEV'/ * ' FOR MUONS : ',E15.4,' GEV'/ * ' FOR ELECTRONS : ',E15.4,' GEV'/ * ' FOR PHOTONS : ',E15.4,' GEV'//) DO 774 I = 1,4 RUNH(20+I) = ELCUT(I) EVTH(60+I) = ELCUT(I) 774 CONTINUE C----------------------------------------------------------------------- C PARAMETERS OF EARTH MAGNETIC FIELD OF MIDDLE EUROPE C +X DIRECTION IS NORTH, +Y DIRECTION IS EAST, +Z DIRECTION IS DOWN BVAL = SQRT( BX**2 + BZ**2 ) C BNORM HAS DIMENSIONS OF MEV/CM BNORM = BVAL * C(25) * 1.D-16 C BNORMC HAS DIMENSIONS OF GEV/CM BNORMC = BNORM * 1.D-3 SINB = BZ / BVAL COSB = BX / BVAL WRITE(MONIOU,*)'EARTH MAGNETIC FIELD STRENGTH IS ',SNGL(BVAL), * ' MICROTESLA' WRITE(MONIOU,*)' WITH INCLINATION ANGLE ', * SNGL(ASIN(SINB)*180./PI),' DEGREES' IF ( BVAL .GE. 10000.D0 ) THEN WRITE(MONIOU,*)'YOU WANT TO MAGNETIZE THE GALAXY ?' WRITE(MONIOU,*)'PLEASE READ THE MANUALS !' STOP ENDIF C LIMITING FACTOR FOR STEP SIZE OF ELECTRON IN MAGNETIC FIELD BLIMIT = 0.2 / BNORM EVTH(71) = BX EVTH(72) = BZ C ANGLE BETWEEN ARRAY X-DIRECTION AND MAGNETIC NORD C POSITIV, IF X-DIRECTION OF ARRAY POINTS TO EASTERN DIRECTION ARRANR = ARRANG * PI / 180.D0 COSANG = COS(ARRANR) SINANG = SIN(ARRANR) EVTH(93) = ARRANR IF ( ARRANG .NE. 0.D0 ) THEN WRITE(MONIOU,*) WRITE(MONIOU,*)'DETECTOR COORDINATE SYSTEM IS ROTATED AWAY ', * 'FROM NORTH BY ',SNGL(ARRANG),' DEGREES' ENDIF C----------------------------------------------------------------------- C DEFINE CERENKOV ARRAY NCERX = MAX( NCERX, 1 ) NCERY = MAX( NCERY, 1 ) ACERX = ABS(ACERX) ACERY = ABS(ACERY) DCERX = MAX( ABS(DCERX), 0.001 ) DCERY = MAX( ABS(DCERY), 0.001 ) XCMAX = (ACERX + (NCERX-1) * DCERX) * 0.5 YCMAX = (ACERY + (NCERY-1) * DCERY) * 0.5 DCERXI = 1./DCERX EPSX = ACERX * 0.5 * DCERXI DCERYI = 1./DCERY EPSY = ACERY * 0.5 * DCERYI IF ( MOD(NCERX,2) .EQ. 0 ) THEN FCERX = -0.5 ELSE FCERX = 0.0 ENDIF IF ( MOD(NCERY,2) .EQ. 0 ) THEN FCERY = -0.5 ELSE FCERY = 0.0 ENDIF WRITE(MONIOU,472) * ACERX*.01,ACERY*.01, DCERX*.01,DCERY*.01, NCERX,NCERY 472 FORMAT(/' CERENKOV ARRAY:'/ * 5X,' CERENKOV STATIONS ARE ',F6.2,' X ',F6.2,' M**2 LARGE'/ * 5X,' THE GRID SPACING IS ',F6.2,' AND ',F6.2,' M',/ * 5X,' THERE ARE ',I3,' X ',I3,' STATIONS IN X/Y DIRECTIONS'/ * 5X,' THE CERENKOV ARRAY IS CENTERED AROUND (0., 0.)'/) C CALCULATE CERENKOV YIELD FACTOR FROM WAVELENGTH BAND IF ( WAVLGL .LT. 100.D0 .OR. WAVLGU .GT. 700.D0 * .OR. WAVLGL .GE. WAVLGU ) THEN WRITE(MONIOU,*)'CERENKOV WAVELENGTH BAND FROM ',SNGL(WAVLGL), * ' TO ',SNGL(WAVLGU),' NANOMETER' WRITE(MONIOU,*)' IS OUT OF VALIDITY RANGE' WRITE(MONIOU,*)'PLEASE READ THE MANUALS' STOP ENDIF WRITE(MONIOU,*)'CERENKOV WAVELENGTH BAND FROM ',SNGL(WAVLGL), * ' TO ',SNGL(WAVLGU),' NANOMETER' C WAVELENGTH IS CONVERTED FROM NM TO CM CYIELD = (1.D0/WAVLGL - 1.D0/WAVLGU) * 2.D7 * PI / C(50) C CALCULATE FACTOR FOR ETA DENSITY NORML.(ETA AT SEA LEVEL = 0.283D-3) ETADSN = 0.283D-3 * CATM(1) / BATM(1) IF ( CERSIZ .GT. 0. ) THEN WRITE(MONIOU,*)'CERENKOV BUNCH SIZE IS SET TO=',CERSIZ ELSE WRITE(MONIOU,*)'CERENKOV BUNCH SIZE IS CALCULATED FOR EACH ', * 'SHOWER' ENDIF IF ( LCERFI ) THEN WRITE(MONIOU,*)'CERENKOV PHOTONS ARE WRITTEN TO SEPARATE FILE' ELSE WRITE(MONIOU,*)'CERENKOV PHOTONS ARE WRITTEN TO PARTICLE ', * 'OUTPUT FILE' ENDIF c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> c Next block of code has been modified, and is passed to MAIN c---------------------------------------------------------------------- cC SCATTERING OF CENTER OF CHERENKOV ARRAY RELATIVE TO SHOWER AXIS c ICERML = MIN(MAX(ICERML,1),20) c XSCATT = ABS(XSCATT) c YSCATT = ABS(YSCATT) c IF ( ICERML .GE. 1 ) THEN c WRITE(MONIOU,5225)ICERML,XSCATT,YSCATT c 5225 FORMAT(' DEFINE MULTIPLE CERENKOV ARRAYS TO USE EACH', c * ' SHOWER SEVERAL TIMES'/ ' USE EACH EVENT ',I2,' TIMES'/ c * ' THE EVENTS ARE SCATTERED QUASI RANDOMLY IN THE RANGE '/ c * ' X = +- ',F10.0,' Y = +- ',F10.0) c DO 4438 I=1,ICERML c CALL SELCOR(CERXOS(I),CERYOS(I)) c WRITE(MONIOU,4437) I,CERXOS(I),CERYOS(I) c 4437 FORMAT(' CORE OF EVENT ',I2,' AT ',2F12.2) c 4438 CONTINUE c XCMAX = XCMAX + XSCATT c YCMAX = YCMAX + YSCATT c ENDIF c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C STORE CERENKOV PARAMETERS IN EVENTHEADER EVTH(86) = NCERX EVTH(87) = NCERY EVTH(88) = DCERX EVTH(89) = DCERY EVTH(90) = ACERX EVTH(91) = ACERY IF ( LCERFI ) THEN EVTH(92) = 1. ELSE EVTH(92) = 0. ENDIF EVTH(96) = WAVLGL EVTH(97) = WAVLGU c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> c Next block of code has been passed to MAIN c---------------------------------------------------------------------- c EVTH(98) = FLOAT(ICERML) c DO 480 I=1,20 c EVTH( 98+I) = CERXOS(I) c EVTH(118+I) = CERYOS(I) c 480 CONTINUE c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C----------------------------------------------------------------------- C FLAG FOR ADDITIONAL MUON INFORMATION IF ( FMUADD ) THEN WRITE(MONIOU,*) WRITE(MONIOU,*)'ADDITIONAL INFORMATION ON MUON ORIGIN IS', * ' WRITTEN TO PARTICLE TAPE' EVTH(94) = 1. ELSE EVTH(94) = 0. ENDIF C PRINTOUT OF INFORMATIONS FOR DEBUGGING IF ( DEBUG ) WRITE(MONIOU,484) MDEBUG 484 FORMAT(/' ATTENTION ! DEBUGGING IS ACTIVE'/ * ' ====> DEBUG INFORMATION WRITTEN TO UNIT ',I3//) c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> c Next block of code is obsolete. c Now it's used "jcio" routines (C) cC----------------------------------------------------------------------- cC OPEN OUTPUT DATA SET FOR RUN c IBL = INDEX(DSN,' ') c DSN(IBL:73) = 'DAT000000' c WRITE(DSN(IBL+3:IBL+8),'(I6)') NRRUN c DO 274 L = IBL+3,IBL+8 c IF ( DSN(L:L) .EQ. ' ' ) DSN(L:L) = '0' c 274 CONTINUE cC OPEN DATASET FOR PARTICLE OUTPUT c OPEN(UNIT=PATAPE,FILE=DSN,STATUS='NEW', c * FORM='UNFORMATTED',ACCESS='SEQUENTIAL') c WRITE(MONIOU,579) DSN c 579 FORMAT(/' PARTICLE OUTPUT TO DIRECTORY : ',A79) cC WRITE RUNHEADER TO OUTPUT BUFFER c CALL TOBUF( RUNH,0 ) c cC OPEN OUTPUT DATA SET FOR CERENKOV PHOTONS c IF ( LCERFI ) THEN c DSN(IBL:73) = 'CER000000' c WRITE(DSN(IBL+3:IBL+8),'(I6)') NRRUN c DO 249 L = IBL+3,IBL+8 c IF (DSN(L:L) .EQ. ' ' ) DSN(L:L) = '0' c 249 CONTINUE c OPEN(UNIT=CETAPE,FILE=DSN,STATUS='NEW', c * FORM='UNFORMATTED',ACCESS='SEQUENTIAL') c WRITE(MONIOU,580) DSN c 580 FORMAT(' CERENKOV OUTPUT TO DIRECTORY : ',A79) c CALL TOBUFC( RUNH,0 ) c ELSE c WRITE(MONIOU,580) DSN c ENDIF cC RESET DSN c DSN(IBL:73) = ' ' c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C OPEN ON EXTERNAL STACK C BLOCKS OF 32640 BYTES = 4080 REAL*8 = 340 PARTICLES OPEN(UNIT=EXST,STATUS='SCRATCH', * FORM='UNFORMATTED',ACCESS='DIRECT',RECL=MAXSTK) C----------------------------------------------------------------------- C WRITE DATA SET FOR INFORMATION BANK IF (FDBASE ) THEN C OPEN OUTPUT DATA SET FOR RUN IBL = INDEX(DSN,' ') DSN(IBL:79) = 'DAT000000.dbase' WRITE(DSN(IBL+3:IBL+8),'(I6)') NRRUN DO 275 L = IBL+3,IBL+8 IF ( DSN(L:L) .EQ. ' ' ) DSN(L:L) = '0' 275 CONTINUE OPEN(UNIT=MDBASE,FILE=DSN,STATUS='NEW') WRITE(MONIOU,581) DSN 581 FORMAT(/' DBASE OUTPUT TO DIRECTORY : ',A79) C RESET DSN DSN(IBL+9:IBL+14) = ' ' LSTDSN(1:3) = 'LST' LSTDSN(4:9) = DSN(IBL+3:IBL+8) VERVEN=FLOAT(IVERVN)/1000. IF ( LLONGI ) THEN ILONG = 1 ELSE ILONG = 0 ENDIF IF ( EVTH(75) .NE. 0. ) THEN ISO = 0 ELSE ISO = 1 ENDIF C SET DPMFLAG (0=VENUS, 1=HDPM, 2=SIBYLL, 3=QGSJET, 4=DPMJET) IF ( EVTH( 76) .NE. 0. ) THEN IDPM = 0 ELSEIF ( EVTH(139) .NE. 0. ) THEN IDPM = 2 ELSEIF ( EVTH(141) .NE. 0. ) THEN IDPM = 3 ELSEIF ( EVTH(143) .NE. 0. ) THEN IDPM = 4 ELSE IDPM = 1 ENDIF C INCREMENT DPMFLAG FOR VARIOUS CROSS SECTIONS C BY (0=HDPM-SIG, 10=VENUSSIG, 20=SIBYLLSIG, 30=QGSSIG, 4=DPMJETSIG) IF ( EVTH(145) .NE. 0 ) THEN IDPM = IDPM + 10 ELSEIF ( EVTH(140) .NE. 0 ) THEN IDPM = IDPM + 20 ELSEIF ( EVTH(142) .NE. 0 ) THEN IDPM = IDPM + 30 ELSEIF ( EVTH(144) .NE. 0 ) THEN IDPM = IDPM + 40 ENDIF MARK = '1' LTHIN = .FALSE. EFRAC = 0.D0 WRITE(MDBASE,666)VERNUM,MARK,MVDATE,VERVEN, $INT(RUNH(3))+19000000,INT(EVTH(80)),INT(EVTH(79)),INT(EVTH(78)), $INT(EVTH(77)),INT(RUNH(2)),INT(PRMPAR(1)), $LLIMIT,ULIMIT,PSLOPE,INT(RUNH(20)), $INT(RUNH(19)),INT(EVTH(76)),INT(EVTH(75)),ISO,IDPM, $NFLAIN,NFLDIF,NFLPI0,NFLPIF, $NFLCHE,NFRAGM,ILONG,THSTEP, $BX,BZ,NOBSLV, $OBSLEV(1),OBSLEV(2),OBSLEV(3), $OBSLEV(4),OBSLEV(5),OBSLEV(6), $OBSLEV(7),OBSLEV(8),OBSLEV(9), $OBSLEV(10),ELCUT(1),ELCUT(2),ELCUT(3), $ELCUT(4),EVTH(81),EVTH(82),EVTH(83), $EVTH(84),FIXHEI,N1STTR,THICK0, $STEPFC,ARRANG,INT(EVTH(94)),NSEQ, $ISEED(1,1),ISEED(2,1),ISEED(3,1),ISEED(1,2), $ISEED(2,2),ISEED(3,2),ISEED(1,3), $ISEED(2,3),ISEED(3,3),0,DSN, $LSTDSN,' JDD300.01',' JDD300.01', $NSHOW,HOST,USER,LTHIN,EFRAC 666 FORMAT('#version#',F6.3,A1,'#versiondate#',I9,'#venusversion#', $F6.3,'#rundate#',I9,/,'#computer#',I2,'#horizont#',I2,'#neutrino#' $,I2,'#cerenkov#',I2,'#runnumber#',I7,/,'#primary#',I5, $'#e_range_l#',E15.7,'#e_range_u#',E15.7,/,'#slope#',E15.7,'#nkg#', $I2,'#egs#',I2,'#venus#',I2,'#gheisha#',I2,'#isobar#',I2,'#hdpm#', $I2,/,'#hadflag1#',I2,'#hadflag2#',I2,'#hadflag3#',I2,'#hadflag4#', $I2,'#hadflag5#',I2,'#hadflag6#',I2,/,'#longi#',I2,'#longistep#', $E15.7,'#magnetx#',E15.7,/,'#magnetz#',E15.7,'#nobslev#',I3,/, $'#obslev1#',E15.7,'#obslev2#',E15.7,'#obslev3#',E15.7,/, $'#obslev4#',E15.7,'#obslev5#',E15.7,'#obslev6#',E15.7,/, $'#obslev7#',E15.7,'#obslev8#',E15.7,'#obslev9#',E15.7,/, $'#obslev10#',E15.7,'#hcut#',E15.7,'#mcut#',E15.7,/,'#ecut#',E15.7, $'#gcut#',E15.7,'#thetal#',E15.7,/,'#thetau#',E15.7,'#phil#',E15.7, $'#phiu#',E15.7,/,'#fixhei#',E15.7,'#n1sttr#',I3,'#fixchi#',E15.7, $/,'#stepfc#',E15.7,'#arrang#',E15.7,'#muaddi#',I2,'#nseq#',I2,/, $'#seq1seed1#',I9,'#seq1seed2#',I9,'#seq1seed3#',I9,/,'#seq2seed1#' $,I9,'#seq2seed2#',I9,'#seq2seed3#',I9,/,'#seq3seed1#',I9, $'#seq3seed2#',I9,'#seq3seed3#',I9,/,'#size#',I10,'#dsn_events#', $A59,/,'#dsn_prtout# ',A9,'#tape_name#',A10,'#backup#',A10,/, $'#howmanyshowers#',I10,'#host#',A20,'#user#',A20,/ $'#thinning#',L4,'#thinninglevel#',E15.7) C RESET DSN DSN(IBL:IBL+14) = ' ' CLOSE(UNIT=MDBASE) ENDIF WRITE(MONIOU,*)'NUMBER OF SHOWERS TO GENERATE =',NSHOW WRITE(MONIOU,*) RETURN END