| 1 |
|
|---|
| 2 | PROGRAM MAIN
|
|---|
| 3 |
|
|---|
| 4 | C-----------------------------------------------------------------------
|
|---|
| 5 | C MAIN PROGRAM
|
|---|
| 6 | C
|
|---|
| 7 | C SIMULATION OF EXTENSIVE AIR SHOWERS
|
|---|
| 8 | C PREPARES INITIALIZATIONS
|
|---|
| 9 | C GENERATES SHOWERS IN THE SHOWER LOOP
|
|---|
| 10 | C TREATES PARTICLES IN THE PARTICLE LOOP
|
|---|
| 11 | C PERFORMS PRINTING OF TABLES AT END OF SHOWER AND AT END OF RUN
|
|---|
| 12 | C-----------------------------------------------------------------------
|
|---|
| 13 |
|
|---|
| 14 | IMPLICIT DOUBLE PRECISION (A-H,O-Z)
|
|---|
| 15 |
|
|---|
| 16 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 17 | parameter (xct=1)
|
|---|
| 18 | parameter (yct=2)
|
|---|
| 19 | parameter (zct=3)
|
|---|
| 20 | parameter (ctthet=4)
|
|---|
| 21 | parameter (ctphi=5)
|
|---|
| 22 | parameter (ctdiam=6)
|
|---|
| 23 | parameter (ctfoc=7)
|
|---|
| 24 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 25 |
|
|---|
| 26 | *KEEP,BAL.
|
|---|
| 27 | COMMON /BAL/ EBAL
|
|---|
| 28 | DOUBLE PRECISION EBAL(10)
|
|---|
| 29 | *KEEP,BUFFS.
|
|---|
| 30 | COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
|
|---|
| 31 | INTEGER MAXBUF,MAXLEN
|
|---|
| 32 | PARAMETER (MAXBUF=39*7)
|
|---|
| 33 | PARAMETER (MAXLEN=12)
|
|---|
| 34 | REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
|
|---|
| 35 | * RUNE(MAXBUF),DATAB(MAXBUF)
|
|---|
| 36 | INTEGER LH
|
|---|
| 37 | CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
|
|---|
| 38 | EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
|
|---|
| 39 | EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
|
|---|
| 40 | *KEEP,CHISTA.
|
|---|
| 41 | COMMON /CHISTA/ IHYCHI,IKACHI,IMUCHI,INNCHI,INUCHI,IPICHI
|
|---|
| 42 | INTEGER IHYCHI(124),IKACHI(124),IMUCHI(124),
|
|---|
| 43 | * INNCHI(124),INUCHI(124),IPICHI(124)
|
|---|
| 44 | *KEEP,CONST.
|
|---|
| 45 | COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
|
|---|
| 46 | DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
|
|---|
| 47 | *KEEP,CURVE.
|
|---|
| 48 | COMMON /CURVE/ CHAPAR,DEP,ERR,NSTP
|
|---|
| 49 | DOUBLE PRECISION CHAPAR(1100),DEP(1100),ERR(1100)
|
|---|
| 50 | INTEGER NSTP
|
|---|
| 51 | *KEEP,ELADPM.
|
|---|
| 52 | COMMON /ELADPM/ ELMEAN,ELMEAA,IELDPM,IELDPA
|
|---|
| 53 | DOUBLE PRECISION ELMEAN(37),ELMEAA(37)
|
|---|
| 54 | INTEGER IELDPM(37,13),IELDPA(37,13)
|
|---|
| 55 | *KEEP,ELASTY.
|
|---|
| 56 | COMMON /ELASTY/ ELAST,IELIS,IELHM,IELNU,IELPI
|
|---|
| 57 | DOUBLE PRECISION ELAST
|
|---|
| 58 | INTEGER IELIS(20),IELHM(20),IELNU(20),IELPI(20)
|
|---|
| 59 | *KEEP,GENER.
|
|---|
| 60 | COMMON /GENER/ GEN,ALEVEL
|
|---|
| 61 | DOUBLE PRECISION GEN,ALEVEL
|
|---|
| 62 | *KEEP,IRET.
|
|---|
| 63 | COMMON /IRET/ IRET1,IRET2
|
|---|
| 64 | INTEGER IRET1,IRET2
|
|---|
| 65 | *KEEP,ISTA.
|
|---|
| 66 | COMMON /ISTA/ IFINET,IFINNU,IFINKA,IFINPI,IFINHY
|
|---|
| 67 | INTEGER IFINET,IFINNU,IFINKA,IFINPI,IFINHY
|
|---|
| 68 | *KEEP,LONGI.
|
|---|
| 69 | COMMON /LONGI/ APLONG,HLONG,PLONG,SPLONG,THSTEP,THSTPI,
|
|---|
| 70 | * NSTEP,LLONGI,FLGFIT
|
|---|
| 71 | DOUBLE PRECISION APLONG(0:1040,9),HLONG(0:1024),PLONG(0:1040,9),
|
|---|
| 72 | * SPLONG(0:1040,9),THSTEP,THSTPI
|
|---|
| 73 | INTEGER NSTEP
|
|---|
| 74 | LOGICAL LLONGI,FLGFIT
|
|---|
| 75 | *KEEP,MPARTI.
|
|---|
| 76 | COMMON /MPARTI/ MPARTO
|
|---|
| 77 | DOUBLE PRECISION MPARTO(10,25),MPHOTO(10),MPOSIT(10),MELECT(10),
|
|---|
| 78 | * MNU(10),MMUP(10),MMUM(10),MPI0(10),MPIP(10),
|
|---|
| 79 | * MPIM(10),MK0L(10),MKPL(10),MKMI(10),MNEUTR(10),
|
|---|
| 80 | * MPROTO(10),MPROTB(10),MK0S(10),MHYP(10),
|
|---|
| 81 | * MNEUTB(10),MDEUT(10),MTRIT(10),MALPHA(10),
|
|---|
| 82 | * MOTHER(10)
|
|---|
| 83 | EQUIVALENCE (MPARTO(1, 1),MPHOTO(1)), (MPARTO(1, 2),MPOSIT(1)),
|
|---|
| 84 | * (MPARTO(1, 3),MELECT(1)), (MPARTO(1, 4),MNU(1)) ,
|
|---|
| 85 | * (MPARTO(1, 5),MMUP(1)) , (MPARTO(1, 6),MMUM(1)) ,
|
|---|
| 86 | * (MPARTO(1, 7),MPI0(1)) , (MPARTO(1, 8),MPIP(1)) ,
|
|---|
| 87 | * (MPARTO(1, 9),MPIM(1)) , (MPARTO(1,10),MK0L(1)) ,
|
|---|
| 88 | * (MPARTO(1,11),MKPL(1)) , (MPARTO(1,12),MKMI(1)) ,
|
|---|
| 89 | * (MPARTO(1,13),MNEUTR(1)), (MPARTO(1,14),MPROTO(1)),
|
|---|
| 90 | * (MPARTO(1,15),MPROTB(1)), (MPARTO(1,16),MK0S(1)) ,
|
|---|
| 91 | * (MPARTO(1,18),MHYP(1)) , (MPARTO(1,19),MDEUT(1)) ,
|
|---|
| 92 | * (MPARTO(1,20),MTRIT(1)) , (MPARTO(1,21),MALPHA(1)),
|
|---|
| 93 | * (MPARTO(1,22),MOTHER(1)), (MPARTO(1,25),MNEUTB(1))
|
|---|
| 94 | *KEEP,MULT.
|
|---|
| 95 | COMMON /MULT/ EKINL,MSMM,MULTMA,MULTOT
|
|---|
| 96 | DOUBLE PRECISION EKINL
|
|---|
| 97 | INTEGER MSMM,MULTMA(37,13),MULTOT(37,13)
|
|---|
| 98 | *KEEP,MUPART.
|
|---|
| 99 | COMMON /MUPART/ AMUPAR,BCUT,CMUON,FMUBRM,FMUORG
|
|---|
| 100 | DOUBLE PRECISION AMUPAR(14),BCUT,CMUON(11)
|
|---|
| 101 | LOGICAL FMUBRM,FMUORG
|
|---|
| 102 | *KEEP,NCOUNT.
|
|---|
| 103 | COMMON /NCOUNT/ NCOUN
|
|---|
| 104 | INTEGER NCOUN(8)
|
|---|
| 105 | *KEEP,NKGI.
|
|---|
| 106 | COMMON /NKGI/ SEL,SELLG,STH,ZEL,ZELLG,ZSL,DIST,
|
|---|
| 107 | * DISX,DISY,DISXY,DISYX,DLAX,DLAY,DLAXY,DLAYX,
|
|---|
| 108 | * OBSATI,RADNKG,RMOL,TLEV,TLEVCM,IALT
|
|---|
| 109 | DOUBLE PRECISION SEL(10),SELLG(10),STH(10),ZEL(10),ZELLG(10),
|
|---|
| 110 | * ZSL(10),DIST(10),
|
|---|
| 111 | * DISX(-10:10),DISY(-10:10),
|
|---|
| 112 | * DISXY(-10:10,2),DISYX(-10:10,2),
|
|---|
| 113 | * DLAX (-10:10,2),DLAY (-10:10,2),
|
|---|
| 114 | * DLAXY(-10:10,2),DLAYX(-10:10,2),
|
|---|
| 115 | * OBSATI(2),RADNKG,RMOL(2),TLEV(10),TLEVCM(10)
|
|---|
| 116 | INTEGER IALT(2)
|
|---|
| 117 | *KEEP,NKGS.
|
|---|
| 118 | COMMON /NKGS/ CZX,CZY,CZXY,CZYX,SAH,SL,ZNE
|
|---|
| 119 | DOUBLE PRECISION CZX(-10:10,2),CZY(-10:10,2),CZXY(-10:10,2),
|
|---|
| 120 | * CZYX(-10:10,2),SAH(10),SL(10),ZNE(10)
|
|---|
| 121 | *KEEP,NPARTI.
|
|---|
| 122 | COMMON /NPARTI/ NPARTO,MUOND
|
|---|
| 123 | DOUBLE PRECISION NPARTO(10,25),NPHOTO(10),NPOSIT(10),NELECT(10),
|
|---|
| 124 | * NNU(10),NMUP(10),NMUM(10),NPI0(10),NPIP(10),
|
|---|
| 125 | * NPIM(10),NK0L(10),NKPL(10),NKMI(10),NNEUTR(10),
|
|---|
| 126 | * NPROTO(10),NPROTB(10),NK0S(10),NHYP(10),
|
|---|
| 127 | * NNEUTB(10),NDEUT(10),NTRIT(10),NALPHA(10),
|
|---|
| 128 | * NOTHER(10),MUOND
|
|---|
| 129 | EQUIVALENCE (NPARTO(1, 1),NPHOTO(1)), (NPARTO(1, 2),NPOSIT(1)),
|
|---|
| 130 | * (NPARTO(1, 3),NELECT(1)), (NPARTO(1, 4),NNU(1)) ,
|
|---|
| 131 | * (NPARTO(1, 5),NMUP(1)) , (NPARTO(1, 6),NMUM(1)) ,
|
|---|
| 132 | * (NPARTO(1, 7),NPI0(1)) , (NPARTO(1, 8),NPIP(1)) ,
|
|---|
| 133 | * (NPARTO(1, 9),NPIM(1)) , (NPARTO(1,10),NK0L(1)) ,
|
|---|
| 134 | * (NPARTO(1,11),NKPL(1)) , (NPARTO(1,12),NKMI(1)) ,
|
|---|
| 135 | * (NPARTO(1,13),NNEUTR(1)), (NPARTO(1,14),NPROTO(1)),
|
|---|
| 136 | * (NPARTO(1,15),NPROTB(1)), (NPARTO(1,16),NK0S(1)) ,
|
|---|
| 137 | * (NPARTO(1,18),NHYP(1)) , (NPARTO(1,19),NDEUT(1)) ,
|
|---|
| 138 | * (NPARTO(1,20),NTRIT(1)) , (NPARTO(1,21),NALPHA(1)),
|
|---|
| 139 | * (NPARTO(1,22),NOTHER(1)), (NPARTO(1,25),NNEUTB(1))
|
|---|
| 140 | *KEEP,OBSPAR.
|
|---|
| 141 | COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
|
|---|
| 142 | * THETPR,PHIPR,NOBSLV
|
|---|
| 143 | DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
|
|---|
| 144 | * THETAP,THETPR(2),PHIP,PHIPR(2)
|
|---|
| 145 | INTEGER NOBSLV
|
|---|
| 146 | *KEEP,PAM.
|
|---|
| 147 | COMMON /PAM/ PAMA,SIGNUM
|
|---|
| 148 | DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
|
|---|
| 149 | *KEEP,PARPAR.
|
|---|
| 150 | COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
|
|---|
| 151 | * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
|
|---|
| 152 | DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
|
|---|
| 153 | * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
|
|---|
| 154 | INTEGER ITYPE,LEVL
|
|---|
| 155 | *KEEP,PARPAE.
|
|---|
| 156 | DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
|
|---|
| 157 | EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
|
|---|
| 158 | * (CURPAR(4), PHI ), (CURPAR(5), H ),
|
|---|
| 159 | * (CURPAR(6), T ), (CURPAR(7), X ),
|
|---|
| 160 | * (CURPAR(8), Y ), (CURPAR(9), CHI ),
|
|---|
| 161 | * (CURPAR(10),BETA), (CURPAR(11),GCM ),
|
|---|
| 162 | * (CURPAR(12),ECM )
|
|---|
| 163 | *KEEP,PBALA.
|
|---|
| 164 | COMMON /PBALA/ PBAL
|
|---|
| 165 | DOUBLE PRECISION PBAL(10)
|
|---|
| 166 | *KEEP,PRIMSP.
|
|---|
| 167 | COMMON /PRIMSP/ PSLOPE,LLIMIT,ULIMIT,LL,UL,SLEX,ISPEC
|
|---|
| 168 | DOUBLE PRECISION PSLOPE,LLIMIT,ULIMIT,LL,UL,SLEX
|
|---|
| 169 | INTEGER ISPEC
|
|---|
| 170 | *KEEP,RANDPA.
|
|---|
| 171 | COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
|
|---|
| 172 | DOUBLE PRECISION FAC,U1,U2
|
|---|
| 173 | REAL RD(3000)
|
|---|
| 174 | INTEGER ISEED(103,10),NSEQ
|
|---|
| 175 | LOGICAL KNOR
|
|---|
| 176 | *KEEP,RECORD.
|
|---|
| 177 | COMMON /RECORD/ IRECOR
|
|---|
| 178 | INTEGER IRECOR
|
|---|
| 179 | *KEEP,REJECT.
|
|---|
| 180 | COMMON /REJECT/ AVNREJ,
|
|---|
| 181 | * ALTMIN,ANEXP,THICKA,THICKD,CUTLN,EONCUT,
|
|---|
| 182 | * FNPRIM
|
|---|
| 183 | DOUBLE PRECISION AVNREJ(10)
|
|---|
| 184 | REAL ALTMIN(10),ANEXP(10),THICKA(10),THICKD(10),
|
|---|
| 185 | * CUTLN,EONCUT
|
|---|
| 186 | LOGICAL FNPRIM
|
|---|
| 187 | *KEEP,RESON.
|
|---|
| 188 | COMMON /RESON/ RDRES,RESRAN,IRESPAR
|
|---|
| 189 | REAL RDRES(2),RESRAN(1000)
|
|---|
| 190 | INTEGER IRESPAR
|
|---|
| 191 |
|
|---|
| 192 | *KEEP,RUNPAR.
|
|---|
| 193 | COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
|
|---|
| 194 | * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
|
|---|
| 195 | * MONIOU,MDEBUG,NUCNUC,
|
|---|
| 196 | * CETAPE,
|
|---|
| 197 | * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
|
|---|
| 198 | * N1STTR,MDBASE,
|
|---|
| 199 | * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
|
|---|
| 200 | * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
|
|---|
| 201 | * ,GHEISH,GHESIG
|
|---|
| 202 | COMMON /RUNPAC/ DSN,HOST,USER
|
|---|
| 203 | DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
|
|---|
| 204 | REAL STEPFC
|
|---|
| 205 | INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
|
|---|
| 206 | * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
|
|---|
| 207 | * N1STTR,MDBASE
|
|---|
| 208 | INTEGER CETAPE
|
|---|
| 209 | CHARACTER*79 DSN
|
|---|
| 210 | CHARACTER*20 HOST,USER
|
|---|
| 211 |
|
|---|
| 212 | LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
|
|---|
| 213 | * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
|
|---|
| 214 | * ,GHEISH,GHESIG
|
|---|
| 215 | *KEEP,STACKF.
|
|---|
| 216 | COMMON /STACKF/ STACK,STACKP,EXST,NSHIFT,NOUREC,ICOUNT,NTO,NFROM
|
|---|
| 217 | INTEGER MAXSTK
|
|---|
| 218 | PARAMETER (MAXSTK = 12*340*2)
|
|---|
| 219 | DOUBLE PRECISION STACK(MAXSTK)
|
|---|
| 220 | INTEGER STACKP,EXST,NSHIFT,NOUREC,ICOUNT,NTO,NFROM
|
|---|
| 221 | *KEEP,STATI.
|
|---|
| 222 | COMMON /STATI/ SABIN,SBBIN,INBIN,IPBIN,IKBIN,IHBIN
|
|---|
| 223 | DOUBLE PRECISION SABIN(37),SBBIN(37)
|
|---|
| 224 | INTEGER INBIN(37),IPBIN(37),IKBIN(37),IHBIN(37)
|
|---|
| 225 | *KEEP,THNVAR.
|
|---|
| 226 | COMMON /THNVAR/ STACKINT,INT_ICOUNT,THINNING
|
|---|
| 227 | INTEGER MAXICOUNT
|
|---|
| 228 | PARAMETER (MAXICOUNT=20000)
|
|---|
| 229 | DOUBLE PRECISION STACKINT(MAXICOUNT,13)
|
|---|
| 230 | INTEGER INT_ICOUNT
|
|---|
| 231 | LOGICAL THINNING
|
|---|
| 232 | *KEEP,VERS.
|
|---|
| 233 | COMMON /VERS/ VERNUM,MVDATE,VERDAT
|
|---|
| 234 | DOUBLE PRECISION VERNUM
|
|---|
| 235 | INTEGER MVDATE
|
|---|
| 236 | CHARACTER*18 VERDAT
|
|---|
| 237 | *KEEP,CEREN1.
|
|---|
| 238 | COMMON /CEREN1/ CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD,
|
|---|
| 239 | * CERSIZ,LCERFI
|
|---|
| 240 | DOUBLE PRECISION CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD
|
|---|
| 241 | REAL CERSIZ
|
|---|
| 242 | LOGICAL LCERFI
|
|---|
| 243 | *KEEP,CEREN2.
|
|---|
| 244 | COMMON /CEREN2/ PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
|
|---|
| 245 | * DCERX,DCERY,ACERX,ACERY,
|
|---|
| 246 | * XCMAX,YCMAX,EPSX,EPSY,
|
|---|
| 247 | * DCERXI,DCERYI,FCERX,FCERY,
|
|---|
| 248 | * XSCATT,YSCATT,CERXOS,CERYOS,
|
|---|
| 249 | * NCERX,NCERY,ICERML
|
|---|
| 250 | REAL PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
|
|---|
| 251 | * DCERX,DCERY,ACERX,ACERY,
|
|---|
| 252 | * XCMAX,YCMAX,EPSX,EPSY,
|
|---|
| 253 | * DCERXI,DCERYI,FCERX,FCERY,
|
|---|
| 254 | * XSCATT,YSCATT,CERXOS(20),CERYOS(20)
|
|---|
| 255 | INTEGER NCERX,NCERY,ICERML
|
|---|
| 256 | *KEEP,CEREN3.
|
|---|
| 257 | COMMON /CEREN3/ CERCNT,DATAB2,LHCER
|
|---|
| 258 | INTEGER MAXBF2
|
|---|
| 259 | PARAMETER (MAXBF2 = 39 * 7)
|
|---|
| 260 | DOUBLE PRECISION CERCNT
|
|---|
| 261 | REAL DATAB2(MAXBF2)
|
|---|
| 262 | INTEGER LHCER
|
|---|
| 263 | *KEEP,CEREN4.
|
|---|
| 264 | COMMON /CEREN4/ NRECER
|
|---|
| 265 | INTEGER NRECER
|
|---|
| 266 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 267 | *keep,certel.
|
|---|
| 268 | common /certel/ cormxd,cord,coralp,ctpars,omega,
|
|---|
| 269 | + photn,photnp,phpt,pht,vphot,
|
|---|
| 270 | + vchi,veta,vzeta,vchim,vetam,vzetam,
|
|---|
| 271 | + lambda,mu,nu,nctels,ncph
|
|---|
| 272 | double precision cormxd,cord,coralp,ctpars(20,7),omega(20,3,3),
|
|---|
| 273 | + photn(3),photnp(3),phpt(3),pht,vphot(3),
|
|---|
| 274 | + vchi(3),veta(3),vzeta(3),vchim,vetam,vzetam,
|
|---|
| 275 | + lambda,mu,nu
|
|---|
| 276 | integer nctels,ncph(5)
|
|---|
| 277 | double precision xg,yg,zg,xgp,ygp,zgp,up,vp,wp,xpcut,ypcut,zpcut
|
|---|
| 278 | equivalence (photn(1) ,xg) ,(photn(2) ,yg) ,(photn(3) ,zg) ,
|
|---|
| 279 | + (photnp(1),xgp) ,(photnp(2),ygp) ,(photnp(3),zgp),
|
|---|
| 280 | + (phpt(1) ,xpcut),(phpt(2) ,ypcut),(phpt(3) ,zpcut),
|
|---|
| 281 | + (vphot(1) ,up) ,(vphot(2) ,vp) ,(vphot(3) ,wp)
|
|---|
| 282 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 283 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 284 | C Modificacion hecha por Aitor (5-feb-98)
|
|---|
| 285 | common /aitor/ aitoth
|
|---|
| 286 | double precision aitoth
|
|---|
| 287 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 288 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 289 | c Angles for the "spinning" of a particle around the
|
|---|
| 290 | c main axis of the CT
|
|---|
| 291 | common /spinang/ spinxi
|
|---|
| 292 | double precision spinxi
|
|---|
| 293 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 294 |
|
|---|
| 295 |
|
|---|
| 296 | *KEND.
|
|---|
| 297 |
|
|---|
| 298 | INTEGER JNBIN(37),JPBIN(37),JKBIN(37),JHBIN(37)
|
|---|
| 299 | DOUBLE PRECISION CHI2,FPARAM(6)
|
|---|
| 300 | DOUBLE PRECISION MPART2(10,25),MPHOT2(10),MPOSI2(10),MELEC2(10),
|
|---|
| 301 | * MNU2(10),MMUP2(10),MMUM2(10),MPI02(10),MPIP2(10),
|
|---|
| 302 | * MPIM2(10),MK0L2(10),MKPL2(10),MKMI2(10),
|
|---|
| 303 | * MNETR2(10),MPROT2(10),MPRTB2(10),MK0S2(10),
|
|---|
| 304 | * MHYP2(10),MNETB2(10),MDEUT2(10),MTRIT2(10),
|
|---|
| 305 | * MALPH2(10),MOTH2(10)
|
|---|
| 306 | EQUIVALENCE (MPART2(1, 1),MPHOT2(1)), (MPART2(1, 2),MPOSI2(1)),
|
|---|
| 307 | * (MPART2(1, 3),MELEC2(1)), (MPART2(1, 4),MNU2(1)) ,
|
|---|
| 308 | * (MPART2(1, 5),MMUP2(1)) , (MPART2(1, 6),MMUM2(1)) ,
|
|---|
| 309 | * (MPART2(1, 7),MPI02(1)) , (MPART2(1, 8),MPIP2(1)) ,
|
|---|
| 310 | * (MPART2(1, 9),MPIM2(1)) , (MPART2(1,10),MK0L2(1)) ,
|
|---|
| 311 | * (MPART2(1,11),MKPL2(1)) , (MPART2(1,12),MKMI2(1)) ,
|
|---|
| 312 | * (MPART2(1,13),MNETR2(1)), (MPART2(1,14),MPROT2(1)),
|
|---|
| 313 | * (MPART2(1,15),MPRTB2(1)), (MPART2(1,16),MK0S2(1)) ,
|
|---|
| 314 | * (MPART2(1,18),MHYP2(1)) , (MPART2(1,19),MDEUT2(1)),
|
|---|
| 315 | * (MPART2(1,20),MTRIT2(1)), (MPART2(1,21),MALPH2(1)),
|
|---|
| 316 | * (MPART2(1,22),MOTH2 (1)), (MPART2(1,25),MNETB2(1))
|
|---|
| 317 | C VARIABLES BEING USED FOR RUNTIME
|
|---|
| 318 | REAL TDIFF
|
|---|
| 319 | INTEGER ILEFTA,ILEFTB,TIME
|
|---|
| 320 | EXTERNAL TIME
|
|---|
| 321 |
|
|---|
| 322 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 323 | double precision ctdiams(20)
|
|---|
| 324 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 325 | double precision theprim, phiprim
|
|---|
| 326 | double precision spinthe, spinphi
|
|---|
| 327 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 328 |
|
|---|
| 329 |
|
|---|
| 330 | C-----------------------------------------------------------------------
|
|---|
| 331 |
|
|---|
| 332 |
|
|---|
| 333 | CERELE = 0.D0
|
|---|
| 334 | CERHAD = 0.D0
|
|---|
| 335 | NRECER = 0
|
|---|
| 336 | C INITIALIZE AND READ RUN STEERING CARDS
|
|---|
| 337 | CALL START
|
|---|
| 338 |
|
|---|
| 339 | IF ( CERSIZ .LE. 0. ) THEN
|
|---|
| 340 | ICRSIZ = 0
|
|---|
| 341 | ELSE
|
|---|
| 342 | ICRSIZ = 1
|
|---|
| 343 | ENDIF
|
|---|
| 344 |
|
|---|
| 345 | C RESET COUNTER FOR WORDS WRITTEN TO TAPE
|
|---|
| 346 | IRECOR = 0
|
|---|
| 347 |
|
|---|
| 348 | C RESET COUNTER FOR AVERAGE HIGHT OF 1ST INTERACTION
|
|---|
| 349 | CHISUM = 0.D0
|
|---|
| 350 | CHISM2 = 0.D0
|
|---|
| 351 |
|
|---|
| 352 | C SET ARRAYS FOR SCALES OF KINETIC ENERGY-INTERACTION TABLE
|
|---|
| 353 | SABIN(1) = 0.D0
|
|---|
| 354 | SBBIN(1) = 0.1D0
|
|---|
| 355 | DO 13 J = 2,37
|
|---|
| 356 | SABIN(J) = 10.D0**((J-4.D0)/3.D0)
|
|---|
| 357 | SBBIN(J) = 10.D0**((J-3.D0)/3.D0)
|
|---|
| 358 | 13 CONTINUE
|
|---|
| 359 |
|
|---|
| 360 | C CHECK AND SET PRIMARY PARAMETERS
|
|---|
| 361 | CALL INPRM
|
|---|
| 362 |
|
|---|
| 363 | do 161 i=1,nctels
|
|---|
| 364 | ctdiams(i) = ctpars(i,ctdiam)
|
|---|
| 365 | 161 continue
|
|---|
| 366 |
|
|---|
| 367 | C INITIALIZE NKG ROUTINES
|
|---|
| 368 | CALL ININKG
|
|---|
| 369 |
|
|---|
| 370 | C RESET COUNTERS FOR NUCLEON, PION AND KAON TABLE FOR ALL SHOWERS
|
|---|
| 371 | C RESET ENERGY-MULTIPLICITY & ENERGY-ELASTICITY MATRIX FOR ALL SHOWERS
|
|---|
| 372 | DO 17 J = 1,37
|
|---|
| 373 | JNBIN(J) = 0
|
|---|
| 374 | JPBIN(J) = 0
|
|---|
| 375 | JKBIN(J) = 0
|
|---|
| 376 | JHBIN(J) = 0
|
|---|
| 377 | ELMEAA(J) = 0.D0
|
|---|
| 378 | DO 17 L = 1,13
|
|---|
| 379 | MULTOT(J,L) = 0
|
|---|
| 380 | IELDPA(J,L) = 0
|
|---|
| 381 | 17 CONTINUE
|
|---|
| 382 |
|
|---|
| 383 | C RESET OTHER ARRAYS FOR ALL SHOWERS
|
|---|
| 384 | DO 99 J = 1,20
|
|---|
| 385 | IELNU(J) = 0
|
|---|
| 386 | IELPI(J) = 0
|
|---|
| 387 | IELIS(J) = 0
|
|---|
| 388 | IELHM(J) = 0
|
|---|
| 389 | 99 CONTINUE
|
|---|
| 390 |
|
|---|
| 391 | C RESET ARRAYS FOR INTERACTION LENGTH STATISTICS
|
|---|
| 392 | DO 90 J = 1,124
|
|---|
| 393 | IHYCHI(J) = 0
|
|---|
| 394 | IKACHI(J) = 0
|
|---|
| 395 | IMUCHI(J) = 0
|
|---|
| 396 | INUCHI(J) = 0
|
|---|
| 397 | IPICHI(J) = 0
|
|---|
| 398 | INNCHI(J) = 0
|
|---|
| 399 | 90 CONTINUE
|
|---|
| 400 |
|
|---|
| 401 | C RESET ARRAY FOR MEAN VALUES AND STANDARD DEVIATION
|
|---|
| 402 | DO 477 K = 1,25
|
|---|
| 403 | DO 477 J = 1,10
|
|---|
| 404 | MPARTO(J,K) = 0.D0
|
|---|
| 405 | MPART2(J,K) = 0.D0
|
|---|
| 406 | 477 CONTINUE
|
|---|
| 407 |
|
|---|
| 408 | C RESET ARRAYS FOR LONGITUDINAL DISTRIBUTION
|
|---|
| 409 | IF ( LLONGI ) THEN
|
|---|
| 410 | DO 478 K = 1,9
|
|---|
| 411 | DO 4781 J = 0,NSTEP
|
|---|
| 412 | APLONG(J,K) = 0.D0
|
|---|
| 413 | SPLONG(J,K) = 0.D0
|
|---|
| 414 | 4781 CONTINUE
|
|---|
| 415 | 478 CONTINUE
|
|---|
| 416 | ENDIF
|
|---|
| 417 |
|
|---|
| 418 | C STEERING OF PRINTOUT OF RANDOM GENERATOR SEEDS
|
|---|
| 419 | IPROUT = MIN(100,NSHOW/20)
|
|---|
| 420 | IPROUT = MAX(1,IPROUT)
|
|---|
| 421 |
|
|---|
| 422 | C TIME AT BEGINNING
|
|---|
| 423 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 424 | c ILEFTA = TIME()
|
|---|
| 425 | ILEFTA = 0
|
|---|
| 426 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 427 | C
|
|---|
| 428 | CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 429 | C
|
|---|
| 430 | C Modified by C. Bigongiari 2001 Jan 16
|
|---|
| 431 | C
|
|---|
| 432 | C
|
|---|
| 433 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 434 | print *,'JCIO::========================================'
|
|---|
| 435 | print *,'JCIO:: Initializing JCIO system for advanced'
|
|---|
| 436 | print *,'JCIO:: saving of data.'
|
|---|
| 437 | print *,'JCIO::========================================'
|
|---|
| 438 | C
|
|---|
| 439 | Cc- initialize jcio system
|
|---|
| 440 | C
|
|---|
| 441 | call jcinitio(dsn,nrrun)
|
|---|
| 442 | Cc- create file run######
|
|---|
| 443 | C call jcstartrun(runh)
|
|---|
| 444 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 445 | C
|
|---|
| 446 | C- Modified JCSTARTRUN creates cer###### and dat###### files !
|
|---|
| 447 | C
|
|---|
| 448 | C ###### is the RUN number !
|
|---|
| 449 | C
|
|---|
| 450 |
|
|---|
| 451 | call jcstartrun(RUNH)
|
|---|
| 452 |
|
|---|
| 453 | C
|
|---|
| 454 | C- write Run Header on cer and dat files
|
|---|
| 455 | C
|
|---|
| 456 | CALL TOBUF(RUNH,0)
|
|---|
| 457 | IF ( LCERFI ) CALL TOBUFC(RUNH,0)
|
|---|
| 458 |
|
|---|
| 459 | C
|
|---|
| 460 | CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 461 | C
|
|---|
| 462 | C-----------------------------------------------------------------------
|
|---|
| 463 | C LOOP OVER SHOWERS
|
|---|
| 464 |
|
|---|
| 465 | DO 2 ISHW = 1,NSHOW
|
|---|
| 466 |
|
|---|
| 467 |
|
|---|
| 468 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 469 | c Next block of code has been modified, and comes from INPRM
|
|---|
| 470 | c----------------------------------------------------------------------
|
|---|
| 471 | C SCATTERING OF CENTER OF CHERENKOV ARRAY RELATIVE TO SHOWER AXIS
|
|---|
| 472 | c>> Actually, XSCATT and YSCATT should be RminSCATT and RmaxSCATT
|
|---|
| 473 | ICERML = MIN(MAX(ICERML,1),20)
|
|---|
| 474 | XSCATT = ABS(XSCATT)
|
|---|
| 475 | YSCATT = ABS(YSCATT)
|
|---|
| 476 | WRITE(MONIOU,5225)ICERML,XSCATT,YSCATT
|
|---|
| 477 | 5225 FORMAT(' ** USING EACH SHOWER SEVERAL TIMES:'/
|
|---|
| 478 | + ' USE EACH EVENT ',I2,' TIMES'/
|
|---|
| 479 | + ' THE EVENTS ARE SCATTERED RANDOMLY IN A SECTOR OF RADII:'/
|
|---|
| 480 | + ' Rmin = ',F10.0,' Rmax = ',F10.0)
|
|---|
| 481 | DO 4438 I=1,ICERML
|
|---|
| 482 | 5226 CALL RMMAR( RD,2,3 )
|
|---|
| 483 | CERXOS(I) = 2.0*YSCATT*(RD(1)-0.5)
|
|---|
| 484 | CERYOS(I) = 2.0*YSCATT*(RD(2)-0.5)
|
|---|
| 485 | R=SQRT(CERXOS(I)**2+CERYOS(I)**2)
|
|---|
| 486 | IF ((R.LT.XSCATT).OR.(R.GT.YSCATT)) GOTO 5226
|
|---|
| 487 | WRITE(MONIOU,4437) I,CERXOS(I),CERYOS(I)
|
|---|
| 488 | 4437 FORMAT(' CORE OF EVENT ',I2,' AT ',2F12.2)
|
|---|
| 489 | 4438 CONTINUE
|
|---|
| 490 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 491 |
|
|---|
| 492 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 493 | c Next block of code comes from INPRM
|
|---|
| 494 | c----------------------------------------------------------------------
|
|---|
| 495 | EVTH(98) = FLOAT(ICERML)
|
|---|
| 496 | DO 480 I=1,20
|
|---|
| 497 | EVTH( 98+I) = CERXOS(I)
|
|---|
| 498 | EVTH(118+I) = CERYOS(I)
|
|---|
| 499 | 480 CONTINUE
|
|---|
| 500 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 501 |
|
|---|
| 502 | SHOWNO = SHOWNO + 1
|
|---|
| 503 | I = ISHW
|
|---|
| 504 | IF ( ISHW .LE. MAXPRT ) THEN
|
|---|
| 505 | FPRINT = .TRUE.
|
|---|
| 506 | ELSE
|
|---|
| 507 | FPRINT = .FALSE.
|
|---|
| 508 | ENDIF
|
|---|
| 509 | CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 510 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 511 | Cc Create cer######,dat######,sta###### files
|
|---|
| 512 | Cc------------------------------------------------------------
|
|---|
| 513 | C call jcnewshower
|
|---|
| 514 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 515 | CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 516 |
|
|---|
| 517 | C RESET COUNTERS
|
|---|
| 518 | DO 447 K = 1,25
|
|---|
| 519 | DO 447 J = 1,10
|
|---|
| 520 | NPARTO(J,K) = 0.D0
|
|---|
| 521 | 447 CONTINUE
|
|---|
| 522 | MUOND = 0.D0
|
|---|
| 523 |
|
|---|
| 524 | C RESET ARRAY FOR LONGITUDINAL DISTRIBUTION PER SHOWER
|
|---|
| 525 | IF ( LLONGI ) THEN
|
|---|
| 526 | DO 479 K = 1,9
|
|---|
| 527 | DO 4791 J = 0,NSTEP
|
|---|
| 528 | PLONG(J,K) = 0.D0
|
|---|
| 529 | 4791 CONTINUE
|
|---|
| 530 | 479 CONTINUE
|
|---|
| 531 | ENDIF
|
|---|
| 532 |
|
|---|
| 533 | NRECS = 0
|
|---|
| 534 | NBLKS = 0
|
|---|
| 535 | DO 922 KKK = 1,10
|
|---|
| 536 | AVNREJ(KKK) = 0.D0
|
|---|
| 537 | 922 CONTINUE
|
|---|
| 538 | IRESPAR = 0
|
|---|
| 539 |
|
|---|
| 540 |
|
|---|
| 541 | C FIRST INTERACTION DATA
|
|---|
| 542 | FIRSTI = .TRUE.
|
|---|
| 543 | IFINET = 0
|
|---|
| 544 | IFINNU = 0
|
|---|
| 545 | IFINKA = 0
|
|---|
| 546 | IFINPI = 0
|
|---|
| 547 | IFINHY = 0
|
|---|
| 548 | ELAST = 0.D0
|
|---|
| 549 |
|
|---|
| 550 | C RESET COUNTERS FOR NUCLEON, PION AND KAON TABLE FOR SHOWER
|
|---|
| 551 | C RESET ENERGY-MULTIPLICITY & ENERGY-ELASTICITY MATRIX FOR SHOWER
|
|---|
| 552 | DO 11 J = 1,37
|
|---|
| 553 | INBIN(J) = 0
|
|---|
| 554 | IPBIN(J) = 0
|
|---|
| 555 | IKBIN(J) = 0
|
|---|
| 556 | IHBIN(J) = 0
|
|---|
| 557 | ELMEAN(J) = 0.D0
|
|---|
| 558 | DO 11 L = 1,13
|
|---|
| 559 | MULTMA(J,L) = 0
|
|---|
| 560 | IELDPM(J,L) = 0
|
|---|
| 561 | 11 CONTINUE
|
|---|
| 562 |
|
|---|
| 563 | DO 12 J = 1,10
|
|---|
| 564 | PBAL(J) = 0.D0
|
|---|
| 565 | EBAL(J) = 0.D0
|
|---|
| 566 | 12 CONTINUE
|
|---|
| 567 |
|
|---|
| 568 | C INITIALIZE PARTICLE STACK
|
|---|
| 569 | CALL ISTACK
|
|---|
| 570 | C RESET STACKINT
|
|---|
| 571 | DO J=1,MAXICOUNT
|
|---|
| 572 | DO K=1,MAXLEN
|
|---|
| 573 | STACKINT(J,K) = 0.D0
|
|---|
| 574 | ENDDO
|
|---|
| 575 | ENDDO
|
|---|
| 576 |
|
|---|
| 577 | C INITIALIZE EVENT HEADER AND END FOR EACH EVENT
|
|---|
| 578 | DO 2123 L = 2,43
|
|---|
| 579 | EVTH(L) = 0.
|
|---|
| 580 | 2123 CONTINUE
|
|---|
| 581 | DO 123 L = 2,MAXBUF
|
|---|
| 582 | EVTE(L) = 0.
|
|---|
| 583 | 123 CONTINUE
|
|---|
| 584 |
|
|---|
| 585 | C SHOWER BEGIN PRINTOUT
|
|---|
| 586 | IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,105) SHOWNO
|
|---|
| 587 | 105 FORMAT ('1',10('='),' SHOWER NO ',I10,' ',47('=')/)
|
|---|
| 588 |
|
|---|
| 589 | C RANDOM GENERATOR STATUS AT BEGINNING OF SHOWER CALCULATION
|
|---|
| 590 | EVTH(13) = NSEQ
|
|---|
| 591 | DO 45 L = 1,NSEQ
|
|---|
| 592 | CALL RMMAQ( ISEED(1,L), L, 'R' )
|
|---|
| 593 | C SEED
|
|---|
| 594 | EVTH(11+L*3) = ISEED(1,L)
|
|---|
| 595 | C NUMBER OF CALLS
|
|---|
| 596 | EVTH(12+L*3) = MOD ( ISEED(2,L), 1000000 )
|
|---|
| 597 | C NUMBER OF MILLIONS
|
|---|
| 598 | EVTH(13+L*3) = ISEED(3,L)*1000 + INT( ISEED(2,L)/1000000 )
|
|---|
| 599 | 45 CONTINUE
|
|---|
| 600 | IF ( FPRINT .OR. DEBUG .OR. MOD(ISHW-1,IPROUT).EQ.0 ) THEN
|
|---|
| 601 | CALL PRTIME(TTIME)
|
|---|
| 602 | WRITE(MONIOU,158) SHOWNO,(L,(ISEED(J,L),J=1,3),L=1,NSEQ)
|
|---|
| 603 | 158 FORMAT(' AND RANDOM NUMBER GENERATOR AT BEGIN OF EVENT :',I8,
|
|---|
| 604 | * /,(' SEQUENCE = ',I2,' SEED = ',I9 ,' CALLS = ',I9,
|
|---|
| 605 | * ' BILLIONS = ',I9))
|
|---|
| 606 | ENDIF
|
|---|
| 607 | C RESET KNOR
|
|---|
| 608 | KNOR = .TRUE.
|
|---|
| 609 |
|
|---|
| 610 | C GET FULL RANDOM GENERATOR STATUS (103 WORDS PER SEQUENCE)
|
|---|
| 611 | CC DO 495 L = 1,NSEQ
|
|---|
| 612 | CC CALL RMMAQ( ISEED(1,L), L, 'RV' )
|
|---|
| 613 | CC WRITE(MONIOU,658) L,(ISEED(J,L),J=1,103)
|
|---|
| 614 | CC658 FORMAT ( ' FULL RANDOM NUMBER GENERATOR STATUS ',
|
|---|
| 615 | CC * 'FOR SEQUENCE ',I2,/(' ',10I11))
|
|---|
| 616 | CC495 CONTINUE
|
|---|
| 617 |
|
|---|
| 618 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 619 | c>> *** ATENTION *** ATENTION *** ATENTION *** ATENTION *** ATENTION >>
|
|---|
| 620 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 621 | c>> >>
|
|---|
| 622 | c>> In the next block (between this ATENTION comments) CORSIKA makes >>
|
|---|
| 623 | c>> three things, in this order: >>
|
|---|
| 624 | c>> >>
|
|---|
| 625 | c>> i. Set ANGLES OF INCIDENCE (different distributions of theta >>
|
|---|
| 626 | c>> for gammas -flat- and hadrons -standard. >>
|
|---|
| 627 | c>> ii. Set HEIGHT for start at THICK0 (normally = 0 => 112.8 Km) >>
|
|---|
| 628 | c>> iii. Set ENERGY of the primary. >>
|
|---|
| 629 | c>> >>
|
|---|
| 630 | c>> (The original order was ii., iii. and i.) >>
|
|---|
| 631 | c>> >>
|
|---|
| 632 | c>> JCG Wed Sep 21 10:49:14 MET DST 1998 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 633 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 634 |
|
|---|
| 635 |
|
|---|
| 636 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 637 | C GET PRIMARY ANGLES OF INCIDENCE
|
|---|
| 638 | IF ( FIXINC ) THEN
|
|---|
| 639 |
|
|---|
| 640 | THETAP = THETPR(1)
|
|---|
| 641 | PHIP = PHIPR(1)
|
|---|
| 642 | PRMPAR(3) = COS(THETAP)
|
|---|
| 643 |
|
|---|
| 644 | ELSE
|
|---|
| 645 |
|
|---|
| 646 | if ( prmpar(1).eq.1 ) then
|
|---|
| 647 |
|
|---|
| 648 | C>> GAMMAS >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 649 | C NOTE!! We will use a FLAT distribution for THETA:
|
|---|
| 650 | C Then, next block (original block) must be commented.
|
|---|
| 651 | c The modificated code follows this block
|
|---|
| 652 | c
|
|---|
| 653 | CALL RMMAR( RD,3,1 )
|
|---|
| 654 | CT1 = THETPR(1)
|
|---|
| 655 | CT2 = THETPR(2)
|
|---|
| 656 | THETAP = RD(2)*(CT2 - CT1) + CT1
|
|---|
| 657 | CTT = COS(THETAP)
|
|---|
| 658 | PRMPAR(3) = CTT
|
|---|
| 659 |
|
|---|
| 660 | else
|
|---|
| 661 |
|
|---|
| 662 | C>> HADRONS AND ELECTRONS (AND ANY OTHER BUT GAMMAS) >>>>>>>>>>>>>>>>
|
|---|
| 663 | c Choose angles at random with equal flux for all directions
|
|---|
| 664 | c with horizontal detector array (see: O.C. Allkofer & P.K.F. Grieder,
|
|---|
| 665 | c Cosmic Rays on Earth, in: Physics Data 25/1, H.Behrens & G.Ebel Ed.,
|
|---|
| 666 | c (Fachinformationszentrum Karlsruhe, Germany, 1983) chpt. 1.1.2)
|
|---|
| 667 | c
|
|---|
| 668 | CALL RMMAR( RD,3,1 )
|
|---|
| 669 | CT1 = SIN(THETPR(1))**2
|
|---|
| 670 | CT2 = SIN(THETPR(2))**2
|
|---|
| 671 | CTT = SQRT( 1.D0 - RD(2)*(CT2 - CT1) - CT1 )
|
|---|
| 672 | PRMPAR(3) = CTT
|
|---|
| 673 | THETAP = ACOS(CTT)
|
|---|
| 674 |
|
|---|
| 675 | endif
|
|---|
| 676 |
|
|---|
| 677 | PHIP = RD(1) * ( PHIPR(2) - PHIPR(1) ) + PHIPR(1)
|
|---|
| 678 |
|
|---|
| 679 | ENDIF
|
|---|
| 680 | PRMPAR(4) = PHIP
|
|---|
| 681 |
|
|---|
| 682 |
|
|---|
| 683 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 684 | C DEFINE HEIGHT FOR START AT THICK0 (IN G/CM**2) (112.8 KM FOR THICK0=0)
|
|---|
| 685 | PRMPAR(5) = HEIGH(THICK0)
|
|---|
| 686 |
|
|---|
| 687 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 688 | C GET PRIMARY ENERGY INTO PRMPAR(2)
|
|---|
| 689 | IF ( ISPEC .EQ. 0 ) THEN
|
|---|
| 690 | PRMPAR(2) = LLIMIT
|
|---|
| 691 | ELSE
|
|---|
| 692 | CALL RMMAR( RD,1,1 )
|
|---|
| 693 | IF ( PSLOPE .NE. -1.D0 ) THEN
|
|---|
| 694 | PRMPAR(2) = ( RD(1)*UL + ( 1.D0-RD(1) )*LL )**SLEX
|
|---|
| 695 | ELSE
|
|---|
| 696 | PRMPAR(2) = LLIMIT * LL**RD(1)
|
|---|
| 697 | ENDIF
|
|---|
| 698 | ENDIF
|
|---|
| 699 |
|
|---|
| 700 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 701 | c>> *** ATENTION *** ATENTION *** ATENTION *** ATENTION *** ATENTION >>
|
|---|
| 702 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 703 |
|
|---|
| 704 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 705 | c>> Modification: this is no longer needed >>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 706 | c>> (Superseeded by Sphere algorithm, see cerenkov.f) >>>>>>>>>>>>>>>>>
|
|---|
| 707 | c>> JCG Wed Sep 21 10:49:14 MET DST 1998 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 708 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 709 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 710 | cc Btw, we now modify the "shadow area" of the telescopes,
|
|---|
| 711 | cc to cover the angle theta.
|
|---|
| 712 | c do 160 i=1,nctels
|
|---|
| 713 | c ctpars(i,ctdiam) = ctdiams(i)/cos(thetap)
|
|---|
| 714 | c write (MONIOU,*)
|
|---|
| 715 | c * 'New region for CT',i,' = ',ctpars(i,ctdiam)
|
|---|
| 716 | c 160 continue
|
|---|
| 717 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 718 |
|
|---|
| 719 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 720 | c Here we calculate the angles spinphi and spinthe, which are the
|
|---|
| 721 | c phi and theta angles of the particle, with respect to the direction
|
|---|
| 722 | c where the CT is pointing to. spinthe is the angular displacement
|
|---|
| 723 | c of the new direction respect to the original (CT); spinphi=0 means
|
|---|
| 724 | c that the new direction is towards the zenith, spinphi=+-180 means
|
|---|
| 725 | c towards the horizont.
|
|---|
| 726 | c See the document "simulation.tex"
|
|---|
| 727 |
|
|---|
| 728 | c First, save the "CT" orientation
|
|---|
| 729 | c (moved from a couple of lines below, marked with [*])
|
|---|
| 730 | EVTH(11) = THETAP
|
|---|
| 731 | EVTH(12) = PHIP
|
|---|
| 732 |
|
|---|
| 733 | CALL RMMAR( RD,3,1 )
|
|---|
| 734 |
|
|---|
| 735 | c Then, calculate the new direction relative to the CT direction
|
|---|
| 736 | spinphi = RD(1)*PI
|
|---|
| 737 | spinthe = RD(2)*spinxi*pi/180
|
|---|
| 738 |
|
|---|
| 739 | c And then, RE-calculate the GLOBAL direction in CORSIKA
|
|---|
| 740 | c We use formulae for spherical triangles
|
|---|
| 741 |
|
|---|
| 742 | if ( THETAP .EQ. 0.0 ) then
|
|---|
| 743 | theprim = spinthe
|
|---|
| 744 | phiprim = spinphi
|
|---|
| 745 | elseif ( spinthe .eq. 0.0 ) then
|
|---|
| 746 | theprim = THETAP
|
|---|
| 747 | phiprim = 0.0
|
|---|
| 748 | else
|
|---|
| 749 | theprim = acos( cos(THETAP)*cos(spinthe)+
|
|---|
| 750 | $ sin(THETAP)*sin(spinthe)*cos(spinphi) )
|
|---|
| 751 | phiprim = asin( sin(spinthe)*sin(spinphi)/sin(theprim) )
|
|---|
| 752 | endif
|
|---|
| 753 | THETAP = theprim
|
|---|
| 754 | EVTH(140) = spinthe
|
|---|
| 755 | if (RD(3).gt.0.5) then
|
|---|
| 756 | PHIP = PHIP - phiprim
|
|---|
| 757 | EVTH(141) = -spinphi
|
|---|
| 758 | else
|
|---|
| 759 | PHIP = PHIP + phiprim
|
|---|
| 760 | EVTH(141) = spinphi
|
|---|
| 761 | endif
|
|---|
| 762 | PRMPAR(3) = COS(THETAP)
|
|---|
| 763 | PRMPAR(4) = PHIP
|
|---|
| 764 |
|
|---|
| 765 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 766 |
|
|---|
| 767 |
|
|---|
| 768 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 769 | c>> Modification (HZA trick) cancelled >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 770 | c>> JCG Wed Sep 21 10:49:14 MET DST 1998 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 771 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 772 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 773 | C Modificacion hecha por Aitor
|
|---|
| 774 | c aitoth = THETAP
|
|---|
| 775 | C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 776 |
|
|---|
| 777 | C CALCULATE COORDINATE CORRECTION FOR TOP OF ATMOSPHERE
|
|---|
| 778 | CALL COORIN( PRMPAR(5) )
|
|---|
| 779 |
|
|---|
| 780 | C COUNTER FOR PARTICLE OUTPUT
|
|---|
| 781 | LH = 0
|
|---|
| 782 | C COUNTER FOR CERENKOV OUTPUT
|
|---|
| 783 | IF ( LCERFI ) LHCER = 0
|
|---|
| 784 | C CALCULATE BUNCH SIZE FOR CERENKOV PHOTONS IF NOT SET IN DATAC
|
|---|
| 785 | IF ( ICRSIZ .EQ. 0 ) THEN
|
|---|
| 786 | CALL GETBUS( NINT(PRMPAR(1)),SNGL(PRMPAR(2)),SNGL(PRMPAR(3)),
|
|---|
| 787 | * CERSIZ )
|
|---|
| 788 | WRITE(MONIOU,*)'CERENKOV BUNCH SIZE IS CALCULATED TO=',CERSIZ
|
|---|
| 789 | ENDIF
|
|---|
| 790 | C GET GAMMA FACTOR FROM ENERGY
|
|---|
| 791 | C FOR GAMMAS PRMPAR(2) STAYS = ENERGY
|
|---|
| 792 | IF ( PRMPAR(1) .NE. 1.D0 )
|
|---|
| 793 | * PRMPAR(2) = PRMPAR(2) / PAMA(NINT(PRMPAR(1)))
|
|---|
| 794 |
|
|---|
| 795 | C SET PRIMARY TO CURRENT PARTICLE
|
|---|
| 796 | DO 3 J = 1,8
|
|---|
| 797 | CURPAR(J) = PRMPAR(J)
|
|---|
| 798 | NCOUN(J) = 0
|
|---|
| 799 | 3 CONTINUE
|
|---|
| 800 | C SET WEIGHT
|
|---|
| 801 |
|
|---|
| 802 | C CALCULATE FIRST INTERACTION POINT IF HADRONIC
|
|---|
| 803 | GEN = 0.D0
|
|---|
| 804 |
|
|---|
| 805 | H = HEIGH(THICK0)
|
|---|
| 806 | CALL BOX2
|
|---|
| 807 | IF ( FIX1I ) THEN
|
|---|
| 808 | CHI = THICK(FIXHEI) / PRMPAR(3)
|
|---|
| 809 | H = FIXHEI
|
|---|
| 810 | FDECAY = .FALSE.
|
|---|
| 811 | ELSE
|
|---|
| 812 | H = HEIGH ( CHI*PRMPAR(3) + THICK0 )
|
|---|
| 813 | ENDIF
|
|---|
| 814 | CHISUM = CHISUM + CHI
|
|---|
| 815 | CHISM2 = CHISM2 + CHI**2
|
|---|
| 816 | ALEVEL = H
|
|---|
| 817 | C INITIALIZE COORDINATE CORRECTIONS FOR HADRONIC PRIMARIES
|
|---|
| 818 | C FOR EM PRIMARIES IT IS DONE IN EGS
|
|---|
| 819 | HH = MAX( H, 0.D0 )
|
|---|
| 820 | IF ( CURPAR(1) .GT. 3.D0 ) CALL COORIN( HH )
|
|---|
| 821 |
|
|---|
| 822 | IF ( FMUADD ) THEN
|
|---|
| 823 | IF ( CURPAR(1) .EQ. 5 .OR. CURPAR(1) .EQ. 6) THEN
|
|---|
| 824 | DO J = 1,MAXLEN
|
|---|
| 825 | AMUPAR(J) = CURPAR(J)
|
|---|
| 826 | ENDDO
|
|---|
| 827 | AMUPAR(5) = PRMPAR(5)
|
|---|
| 828 | IF(DEBUG)WRITE(MDEBUG,*)'MAIN : MUON STORED IN AMUPAR'
|
|---|
| 829 | FMUORG = .TRUE.
|
|---|
| 830 | ENDIF
|
|---|
| 831 | ENDIF
|
|---|
| 832 |
|
|---|
| 833 | C SET TARGET FLAG IF SELECTED FOR FIRST INTERACTION
|
|---|
| 834 | IF ( N1STTR .GT. 0 ) THEN
|
|---|
| 835 | FIXTAR = .TRUE.
|
|---|
| 836 | FDECAY = .FALSE.
|
|---|
| 837 | EVTH(6) = REAL(N1STTR)
|
|---|
| 838 | ELSE
|
|---|
| 839 | FIXTAR = .FALSE.
|
|---|
| 840 | EVTH(6) = 0.0
|
|---|
| 841 | ENDIF
|
|---|
| 842 |
|
|---|
| 843 | C INITIALIZE ARRAYS FOR NKG FOR EACH SHOWER
|
|---|
| 844 | IF ( FNKG ) CALL STANKG
|
|---|
| 845 |
|
|---|
| 846 | C STORE FIRST PARTICLE IN HEADER AND PRINT IT OUT
|
|---|
| 847 | EVTH( 2) = REAL(SHOWNO)
|
|---|
| 848 | EVTH( 3) = CURPAR(1)
|
|---|
| 849 | IF ( CURPAR(1) .EQ. 1.D0 ) THEN
|
|---|
| 850 | C PRIMARY ENERGY FOR PHOTONS
|
|---|
| 851 | E00 = GAMMA
|
|---|
| 852 | E00PN = GAMMA
|
|---|
| 853 | INUCL = 1
|
|---|
| 854 | ELSE
|
|---|
| 855 | E00 = GAMMA * PAMA(NINT(CURPAR(1)))
|
|---|
| 856 | INUCL = INT(MAX(1.D0,CURPAR(1)/100.D0))
|
|---|
| 857 | E00PN = E00 / INUCL
|
|---|
| 858 | ENDIF
|
|---|
| 859 | EVTH(147) = 0.
|
|---|
| 860 |
|
|---|
| 861 | IF ( FEGS ) THEN
|
|---|
| 862 | C PARAMETER FOR ELECTRON AND PHOTON REJECT (CONVERT ENERGY TO MEV)
|
|---|
| 863 | EONCUT = .5E-9*SQRT(E00*1000.D0)
|
|---|
| 864 | CUTLN = LOG(EONCUT)
|
|---|
| 865 | ENDIF
|
|---|
| 866 | EVTH( 4) = E00
|
|---|
| 867 | EVTH( 5) = THICK0
|
|---|
| 868 | EVTH( 7) = H
|
|---|
| 869 | PTOT0 = SQRT( E00**2 - PAMA(NINT(CURPAR(1)))**2 )
|
|---|
| 870 | PTOT0N = PTOT0 / INUCL
|
|---|
| 871 | ST = SQRT(1.D0-COSTHE**2)
|
|---|
| 872 | EVTH( 8) = PTOT0 * ST * COS(PHI)
|
|---|
| 873 | EVTH( 9) = PTOT0 * ST * SIN(PHI)
|
|---|
| 874 | EVTH(10) = PTOT0 * COSTHE
|
|---|
| 875 |
|
|---|
| 876 |
|
|---|
| 877 | c
|
|---|
| 878 | c [*] one block from here sent above
|
|---|
| 879 | c
|
|---|
| 880 |
|
|---|
| 881 | EVTH(85) = CERSIZ
|
|---|
| 882 |
|
|---|
| 883 | IF ( CURPAR(1) .GT. 3.D0 ) THEN
|
|---|
| 884 | IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,102) (CURPAR(J),J = 1,8)
|
|---|
| 885 | 102 FORMAT (/' PRIMARY PARAMETERS AT FIRST INTERACTION POINT'/
|
|---|
| 886 | * 16X,1P,8E10.3)
|
|---|
| 887 | ELSE
|
|---|
| 888 | IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,132)
|
|---|
| 889 | 132 FORMAT (/' PRIMARY PARTICLE IS ELECTROMAGNETIC')
|
|---|
| 890 | ENDIF
|
|---|
| 891 |
|
|---|
| 892 | C WRITE EVENT HEADER INTO BUFFER
|
|---|
| 893 | C FOR EM PARTICLES EVTH IS WRITTEN TO BUFFER IN EGS (IF ACTIVE)
|
|---|
| 894 | IF ( EVTH(3) .GT. 3.0 .OR. .NOT. FEGS ) THEN
|
|---|
| 895 | CALL TOBUF ( EVTH,0 )
|
|---|
| 896 | IF ( LCERFI ) CALL TOBUFC( EVTH,0 )
|
|---|
| 897 | ENDIF
|
|---|
| 898 |
|
|---|
| 899 | C PRINT HEADER FOR HIGH ENERGY PARTICLES
|
|---|
| 900 | IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,103)
|
|---|
| 901 | 103 FORMAT(/' TYPE GAMMA COSTHETA ',
|
|---|
| 902 | * ' PHI HEIGHT TIME X-CM Y-CM ',
|
|---|
| 903 | * ' GEN LEVEL E ON STACK'/)
|
|---|
| 904 | NOPART = 0
|
|---|
| 905 |
|
|---|
| 906 | IF ( CURPAR(1) .LE. 3.D0 .OR.
|
|---|
| 907 | * (CURPAR(1) .GE. 5.D0 .AND. CURPAR(1) .LE. 7.D0) ) THEN
|
|---|
| 908 | C GIVE PARTICLE TO EGS OR NKG IF ELECTROMAGNETIC
|
|---|
| 909 | C AND TAKE THEN NEXT PARTICLE FROM STACK
|
|---|
| 910 | C FLAG FOR NO PRIMARY INTERACTION IS SET FOR ALL BUT ELM. PRIMARIES
|
|---|
| 911 | IF ( CURPAR(1) .LE. 3.D0 ) THEN
|
|---|
| 912 | FNPRIM = .FALSE.
|
|---|
| 913 | ELSE
|
|---|
| 914 | FNPRIM = .TRUE.
|
|---|
| 915 | H = PRMPAR(5)
|
|---|
| 916 | ENDIF
|
|---|
| 917 | CALL BOX3
|
|---|
| 918 | IF ( FEGS ) THEN
|
|---|
| 919 | CHISUM = CHISUM + THICK(DBLE(EVTH(7)))
|
|---|
| 920 | CHISM2 = CHISM2 + THICK(DBLE(EVTH(7)))**2
|
|---|
| 921 | ENDIF
|
|---|
| 922 | FIRSTI = .FALSE.
|
|---|
| 923 | GOTO 4
|
|---|
| 924 |
|
|---|
| 925 | ELSE
|
|---|
| 926 | C HADRONIC PARTICLES
|
|---|
| 927 | FNPRIM = .TRUE.
|
|---|
| 928 |
|
|---|
| 929 | C FILL LONGITUDINAL DISTRIBUTION FOR THE PRIMARY PARTICLE
|
|---|
| 930 | C THE PARTICLE IS TRACKED FROM THICK0 DOWN TO THICK0+CHI*PRMPAR(3)
|
|---|
| 931 | C COUNT THE PARTICLES FOR THE LONGITUDINAL DEVELOPMENT
|
|---|
| 932 | IF ( LLONGI ) THEN
|
|---|
| 933 | LPCT1 = INT( THICK0 * THSTPI )
|
|---|
| 934 | LPCT2 = INT( (THICK0 + PRMPAR(3)*CHI) * THSTPI )
|
|---|
| 935 | LPCT2 = MIN(NSTEP,LPCT2)
|
|---|
| 936 | C GAMMAS, ELECTRONS AND POSITRONS ARE NOT TRANSPORTED HERE, SEE EGS
|
|---|
| 937 | C MUONS ARE TRANSPORTED IN MUTRAC
|
|---|
| 938 | C HADRONS
|
|---|
| 939 | IF ( ITYPE .GE. 7 .AND. ITYPE .LE. 41 ) THEN
|
|---|
| 940 | DO 5004 L = LPCT1,LPCT2
|
|---|
| 941 | PLONG(L,6) = PLONG(L,6) + 1.
|
|---|
| 942 | 5004 CONTINUE
|
|---|
| 943 | C CHARGED HADRONS
|
|---|
| 944 | IF ( SIGNUM(ITYPE) .NE. 0.D0 ) THEN
|
|---|
| 945 | DO 5005 L = LPCT1,LPCT2
|
|---|
| 946 | PLONG(L,7) = PLONG(L,7) + 1.
|
|---|
| 947 | 5005 CONTINUE
|
|---|
| 948 | ENDIF
|
|---|
| 949 | C NUCLEI
|
|---|
| 950 | ELSEIF ( ITYPE .GT. 100 ) THEN
|
|---|
| 951 | DO 5006 L = LPCT1,LPCT2
|
|---|
| 952 | PLONG(L,8) = PLONG(L,8) + 1.
|
|---|
| 953 | 5006 CONTINUE
|
|---|
| 954 | ENDIF
|
|---|
| 955 | ENDIF
|
|---|
| 956 |
|
|---|
| 957 | C CHECK OBSERVATION LEVEL PASSAGE AND UPDATE PARTICLE COORDINATES
|
|---|
| 958 | HNEW = H
|
|---|
| 959 | C FOR UPDATE WE NEED THE START ALTITUDE H
|
|---|
| 960 | H = HEIGH(THICK0)
|
|---|
| 961 | DO 251 J = 1,NOBSLV
|
|---|
| 962 | C JUMP INTO NORMAL PARTICLE TREATMENT FOR HADRONS
|
|---|
| 963 | IF ( HNEW .GT. OBSLEV(J) ) THEN
|
|---|
| 964 | H = HNEW
|
|---|
| 965 | GOTO 6
|
|---|
| 966 | ENDIF
|
|---|
| 967 | IF ( H .LT. OBSLEV(J) ) GOTO 251
|
|---|
| 968 | C REMEMBER NUMBER OF LEVEL FOR OUTPUT
|
|---|
| 969 | LEVL = J
|
|---|
| 970 | CALL UPDATE( OBSLEV(J),THCKOB(J),J )
|
|---|
| 971 | IF (DEBUG) WRITE(MDEBUG,256) J,IRET1,IRET2
|
|---|
| 972 | 256 FORMAT(' MAIN : LEVEL ',I5,' IRET1,2=',2I5)
|
|---|
| 973 | C IF PARTICLE IS NOT CUTTED, BRING IT TO OUTPUT
|
|---|
| 974 | IF ( IRET2 .EQ. 0 ) THEN
|
|---|
| 975 | CALL OUTPUT
|
|---|
| 976 | ENDIF
|
|---|
| 977 | 251 CONTINUE
|
|---|
| 978 | IF (DEBUG) WRITE(MDEBUG,*)
|
|---|
| 979 | * 'MAIN : PRIMARY REACHED LOWEST OBSERVATION LEVEL'
|
|---|
| 980 | GOTO 4
|
|---|
| 981 | ENDIF
|
|---|
| 982 |
|
|---|
| 983 |
|
|---|
| 984 | C-----------------------------------------------------------------------
|
|---|
| 985 | C NORMAL CYCLE
|
|---|
| 986 | 7 CONTINUE
|
|---|
| 987 |
|
|---|
| 988 | C IF ENERGY TOO SMALL TAKE NEXT PARTICLE
|
|---|
| 989 | IF ( GAMMA .LE. 1.D0 ) THEN
|
|---|
| 990 | IF ( CURPAR(1) .NE. 1.D0 ) THEN
|
|---|
| 991 | IF ( CURPAR(1).EQ.5.D0 .OR. CURPAR(1).EQ.6.D0 )
|
|---|
| 992 | * FMUORG = .FALSE.
|
|---|
| 993 | GOTO 4
|
|---|
| 994 | ENDIF
|
|---|
| 995 | C SPECIAL TREATMENT FOR PHOTONS
|
|---|
| 996 | ITYPE = 1
|
|---|
| 997 | CHI = 0.D0
|
|---|
| 998 | GOTO 5
|
|---|
| 999 | ENDIF
|
|---|
| 1000 |
|
|---|
| 1001 | C DETERMINE PLACE OF NEXT INTERACTION
|
|---|
| 1002 | CALL BOX2
|
|---|
| 1003 |
|
|---|
| 1004 | C CHECK PASSAGE THROUGH OBSERVATION LEVELS AND TRACK PARTICLES TO THE
|
|---|
| 1005 | C PLACE OF INTERACTION
|
|---|
| 1006 | 5 CONTINUE
|
|---|
| 1007 | IRET1 = 0
|
|---|
| 1008 | CALL BOX3
|
|---|
| 1009 | IF ( IRET1 .NE. 0 ) GOTO 4
|
|---|
| 1010 |
|
|---|
| 1011 | 6 CONTINUE
|
|---|
| 1012 | IRET1 = 0
|
|---|
| 1013 | MSMM = 0
|
|---|
| 1014 |
|
|---|
| 1015 | C INCREMENT PARTICLE GENERATION AND PROCESS NUCLEAR INTERACTION
|
|---|
| 1016 | GEN = GEN + 1.D0
|
|---|
| 1017 | C INITIALIZE INTERMEDIATE STACK FOR ONE REACTION
|
|---|
| 1018 | CALL TSTINI
|
|---|
| 1019 | CALL NUCINT
|
|---|
| 1020 | C TRANSFER INTERMEDIATE STACK FOR ONE REACTION
|
|---|
| 1021 | CALL TSTEND
|
|---|
| 1022 |
|
|---|
| 1023 | C ENERGY - MULTIPLICITY STATISTICS
|
|---|
| 1024 | IF ( EKINL .LE. 0.1D0 ) THEN
|
|---|
| 1025 | MEN = 1
|
|---|
| 1026 | ELSE
|
|---|
| 1027 | MEN = 4.D0 + 3.D0 * LOG10(EKINL)
|
|---|
| 1028 | MEN = MIN( MEN, 37 )
|
|---|
| 1029 | ENDIF
|
|---|
| 1030 | IF ( MSMM .LE. 1 ) THEN
|
|---|
| 1031 | MMU = 1
|
|---|
| 1032 | ELSE
|
|---|
| 1033 | MMU = 1.D0 + 3.D0 * LOG10(DBLE(MSMM))
|
|---|
| 1034 | MMU = MIN( MMU, 13 )
|
|---|
| 1035 | ENDIF
|
|---|
| 1036 | MULTMA(MEN,MMU) = MULTMA(MEN,MMU) + 1
|
|---|
| 1037 | MULTOT(MEN,MMU) = MULTOT(MEN,MMU) + 1
|
|---|
| 1038 | IF ( DEBUG ) WRITE(MDEBUG,*) 'MAIN : EKINL,MSMM=',
|
|---|
| 1039 | * SNGL(EKINL),MSMM
|
|---|
| 1040 |
|
|---|
| 1041 | IF ( IRET1 .EQ. 0 ) THEN
|
|---|
| 1042 | IF ( DEBUG ) WRITE(MDEBUG,666) (CURPAR(II),II=1,11)
|
|---|
| 1043 | 666 FORMAT(' MAIN : CURPAR=',1P,11E10.3)
|
|---|
| 1044 | GOTO 7
|
|---|
| 1045 | ENDIF
|
|---|
| 1046 |
|
|---|
| 1047 | C GET NEXT PARTICLE FROM STACK, IF IRET=1 ALL PARTICLES ARE DONE
|
|---|
| 1048 | 4 CONTINUE
|
|---|
| 1049 | IRET1 = 0
|
|---|
| 1050 | CALL FSTACK
|
|---|
| 1051 | IF ( FMUADD ) THEN
|
|---|
| 1052 | IF ( (CURPAR(1) .EQ. 5 .OR. CURPAR(1) .EQ. 6)
|
|---|
| 1053 | * .AND. IRET1 .EQ. 0 .AND. .NOT. FMUORG ) THEN
|
|---|
| 1054 | DO J = 1,MAXLEN
|
|---|
| 1055 | AMUPAR(J) = CURPAR(J)
|
|---|
| 1056 | ENDDO
|
|---|
| 1057 | IF(DEBUG)WRITE(MDEBUG,*)'MAIN : MUON STORED IN AMUPAR'
|
|---|
| 1058 | FMUORG = .TRUE.
|
|---|
| 1059 | ENDIF
|
|---|
| 1060 | ENDIF
|
|---|
| 1061 |
|
|---|
| 1062 | IF ( IRET1 .EQ. 0 ) GOTO 7
|
|---|
| 1063 |
|
|---|
| 1064 | C-----------------------------------------------------------------------
|
|---|
| 1065 | C FINISH SHOWER AND PRINT INFORMATION
|
|---|
| 1066 | CALL OUTEND
|
|---|
| 1067 |
|
|---|
| 1068 | * IF ( DEBUG ) WRITE(MDEBUG,442) NPARTO
|
|---|
| 1069 | *442 FORMAT(' MAIN : NPARTO='/(' ',10F10.0))
|
|---|
| 1070 |
|
|---|
| 1071 | IF ( FPRINT .OR. DEBUG ) THEN
|
|---|
| 1072 | IOBSLV = MIN( 5, NOBSLV )
|
|---|
| 1073 | WRITE(MONIOU,54) (K,K=1,IOBSLV)
|
|---|
| 1074 | 54 FORMAT (/' PARTICLES AT DETECTOR LEVEL :'/
|
|---|
| 1075 | * ' FOR LEVEL ', 5I13)
|
|---|
| 1076 | WRITE(MONIOU,55) (OBSLEV(K),K=1,IOBSLV)
|
|---|
| 1077 | 55 FORMAT ( ' HEIGHT IN CM ',1P, 5E13.3/)
|
|---|
| 1078 | WRITE(MONIOU,776) 'PROTONS ',(NPROTO(K),K=1,IOBSLV)
|
|---|
| 1079 | WRITE(MONIOU,776) 'ANTIPROTONS ',(NPROTB(K),K=1,IOBSLV)
|
|---|
| 1080 | WRITE(MONIOU,776) 'NEUTRONS ',(NNEUTR(K),K=1,IOBSLV)
|
|---|
| 1081 | WRITE(MONIOU,776) 'ANTINEUTRONS ',(NNEUTB(K),K=1,IOBSLV)
|
|---|
| 1082 | WRITE(MONIOU,776) 'PHOTONS ',(NPHOTO(K),K=1,IOBSLV)
|
|---|
| 1083 | WRITE(MONIOU,776) 'ELECTRONS ',(NELECT(K),K=1,IOBSLV)
|
|---|
| 1084 | WRITE(MONIOU,776) 'POSITRONS ',(NPOSIT(K),K=1,IOBSLV)
|
|---|
| 1085 | WRITE(MONIOU,776) 'NEUTRINOS ',(NNU (K),K=1,IOBSLV)
|
|---|
| 1086 | WRITE(MONIOU,776) 'MU - ',(NMUM (K),K=1,IOBSLV)
|
|---|
| 1087 | WRITE(MONIOU,776) 'MU + ',(NMUP (K),K=1,IOBSLV)
|
|---|
| 1088 | WRITE(MONIOU,776) 'PI 0 ',(NPI0 (K),K=1,IOBSLV)
|
|---|
| 1089 | WRITE(MONIOU,776) 'PI - ',(NPIM (K),K=1,IOBSLV)
|
|---|
| 1090 | WRITE(MONIOU,776) 'PI + ',(NPIP (K),K=1,IOBSLV)
|
|---|
| 1091 | WRITE(MONIOU,776) 'K0L ',(NK0L (K),K=1,IOBSLV)
|
|---|
| 1092 | WRITE(MONIOU,776) 'K0S ',(NK0S (K),K=1,IOBSLV)
|
|---|
| 1093 | WRITE(MONIOU,776) 'K - ',(NKMI (K),K=1,IOBSLV)
|
|---|
| 1094 | WRITE(MONIOU,776) 'K + ',(NKPL (K),K=1,IOBSLV)
|
|---|
| 1095 | WRITE(MONIOU,776) 'STR. BARYONS ',(NHYP (K),K=1,IOBSLV)
|
|---|
| 1096 | WRITE(MONIOU,776) 'DEUTERONS ',(NDEUT (K),K=1,IOBSLV)
|
|---|
| 1097 | WRITE(MONIOU,776) 'TRITONS ',(NTRIT (K),K=1,IOBSLV)
|
|---|
| 1098 | WRITE(MONIOU,776) 'ALPHAS ',(NALPHA(K),K=1,IOBSLV)
|
|---|
| 1099 | WRITE(MONIOU,776) 'OTHER PARTIC.',(NOTHER(K),K=1,IOBSLV)
|
|---|
| 1100 | WRITE(MONIOU,*)
|
|---|
| 1101 | WRITE(MONIOU,776) 'DECAYED MUONS',MUOND
|
|---|
| 1102 | 776 FORMAT(' NO OF ',A13, '= ',5F13.0)
|
|---|
| 1103 |
|
|---|
| 1104 | IF ( NOBSLV .GT. 5 ) THEN
|
|---|
| 1105 | IOBSLV = NOBSLV
|
|---|
| 1106 | WRITE(MONIOU,54) (K,K=6,IOBSLV)
|
|---|
| 1107 | WRITE(MONIOU,55) (OBSLEV(K),K=6,IOBSLV)
|
|---|
| 1108 | WRITE(MONIOU,776) 'PROTONS ',(NPROTO(K),K=6,IOBSLV)
|
|---|
| 1109 | WRITE(MONIOU,776) 'ANTIPROTONS ',(NPROTB(K),K=6,IOBSLV)
|
|---|
| 1110 | WRITE(MONIOU,776) 'NEUTRONS ',(NNEUTR(K),K=6,IOBSLV)
|
|---|
| 1111 | WRITE(MONIOU,776) 'ANTINEUTRONS ',(NNEUTB(K),K=6,IOBSLV)
|
|---|
| 1112 | WRITE(MONIOU,776) 'PHOTONS ',(NPHOTO(K),K=6,IOBSLV)
|
|---|
| 1113 | WRITE(MONIOU,776) 'ELECTRONS ',(NELECT(K),K=6,IOBSLV)
|
|---|
| 1114 | WRITE(MONIOU,776) 'POSITRONS ',(NPOSIT(K),K=6,IOBSLV)
|
|---|
| 1115 | WRITE(MONIOU,776) 'NEUTRINOS ',(NNU (K),K=6,IOBSLV)
|
|---|
| 1116 | WRITE(MONIOU,776) 'MU - ',(NMUM (K),K=6,IOBSLV)
|
|---|
| 1117 | WRITE(MONIOU,776) 'MU + ',(NMUP (K),K=6,IOBSLV)
|
|---|
| 1118 | WRITE(MONIOU,776) 'PI 0 ',(NPI0 (K),K=6,IOBSLV)
|
|---|
| 1119 | WRITE(MONIOU,776) 'PI - ',(NPIM (K),K=6,IOBSLV)
|
|---|
| 1120 | WRITE(MONIOU,776) 'PI + ',(NPIP (K),K=6,IOBSLV)
|
|---|
| 1121 | WRITE(MONIOU,776) 'K0L ',(NK0L (K),K=6,IOBSLV)
|
|---|
| 1122 | WRITE(MONIOU,776) 'K0S ',(NK0S (K),K=6,IOBSLV)
|
|---|
| 1123 | WRITE(MONIOU,776) 'K - ',(NKMI (K),K=6,IOBSLV)
|
|---|
| 1124 | WRITE(MONIOU,776) 'K + ',(NKPL (K),K=6,IOBSLV)
|
|---|
| 1125 | WRITE(MONIOU,776) 'STR. BARYONS ',(NHYP (K),K=6,IOBSLV)
|
|---|
| 1126 | WRITE(MONIOU,776) 'DEUTERONS ',(NDEUT (K),K=6,IOBSLV)
|
|---|
| 1127 | WRITE(MONIOU,776) 'TRITONS ',(NTRIT (K),K=6,IOBSLV)
|
|---|
| 1128 | WRITE(MONIOU,776) 'ALPHAS ',(NALPHA(K),K=6,IOBSLV)
|
|---|
| 1129 | WRITE(MONIOU,776) 'OTHER PARTIC.',(NOTHER(K),K=6,IOBSLV)
|
|---|
| 1130 | WRITE(MONIOU,*)
|
|---|
| 1131 | ENDIF
|
|---|
| 1132 | ENDIF
|
|---|
| 1133 |
|
|---|
| 1134 | C ADD UP FOR MEAN VALUES
|
|---|
| 1135 | DO 779 K = 1,25
|
|---|
| 1136 | DO 779 J = 1,10
|
|---|
| 1137 | MPARTO(J,K) = MPARTO(J,K) + NPARTO(J,K)
|
|---|
| 1138 | MPART2(J,K) = MPART2(J,K) + NPARTO(J,K)**2
|
|---|
| 1139 | 779 CONTINUE
|
|---|
| 1140 | EVTE(2) = SHOWNO
|
|---|
| 1141 | DO 335 K = 1,NOBSLV
|
|---|
| 1142 | EVTE(3) = EVTE(3) + NPHOTO(K)
|
|---|
| 1143 | EVTE(4) = EVTE(4) + NELECT(K) + NPOSIT(K)
|
|---|
| 1144 | EVTE(5) = EVTE(5) + NPROTO(K) + NPROTB(K) + NNEUTR(K) +
|
|---|
| 1145 | * NNEUTB(K) + NPI0(K) + NPIM(K) + NPIP(K) + NK0L(K) +
|
|---|
| 1146 | * NK0S(K) + NKMI(K) + NKPL(K) + NHYP(K) +
|
|---|
| 1147 | * NDEUT(K) + NTRIT(K) + NALPHA(K) + NOTHER(K)
|
|---|
| 1148 | EVTE(6) = EVTE(6) + NMUP(K) + NMUM(K)
|
|---|
| 1149 | 335 CONTINUE
|
|---|
| 1150 | EVTE(7) = NOPART
|
|---|
| 1151 |
|
|---|
| 1152 | IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,110)
|
|---|
| 1153 | * IFINNU,IFINPI,IFINET,IFINKA,IFINHY,
|
|---|
| 1154 | * IFINNU+IFINPI+IFINET+IFINKA+IFINHY,ELAST
|
|---|
| 1155 | 110 FORMAT(/' NO OF NUCLEONS PRODUCED IN FIRST INTERACTION =',I10/
|
|---|
| 1156 | * ' NO OF PIONS PRODUCED IN FIRST INTERACTION =',I10/
|
|---|
| 1157 | * ' NO OF ETAS PRODUCED IN FIRST INTERACTION =',I10/
|
|---|
| 1158 | * ' NO OF KAONS PRODUCED IN FIRST INTERACTION =',I10/
|
|---|
| 1159 | * ' NO OF S.BARYONS PRODUCED IN FIRST INTERACTION =',I10/
|
|---|
| 1160 | * ' TOTAL MULTIPLICITY OF FIRST INTERACTION =',I10/
|
|---|
| 1161 | * ' ELASTICITY OF FIRST INTERACTION =',F10.4)
|
|---|
| 1162 |
|
|---|
| 1163 | C PRINT OUT NKG RESULT FOR ONE SHOWER IF SELECTED
|
|---|
| 1164 | IF ( FNKG ) CALL AVAGE
|
|---|
| 1165 |
|
|---|
| 1166 | IF ( LLONGI ) THEN
|
|---|
| 1167 | C TREAT LONGITUDINAL DISTRIBUTIONS
|
|---|
| 1168 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1169 | c calculated here again, 'cos it's rewrite I dont know where
|
|---|
| 1170 | LPCT1 = INT( THICK0 * THSTPI )
|
|---|
| 1171 | LPCT2 = INT( (THICK0 + PRMPAR(3)*CHI) * THSTPI )
|
|---|
| 1172 | LPCT2 = MIN(NSTEP,LPCT2)
|
|---|
| 1173 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1174 | DO 980 J = LPCT1,NSTEP
|
|---|
| 1175 | C ADD ELECTRONS, POSITRONS, MUONS AND NUCLEI TO THE CHARGED PARTICLES
|
|---|
| 1176 | PLONG(J,7) = PLONG(J,7) + PLONG(J,2) + PLONG(J,3)
|
|---|
| 1177 | * + PLONG(J,4) + PLONG(J,5) + PLONG(J,8)
|
|---|
| 1178 | C ADD UP FOR MEAN VALUES OF LONGITUDINAL DISTRIBUTION
|
|---|
| 1179 | DO 979 K = 1,9
|
|---|
| 1180 | APLONG(J,K) = APLONG(J,K) + PLONG(J,K)
|
|---|
| 1181 | SPLONG(J,K) = SPLONG(J,K) + PLONG(J,K)**2
|
|---|
| 1182 | 979 CONTINUE
|
|---|
| 1183 | 980 CONTINUE
|
|---|
| 1184 |
|
|---|
| 1185 | C PRINT LONGITUDINAL DISTRIBUTIONS PER SHOWER
|
|---|
| 1186 | IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,910) THSTEP,
|
|---|
| 1187 | * 'GAMMAS','POSITRONS','ELECTRONS','MU-','MU+','HADRONS',
|
|---|
| 1188 | * 'CHARGED','NUCLEI','CERENKOV',
|
|---|
| 1189 | * (J*THSTEP,(PLONG(J,K),K=1,9),J=LPCT1,NSTEP)
|
|---|
| 1190 | 910 FORMAT(/' ---------- LONGITUDINAL DISTRIBUTION IN STEPS OF ',
|
|---|
| 1191 | * F5.0,' G/CM**2 ----------------'/
|
|---|
| 1192 | * ' DEPTH ',3A12,3A11,A12,A11,A12/
|
|---|
| 1193 | * (' ',F6.0,F13.0,2F12.0,3F11.0,F12.0,F11.0,1P,E12.5,0P) )
|
|---|
| 1194 | CJOK ADAPTED FOR HEAT CALCULATION
|
|---|
| 1195 | C910 FORMAT(/
|
|---|
| 1196 | C * ' LONGITUDINAL DISTRIBUTION IN STEPS OF ',F5.0,' G/CM**2'
|
|---|
| 1197 | C * /' ',92('=')/' DEPTH',8A10,A12/1P
|
|---|
| 1198 | C * (' ',0P,F6.0,1P,9E11.4))
|
|---|
| 1199 | CJOK
|
|---|
| 1200 |
|
|---|
| 1201 | IF ( FLGFIT ) THEN
|
|---|
| 1202 | C PERFORM FIT TO THE LONGITUDINAL DISTRIBUTION OF ALL CHARGED PARTICLES
|
|---|
| 1203 | C IF EGS IS SELECTED THIS IS THE DISTRIBUTION WHICH IS TO BE TAKEN
|
|---|
| 1204 | IF ( FEGS ) THEN
|
|---|
| 1205 | DO 930 J=0,NSTEP-LPCT1
|
|---|
| 1206 | DEP(J+1) = (J+LPCT1)*THSTEP
|
|---|
| 1207 | CHAPAR(J+1) = PLONG(J+LPCT1,7)
|
|---|
| 1208 | 930 CONTINUE
|
|---|
| 1209 | NSTP = NSTEP + 1 - LPCT1
|
|---|
| 1210 | WRITE(MONIOU,8229) 'ALL CHARGED PARTICLES'
|
|---|
| 1211 | 8229 FORMAT(/' FIT OF THE CURVE ',
|
|---|
| 1212 | * ' N(T) = P1*((T-P2)/(P3-P2))**((P3-T)/(P4+P5*T+P6*T**2))'/
|
|---|
| 1213 | * ' TO LONGITUDINAL DISTRIBUTION OF ',A35)
|
|---|
| 1214 | C IF NKG IS SELECTED ONLY THE ELECTRON DISTRIBUTION IS AVAILABLE
|
|---|
| 1215 | ELSEIF ( FNKG ) THEN
|
|---|
| 1216 | DEP(1) = 0.D0
|
|---|
| 1217 | CHAPAR(1) = 0.D0
|
|---|
| 1218 | DO 931 J = 1,IALT(1)
|
|---|
| 1219 | DEP(J+1) = TLEV(J)
|
|---|
| 1220 | CHAPAR(J+1) = SL(J)
|
|---|
| 1221 | 931 CONTINUE
|
|---|
| 1222 | NSTP = IALT(1) + 1
|
|---|
| 1223 | WRITE(MONIOU,8229) 'NKG ELECTRONS'
|
|---|
| 1224 | C IF NONE IS SELECTED IT DOES NOT REALLY MAKE SENSE TO FIT
|
|---|
| 1225 | C BUT LET'S TAKE THEN ALL CHARGED WHICH ARE MUONS AND HADRONS
|
|---|
| 1226 | ELSE
|
|---|
| 1227 | DO 932 J=0,NSTEP-LPCT1
|
|---|
| 1228 | DEP(J+1) = (J+LPCT1)*THSTEP
|
|---|
| 1229 | CHAPAR(J+1) = PLONG(J+LPCT1,7)
|
|---|
| 1230 | 932 CONTINUE
|
|---|
| 1231 | NSTP = NSTEP + 1 - LPCT1
|
|---|
| 1232 | WRITE(MONIOU,8229) 'MUONS AND CHARGED HADRONS'
|
|---|
| 1233 | ENDIF
|
|---|
| 1234 | IF ( NSTP .GT. 6 ) THEN
|
|---|
| 1235 | C THERE ARE MORE THAN 6 STEP VALUES, A FIT SHOULD BE POSSIBLE.
|
|---|
| 1236 | C DO THE FIT: NPAR AND FPARAM GIVE THE NUMBER OF PARAMETERS USED
|
|---|
| 1237 | C AND THE FINAL VALUES FOR THE PARAMETERS. CHISQ GIVES THE CHI**2/DOF
|
|---|
| 1238 | C FOR THE FIT.
|
|---|
| 1239 | CALL LONGFT(FPARAM,CHI2)
|
|---|
| 1240 | WRITE(MONIOU,8230) FPARAM,CHI2,CHI2/SQRT(FPARAM(1))*100.D0
|
|---|
| 1241 | 8230 FORMAT(' PARAMETERS = ',1P,6E12.4,0P/
|
|---|
| 1242 | * ' CHI**2/DOF = ',F10.1/
|
|---|
| 1243 | * ' AV. DEVIATION IN % = ',F10.4)
|
|---|
| 1244 | C STORE RESULT IN END EVENT BLOCK
|
|---|
| 1245 | DO 933 K = 1,6
|
|---|
| 1246 | EVTE(255+K) = FPARAM(K)
|
|---|
| 1247 | 933 CONTINUE
|
|---|
| 1248 | EVTE(262) = CHI2
|
|---|
| 1249 | ELSE
|
|---|
| 1250 | WRITE(MONIOU,*) 'NO LONGI. FIT POSSIBLE, ',
|
|---|
| 1251 | * ' NSTP = ',NSTP,' TOO SMALL.'
|
|---|
| 1252 | DO 934 K = 1,6
|
|---|
| 1253 | EVTE(255+K) = 0.
|
|---|
| 1254 | 934 CONTINUE
|
|---|
| 1255 | EVTE(262) = 0.
|
|---|
| 1256 | ENDIF
|
|---|
| 1257 | ENDIF
|
|---|
| 1258 | ENDIF
|
|---|
| 1259 |
|
|---|
| 1260 | C
|
|---|
| 1261 | CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 1262 | C
|
|---|
| 1263 | C Modified by C. Bigongiari 2001 Jan 16
|
|---|
| 1264 | C
|
|---|
| 1265 | C
|
|---|
| 1266 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1267 | c Saves statistics to sta###### file
|
|---|
| 1268 | C call jcstadata(EVTH,EVTE,
|
|---|
| 1269 | C + NPROTO,NPROTB,NNEUTR,NNEUTB,NPHOTO,NELECT,NPOSIT,
|
|---|
| 1270 | C + NNU ,NMUM ,NMUP ,NPI0 ,NPIM ,NPIP ,NK0L ,
|
|---|
| 1271 | C + NK0S ,NKMI ,NKPL ,NHYP ,NDEUT ,NTRIT ,NALPHA,
|
|---|
| 1272 | C + NOTHER,IFINNU,IFINPI,IFINET,IFINKA,IFINHY,
|
|---|
| 1273 | C + CERELE,CERHAD,PLONG,LPCT1,NSTEP,THSTEP)
|
|---|
| 1274 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1275 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1276 | C WRITE SHOWER END TO OUTPUT BUFFER
|
|---|
| 1277 | Cc CALL TOBUF( EVTE,0 )
|
|---|
| 1278 | C CALL TOBUF( EVTE,1 )
|
|---|
| 1279 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1280 | C IF ( LCERFI ) THEN
|
|---|
| 1281 | C CALL OUTND2
|
|---|
| 1282 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1283 | Cc CALL TOBUFC( EVTE,0 )
|
|---|
| 1284 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1285 | C ENDIF
|
|---|
| 1286 |
|
|---|
| 1287 | C
|
|---|
| 1288 | C WRITE SHOWER END TO OUTPUT BUFFER
|
|---|
| 1289 | CALL TOBUF( EVTE,0 )
|
|---|
| 1290 | IF ( LCERFI ) THEN
|
|---|
| 1291 | CALL OUTND2
|
|---|
| 1292 | CALL TOBUFC( EVTE,0 )
|
|---|
| 1293 | ENDIF
|
|---|
| 1294 |
|
|---|
| 1295 | C
|
|---|
| 1296 | CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 1297 | C
|
|---|
| 1298 |
|
|---|
| 1299 | IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,*)
|
|---|
| 1300 | * 'CERENKOV PH. FROM ELECTRONS = ',SNGL(CERELE),
|
|---|
| 1301 | * ' CERENKOV PH. FROM HADRONS = ',SNGL(CERHAD)
|
|---|
| 1302 | CERELE = 0.D0
|
|---|
| 1303 | CERHAD = 0.D0
|
|---|
| 1304 | NRECER = 0
|
|---|
| 1305 |
|
|---|
| 1306 | IF ( FPRINT .OR. DEBUG ) WRITE(MONIOU,210) SHOWNO
|
|---|
| 1307 | 210 FORMAT(/' END OF SHOWER NO ',I10)
|
|---|
| 1308 |
|
|---|
| 1309 | DO 19 J = 1,37
|
|---|
| 1310 | JNBIN(J) = JNBIN(J) + INBIN(J)
|
|---|
| 1311 | JPBIN(J) = JPBIN(J) + IPBIN(J)
|
|---|
| 1312 | JKBIN(J) = JKBIN(J) + IKBIN(J)
|
|---|
| 1313 | JHBIN(J) = JHBIN(J) + IHBIN(J)
|
|---|
| 1314 | 19 CONTINUE
|
|---|
| 1315 |
|
|---|
| 1316 |
|
|---|
| 1317 | 2 CONTINUE
|
|---|
| 1318 |
|
|---|
| 1319 | C END OF SHOWER LOOP
|
|---|
| 1320 |
|
|---|
| 1321 | C-----------------------------------------------------------------------
|
|---|
| 1322 | 992 CONTINUE
|
|---|
| 1323 |
|
|---|
| 1324 | C RESET NUMBER OF SHOWERS TO CORRECT VALUE
|
|---|
| 1325 | ISHW = I
|
|---|
| 1326 |
|
|---|
| 1327 | RUNE(3) = REAL(ISHW)
|
|---|
| 1328 | C WRITE RUN END TO OUTPUT BUFFER AND FINISH OUTPUT
|
|---|
| 1329 | C
|
|---|
| 1330 | CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 1331 | C
|
|---|
| 1332 | C Modified by C. Bigongiari 2001 Jan 16
|
|---|
| 1333 | C
|
|---|
| 1334 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1335 | Cc CALL TOBUF ( RUNE,1 )
|
|---|
| 1336 | C call jcendrun(rune)
|
|---|
| 1337 | Cc IF ( LCERFI ) CALL TOBUFC( RUNE,1 )
|
|---|
| 1338 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1339 | C
|
|---|
| 1340 | C- write Run End
|
|---|
| 1341 | C
|
|---|
| 1342 |
|
|---|
| 1343 | CALL TOBUF ( RUNE,1 )
|
|---|
| 1344 | IF ( LCERFI ) CALL TOBUFC( RUNE,1)
|
|---|
| 1345 |
|
|---|
| 1346 | C
|
|---|
| 1347 | CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 1348 | C
|
|---|
| 1349 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1350 | C TIME SINCE BEGINNING
|
|---|
| 1351 | c ILEFTB = TIME()
|
|---|
| 1352 | ILEFTB = 1
|
|---|
| 1353 | c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1354 | TDIFF = ILEFTB - ILEFTA
|
|---|
| 1355 |
|
|---|
| 1356 | C MEAN VALUE FOR FIRST INTERACTION ALTITUDE (G/CM**2)
|
|---|
| 1357 | IF ( ISHW .GT. 1 ) THEN
|
|---|
| 1358 | CHISM2 = SQRT( ABS(CHISM2-CHISUM**2/ISHW) / (ISHW-1) )
|
|---|
| 1359 | CHISUM = CHISUM / ISHW
|
|---|
| 1360 | ELSE
|
|---|
| 1361 | CHISM2 = 0.D0
|
|---|
| 1362 | ENDIF
|
|---|
| 1363 |
|
|---|
| 1364 | C OUTPUTS FOR ALL SHOWERS
|
|---|
| 1365 | WRITE(MONIOU,201) ISHW,TDIFF,TDIFF/ISHW,IRECOR,IRECOR/ISHW,
|
|---|
| 1366 | * CHISUM,CHISM2
|
|---|
| 1367 | 201 FORMAT('1',10('='),' RUN SUMMARY ',56('=')//
|
|---|
| 1368 | * ' NUMBER OF GENERATED EVENTS = ',I10,/
|
|---|
| 1369 | * ' TOTAL TIME USED = ',E10.3,' SEC'/
|
|---|
| 1370 | * ' TIME PER EVENT = ',E10.3,' SEC'/
|
|---|
| 1371 | * ' TOTAL SPACE ON PATAPE USED = ',I10,' WORDS'/
|
|---|
| 1372 | * ' SPACE PER EVENT ON PATAPE = ',I10,' WORDS'/
|
|---|
| 1373 | * ' AVERAGE HEIGHT OF 1ST INT. = ',F10.3,' +-',F10.3,' G/CM**2'/)
|
|---|
| 1374 |
|
|---|
| 1375 | C ENERGY - MULTIPLICITY MATRIX FOR ALL SHOWERS
|
|---|
| 1376 | WRITE(MONIOU,209) (K,K=1,13),
|
|---|
| 1377 | * (J,(MULTOT(J,K),K=1,13),10**((J-4.)/3.),10**((J-3.)/3.),J=1,37),
|
|---|
| 1378 | * 1,(INT(10**((K-1.)/3.)+1),K = 2,13),
|
|---|
| 1379 | * 2,(INT(10**((K )/3.) ),K = 2,13)
|
|---|
| 1380 | 209 FORMAT(//' ENERGY - MULTIPLICITY MATRIX FOR ALL SHOWERS'/
|
|---|
| 1381 | * ' ENERGY RUNS VERTICALLY, MULTIPLICITY HORIZONTALLY'//,
|
|---|
| 1382 | * ' ',6X,5I9,3I8,5I7' ENERGY RANGE (GEV)'/
|
|---|
| 1383 | * 37(/' ',I4,1X,I10,4I9,3I8,5I7,1X,1P,2E10.1,0P)//
|
|---|
| 1384 | * ' MULT. ',5I9,3I8,5I7,4X,'LOWER BIN LIMIT'/
|
|---|
| 1385 | * ' RANGE ',5I9,3I8,5I7,4X,'UPPER BIN LIMIT')
|
|---|
| 1386 |
|
|---|
| 1387 | C GET MEAN OF ELASTICITY FOR ENERGY BINS
|
|---|
| 1388 | DO 3377 J = 1,37
|
|---|
| 1389 | NELMEA = 0
|
|---|
| 1390 | DO 3378 K = 1,10
|
|---|
| 1391 | NELMEA = NELMEA + IELDPA(J,K)
|
|---|
| 1392 | 3378 CONTINUE
|
|---|
| 1393 | IF ( NELMEA .NE. 0 ) ELMEAA(J) = ELMEAA(J) / NELMEA
|
|---|
| 1394 | 3377 CONTINUE
|
|---|
| 1395 |
|
|---|
| 1396 | C PRINT ENERGY - ELASTICITY MATRIX FOR ALL SHOWERS
|
|---|
| 1397 | WRITE(MONIOU,408) (K,K=1,10), (J,(IELDPA(J,K),K=1,10),
|
|---|
| 1398 | * ELMEAA(J),10**((J-4.D0)/3.D0),10**((J-3.)/3.D0),J=1,37),
|
|---|
| 1399 | * ((K-1)*0.1D0,K=1,10),(K*0.1D0,K=1,10)
|
|---|
| 1400 | 408 FORMAT (//' ENERGY - ELASTICITY MATRIX FOR ALL SHOWERS'/
|
|---|
| 1401 | * ' ENERGY RUNS VERTICALLY, ELASTICITY HORIZONTALLY'//
|
|---|
| 1402 | * ' ',5X,10I9,' MEAN EL. ENERGY RANGE (GEV)'/
|
|---|
| 1403 | * 37(/' ',I4,1X,10I9,2X,1P,E10.3,2E10.1,0P)//
|
|---|
| 1404 | * ' ELA. ',10F9.2,5X,'LOWER BIN LIMIT'/
|
|---|
| 1405 | * ' RANGE',10F9.2,5X,'UPPER BIN LIMIT')
|
|---|
| 1406 |
|
|---|
| 1407 | WRITE(MONIOU,204)
|
|---|
| 1408 | 204 FORMAT (//' INTERACTIONS PER KINETIC ENERGY INTERVAL FOR ALL ',
|
|---|
| 1409 | * 'SHOWERS'//' BIN LOWER LIMIT UPPER LIMIT ',
|
|---|
| 1410 | * 'NUCLEON PIONS KAONS S.BARYONS TOTAL'/
|
|---|
| 1411 | * 12X,'IN GEV',9X,'IN GEV',7X,
|
|---|
| 1412 | * ' EVENTS EVENTS EVENTS EVENTS '//)
|
|---|
| 1413 | WRITE(MONIOU,207) (I,SABIN(I),SBBIN(I),JNBIN(I),JPBIN(I),JKBIN(I)
|
|---|
| 1414 | * ,JHBIN(I),JNBIN(I)+JPBIN(I)+JKBIN(I)+JHBIN(I),I=1,37)
|
|---|
| 1415 | 207 FORMAT(' ',I5,1P,2E15.4,0P,I12,3I10,I11)
|
|---|
| 1416 |
|
|---|
| 1417 | IF ( .NOT.GHEISH ) THEN
|
|---|
| 1418 | C PRINT ELASTICITY STATISTICS
|
|---|
| 1419 | WRITE(MONIOU,89) (I,(I-1)*.05,I*.05,
|
|---|
| 1420 | * IELIS(I),IELHM(I),IELNU(I),IELPI(I),I = 1,20)
|
|---|
| 1421 | 89 FORMAT (//' ELASTICITY STATISTICS '//
|
|---|
| 1422 | * ' BIN LOW HIGH EDGE FOR ISOBARS HEAVY MESONS',
|
|---|
| 1423 | * ' SINGLE NUCLEONS AND PIONS'/
|
|---|
| 1424 | * (' ',I3,' ',F4.2,' ',F4.2,' ',4I17))
|
|---|
| 1425 | ENDIF
|
|---|
| 1426 |
|
|---|
| 1427 | C CALCULATE MEAN VALUES AND STANDARD DEVIATIONS OF PARTICLE NUMBERS
|
|---|
| 1428 | IF ( ISHW .GT. 1 ) THEN
|
|---|
| 1429 | DO 879 K = 1,25
|
|---|
| 1430 | DO 879 J = 1,NOBSLV
|
|---|
| 1431 | MPART2(J,K) = SQRT( abs(MPART2(J,K)-MPARTO(J,K)**2/ISHW)
|
|---|
| 1432 | * /(ISHW-1) )
|
|---|
| 1433 | MPARTO(J,K) = MPARTO(J,K)/ISHW
|
|---|
| 1434 | 879 CONTINUE
|
|---|
| 1435 | ELSE
|
|---|
| 1436 | DO 880 K = 1,25
|
|---|
| 1437 | DO 880 J = 1,NOBSLV
|
|---|
| 1438 | MPART2(J,K) = 0.D0
|
|---|
| 1439 | 880 CONTINUE
|
|---|
| 1440 | ENDIF
|
|---|
| 1441 |
|
|---|
| 1442 | C PRINT MEAN VALUES AND STANDARD DEVIATIONS OF PARTICLE NUMBERS
|
|---|
| 1443 | IOBSLV = MIN( 3, NOBSLV )
|
|---|
| 1444 | WRITE(MONIOU,854) (K,K=1,IOBSLV)
|
|---|
| 1445 | 854 FORMAT (/ ' AVERAGE NUMBER OF PARTICLES PER EVENT :'/
|
|---|
| 1446 | * ' FROM LEVEL NUMBER ', 3(10X,I10,10X) )
|
|---|
| 1447 | WRITE(MONIOU,855) (OBSLEV(K),K=1,IOBSLV)
|
|---|
| 1448 | 855 FORMAT ( ' HEIGHT IN CM',1P,3(20X,E10.3)/)
|
|---|
| 1449 |
|
|---|
| 1450 | WRITE(MONIOU,778)'PROTONS ',(MPROTO(K),MPROT2(K),K=1,IOBSLV)
|
|---|
| 1451 | WRITE(MONIOU,778)'ANTIPROTONS ',(MPROTB(K),MPRTB2(K),K=1,IOBSLV)
|
|---|
| 1452 | WRITE(MONIOU,778)'NEUTRONS ',(MNEUTR(K),MNETR2(K),K=1,IOBSLV)
|
|---|
| 1453 | WRITE(MONIOU,778)'ANTINEUTRONS',(MNEUTB(K),MNETB2(K),K=1,IOBSLV)
|
|---|
| 1454 | WRITE(MONIOU,778)'PHOTONS ',(MPHOTO(K),MPHOT2(K),K=1,IOBSLV)
|
|---|
| 1455 | WRITE(MONIOU,778)'ELECTRONS ',(MELECT(K),MELEC2(K),K=1,IOBSLV)
|
|---|
| 1456 | WRITE(MONIOU,778)'POSITRONS ',(MPOSIT(K),MPOSI2(K),K=1,IOBSLV)
|
|---|
| 1457 | WRITE(MONIOU,778)'NEUTRINOS ',(MNU (K),MNU2 (K),K=1,IOBSLV)
|
|---|
| 1458 | WRITE(MONIOU,778)'MU - ',(MMUM (K),MMUM2 (K),K=1,IOBSLV)
|
|---|
| 1459 | WRITE(MONIOU,778)'MU + ',(MMUP (K),MMUP2 (K),K=1,IOBSLV)
|
|---|
| 1460 | WRITE(MONIOU,778)'PI 0 ',(MPI0 (K),MPI02 (K),K=1,IOBSLV)
|
|---|
| 1461 | WRITE(MONIOU,778)'PI - ',(MPIM (K),MPIM2 (K),K=1,IOBSLV)
|
|---|
| 1462 | WRITE(MONIOU,778)'PI + ',(MPIP (K),MPIP2 (K),K=1,IOBSLV)
|
|---|
| 1463 | WRITE(MONIOU,778)'K0L ',(MK0L (K),MK0L2 (K),K=1,IOBSLV)
|
|---|
| 1464 | WRITE(MONIOU,778)'K0S ',(MK0S (K),MK0S2 (K),K=1,IOBSLV)
|
|---|
| 1465 | WRITE(MONIOU,778)'K - ',(MKMI (K),MKMI2 (K),K=1,IOBSLV)
|
|---|
| 1466 | WRITE(MONIOU,778)'K + ',(MKPL (K),MKPL2 (K),K=1,IOBSLV)
|
|---|
| 1467 | WRITE(MONIOU,778)'STR. BARYONS',(MHYP (K),MHYP2 (K),K=1,IOBSLV)
|
|---|
| 1468 | WRITE(MONIOU,778)'DEUTERONS ',(MDEUT (K),MDEUT2(K),K=1,IOBSLV)
|
|---|
| 1469 | WRITE(MONIOU,778)'TRITONS ',(MTRIT (K),MTRIT2(K),K=1,IOBSLV)
|
|---|
| 1470 | WRITE(MONIOU,778)'ALPHAS ',(MALPHA(K),MALPH2(K),K=1,IOBSLV)
|
|---|
| 1471 | WRITE(MONIOU,778)'OTHER PART. ',(MOTHER(K),MOTH2 (K),K=1,IOBSLV)
|
|---|
| 1472 | WRITE(MONIOU,*)
|
|---|
| 1473 | 778 FORMAT(' NO OF ',A12,' = ',3(F13.1,' +-',F13.1,' '))
|
|---|
| 1474 |
|
|---|
| 1475 | IF ( NOBSLV .GT. 3 ) THEN
|
|---|
| 1476 | IOBSLV = MIN( 6, NOBSLV )
|
|---|
| 1477 | WRITE(MONIOU,854) (K,K=4,IOBSLV)
|
|---|
| 1478 | WRITE(MONIOU,855) (OBSLEV(K),K=4,IOBSLV)
|
|---|
| 1479 | WRITE(MONIOU,778)'PROTONS ',(MPROTO(K),MPROT2(K),K=4,IOBSLV)
|
|---|
| 1480 | WRITE(MONIOU,778)'ANTIPROTONS ',(MPROTB(K),MPRTB2(K),K=4,IOBSLV)
|
|---|
| 1481 | WRITE(MONIOU,778)'NEUTRONS ',(MNEUTR(K),MNETR2(K),K=4,IOBSLV)
|
|---|
| 1482 | WRITE(MONIOU,778)'ANTINEUTRONS',(MNEUTB(K),MNETB2(K),K=4,IOBSLV)
|
|---|
| 1483 | WRITE(MONIOU,778)'PHOTONS ',(MPHOTO(K),MPHOT2(K),K=4,IOBSLV)
|
|---|
| 1484 | WRITE(MONIOU,778)'ELECTRONS ',(MELECT(K),MELEC2(K),K=4,IOBSLV)
|
|---|
| 1485 | WRITE(MONIOU,778)'POSITRONS ',(MPOSIT(K),MPOSI2(K),K=4,IOBSLV)
|
|---|
| 1486 | WRITE(MONIOU,778)'NEUTRINOS ',(MNU (K),MNU2 (K),K=4,IOBSLV)
|
|---|
| 1487 | WRITE(MONIOU,778)'MU - ',(MMUM (K),MMUM2 (K),K=4,IOBSLV)
|
|---|
| 1488 | WRITE(MONIOU,778)'MU + ',(MMUP (K),MMUP2 (K),K=4,IOBSLV)
|
|---|
| 1489 | WRITE(MONIOU,778)'PI 0 ',(MPI0 (K),MPI02 (K),K=4,IOBSLV)
|
|---|
| 1490 | WRITE(MONIOU,778)'PI - ',(MPIM (K),MPIM2 (K),K=4,IOBSLV)
|
|---|
| 1491 | WRITE(MONIOU,778)'PI + ',(MPIP (K),MPIP2 (K),K=4,IOBSLV)
|
|---|
| 1492 | WRITE(MONIOU,778)'K0L ',(MK0L (K),MK0L2 (K),K=4,IOBSLV)
|
|---|
| 1493 | WRITE(MONIOU,778)'K0S ',(MK0S (K),MK0S2 (K),K=4,IOBSLV)
|
|---|
| 1494 | WRITE(MONIOU,778)'K - ',(MKMI (K),MKMI2 (K),K=4,IOBSLV)
|
|---|
| 1495 | WRITE(MONIOU,778)'K + ',(MKPL (K),MKPL2 (K),K=4,IOBSLV)
|
|---|
| 1496 | WRITE(MONIOU,778)'STR. BARYONS',(MHYP (K),MHYP2 (K),K=4,IOBSLV)
|
|---|
| 1497 | WRITE(MONIOU,778)'DEUTERONS ',(MDEUT (K),MDEUT2(K),K=4,IOBSLV)
|
|---|
| 1498 | WRITE(MONIOU,778)'TRITONS ',(MTRIT (K),MTRIT2(K),K=4,IOBSLV)
|
|---|
| 1499 | WRITE(MONIOU,778)'ALPHAS ',(MALPHA(K),MALPH2(K),K=4,IOBSLV)
|
|---|
| 1500 | WRITE(MONIOU,778)'OTHER PART. ',(MOTHER(K),MOTH2 (K),K=4,IOBSLV)
|
|---|
| 1501 | WRITE(MONIOU,*)
|
|---|
| 1502 |
|
|---|
| 1503 | IF ( NOBSLV .GT. 6 ) THEN
|
|---|
| 1504 | IOBSLV = MIN( 9, NOBSLV )
|
|---|
| 1505 | WRITE(MONIOU,854) (K,K=7,IOBSLV)
|
|---|
| 1506 | WRITE(MONIOU,855) (OBSLEV(K),K=7,IOBSLV)
|
|---|
| 1507 | WRITE(MONIOU,778)'PROTONS ',(MPROTO(K),MPROT2(K),K=7,IOBSLV)
|
|---|
| 1508 | WRITE(MONIOU,778)'ANTIPROTONS ',(MPROTB(K),MPRTB2(K),K=7,IOBSLV)
|
|---|
| 1509 | WRITE(MONIOU,778)'NEUTRONS ',(MNEUTR(K),MNETR2(K),K=7,IOBSLV)
|
|---|
| 1510 | WRITE(MONIOU,778)'ANTINEUTRONS',(MNEUTB(K),MNETB2(K),K=7,IOBSLV)
|
|---|
| 1511 | WRITE(MONIOU,778)'PHOTONS ',(MPHOTO(K),MPHOT2(K),K=7,IOBSLV)
|
|---|
| 1512 | WRITE(MONIOU,778)'ELECTRONS ',(MELECT(K),MELEC2(K),K=7,IOBSLV)
|
|---|
| 1513 | WRITE(MONIOU,778)'POSITRONS ',(MPOSIT(K),MPOSI2(K),K=7,IOBSLV)
|
|---|
| 1514 | WRITE(MONIOU,778)'NEUTRINOS ',(MNU (K),MNU2 (K),K=7,IOBSLV)
|
|---|
| 1515 | WRITE(MONIOU,778)'MU - ',(MMUM (K),MMUM2 (K),K=7,IOBSLV)
|
|---|
| 1516 | WRITE(MONIOU,778)'MU + ',(MMUP (K),MMUP2 (K),K=7,IOBSLV)
|
|---|
| 1517 | WRITE(MONIOU,778)'PI 0 ',(MPI0 (K),MPI02 (K),K=7,IOBSLV)
|
|---|
| 1518 | WRITE(MONIOU,778)'PI - ',(MPIM (K),MPIM2 (K),K=7,IOBSLV)
|
|---|
| 1519 | WRITE(MONIOU,778)'PI + ',(MPIP (K),MPIP2 (K),K=7,IOBSLV)
|
|---|
| 1520 | WRITE(MONIOU,778)'K0L ',(MK0L (K),MK0L2 (K),K=7,IOBSLV)
|
|---|
| 1521 | WRITE(MONIOU,778)'K0S ',(MK0S (K),MK0S2 (K),K=7,IOBSLV)
|
|---|
| 1522 | WRITE(MONIOU,778)'K - ',(MKMI (K),MKMI2 (K),K=7,IOBSLV)
|
|---|
| 1523 | WRITE(MONIOU,778)'K + ',(MKPL (K),MKPL2 (K),K=7,IOBSLV)
|
|---|
| 1524 | WRITE(MONIOU,778)'STR. BARYONS',(MHYP (K),MHYP2 (K),K=7,IOBSLV)
|
|---|
| 1525 | WRITE(MONIOU,778)'DEUTERONS ',(MDEUT (K),MDEUT2(K),K=7,IOBSLV)
|
|---|
| 1526 | WRITE(MONIOU,778)'TRITONS ',(MTRIT (K),MTRIT2(K),K=7,IOBSLV)
|
|---|
| 1527 | WRITE(MONIOU,778)'ALPHAS ',(MALPHA(K),MALPH2(K),K=7,IOBSLV)
|
|---|
| 1528 | WRITE(MONIOU,778)'OTHER PART. ',(MOTHER(K),MOTH2 (K),K=7,IOBSLV)
|
|---|
| 1529 | WRITE(MONIOU,*)
|
|---|
| 1530 |
|
|---|
| 1531 | IF ( NOBSLV .GT. 9 ) THEN
|
|---|
| 1532 | IOBSLV = MIN( 10, NOBSLV )
|
|---|
| 1533 | WRITE(MONIOU,854) (K,K=9,IOBSLV)
|
|---|
| 1534 | WRITE(MONIOU,855) (OBSLEV(K),K=9,IOBSLV)
|
|---|
| 1535 | WRITE(MONIOU,778)'PROTONS ',(MPROTO(K),MPROT2(K),K=9,IOBSLV)
|
|---|
| 1536 | WRITE(MONIOU,778)'ANTIPROTONS ',(MPROTB(K),MPRTB2(K),K=9,IOBSLV)
|
|---|
| 1537 | WRITE(MONIOU,778)'NEUTRONS ',(MNEUTR(K),MNETR2(K),K=9,IOBSLV)
|
|---|
| 1538 | WRITE(MONIOU,778)'ANTINEUTRONS',(MNEUTB(K),MNETB2(K),K=9,IOBSLV)
|
|---|
| 1539 | WRITE(MONIOU,778)'PHOTONS ',(MPHOTO(K),MPHOT2(K),K=9,IOBSLV)
|
|---|
| 1540 | WRITE(MONIOU,778)'ELECTRONS ',(MELECT(K),MELEC2(K),K=9,IOBSLV)
|
|---|
| 1541 | WRITE(MONIOU,778)'POSITRONS ',(MPOSIT(K),MPOSI2(K),K=9,IOBSLV)
|
|---|
| 1542 | WRITE(MONIOU,778)'NEUTRINOS ',(MNU (K),MNU2 (K),K=9,IOBSLV)
|
|---|
| 1543 | WRITE(MONIOU,778)'MU - ',(MMUM (K),MMUM2 (K),K=9,IOBSLV)
|
|---|
| 1544 | WRITE(MONIOU,778)'MU + ',(MMUP (K),MMUP2 (K),K=9,IOBSLV)
|
|---|
| 1545 | WRITE(MONIOU,778)'PI 0 ',(MPI0 (K),MPI02 (K),K=9,IOBSLV)
|
|---|
| 1546 | WRITE(MONIOU,778)'PI - ',(MPIM (K),MPIM2 (K),K=9,IOBSLV)
|
|---|
| 1547 | WRITE(MONIOU,778)'PI + ',(MPIP (K),MPIP2 (K),K=9,IOBSLV)
|
|---|
| 1548 | WRITE(MONIOU,778)'K0L ',(MK0L (K),MK0L2 (K),K=9,IOBSLV)
|
|---|
| 1549 | WRITE(MONIOU,778)'K0S ',(MK0S (K),MK0S2 (K),K=9,IOBSLV)
|
|---|
| 1550 | WRITE(MONIOU,778)'K - ',(MKMI (K),MKMI2 (K),K=9,IOBSLV)
|
|---|
| 1551 | WRITE(MONIOU,778)'K + ',(MKPL (K),MKPL2 (K),K=9,IOBSLV)
|
|---|
| 1552 | WRITE(MONIOU,778)'STR. BARYONS',(MHYP (K),MHYP2 (K),K=9,IOBSLV)
|
|---|
| 1553 | WRITE(MONIOU,778)'DEUTERONS ',(MDEUT (K),MDEUT2(K),K=9,IOBSLV)
|
|---|
| 1554 | WRITE(MONIOU,778)'TRITONS ',(MTRIT (K),MTRIT2(K),K=9,IOBSLV)
|
|---|
| 1555 | WRITE(MONIOU,778)'ALPHAS ',(MALPHA(K),MALPH2(K),K=9,IOBSLV)
|
|---|
| 1556 | WRITE(MONIOU,778)'OTHER PART. ',(MOTHER(K),MOTH2 (K),K=9,IOBSLV)
|
|---|
| 1557 | WRITE(MONIOU,*)
|
|---|
| 1558 | ENDIF
|
|---|
| 1559 |
|
|---|
| 1560 | ENDIF
|
|---|
| 1561 | ENDIF
|
|---|
| 1562 |
|
|---|
| 1563 | C PRINT OUT NKG RESULT FOR ALL SHOWERS IF SELECTED
|
|---|
| 1564 | IF ( FNKG ) CALL MITAGE
|
|---|
| 1565 |
|
|---|
| 1566 | C CALCULATE MEAN VALUES AND SIGMAS OF LONGITUDINAL DISTRIBUTION
|
|---|
| 1567 | IF ( LLONGI ) THEN
|
|---|
| 1568 | IF ( ISHW .GT. 1 ) THEN
|
|---|
| 1569 | DO 790 K = 1,9
|
|---|
| 1570 | DO 789 J = LPCT1,NSTEP
|
|---|
| 1571 | SPLONG(J,K) = SQRT( abs(SPLONG(J,K)-APLONG(J,K)**2/ISHW)
|
|---|
| 1572 | * /(ISHW-1) )
|
|---|
| 1573 | APLONG(J,K) = APLONG(J,K)/ISHW
|
|---|
| 1574 | 789 CONTINUE
|
|---|
| 1575 | 790 CONTINUE
|
|---|
| 1576 | ELSE
|
|---|
| 1577 | DO 990 K = 1,9
|
|---|
| 1578 | DO 989 J = LPCT1,NSTEP
|
|---|
| 1579 | SPLONG(J,K) = 0.D0
|
|---|
| 1580 | 989 CONTINUE
|
|---|
| 1581 | 990 CONTINUE
|
|---|
| 1582 | ENDIF
|
|---|
| 1583 |
|
|---|
| 1584 | C PRINT AVERAGE LONGITUDINAL DISTRIBUTIONS
|
|---|
| 1585 | WRITE(MONIOU,911) THSTEP,
|
|---|
| 1586 | * 'GAMMAS ','POSITRONS','ELECTRONS','MU- ','MU+ ',
|
|---|
| 1587 | * (J*THSTEP,(APLONG(J,K),SPLONG(J,K),K=1,5),J=LPCT1,NSTEP)
|
|---|
| 1588 | 911 FORMAT(/' AVERAGE LONGITUDINAL DISTRIBUTION IN STEPS OF ',
|
|---|
| 1589 | * F5.0,' G/CM**2 '/' ',131('=')/
|
|---|
| 1590 | * ' DEPTH',8X,3(A10,16X),A9,15X,A9 //
|
|---|
| 1591 | * (' ',F5.0,2X,1P,E11.4,'+-',E11.4,0P,1X,F12.0,'+-',F11.0,
|
|---|
| 1592 | * 1X,F12.0,'+-',E11.4,1X,F11.1,'+-',F10.1,
|
|---|
| 1593 | * 1X,F11.1,'+-',F10.1 ))
|
|---|
| 1594 | WRITE(MONIOU,912) THSTEP,
|
|---|
| 1595 | * 'HADRONS','CHARGED','NUCLEI','CERENKOV',
|
|---|
| 1596 | * (J*THSTEP,(APLONG(J,K),SPLONG(J,K),K=6,9),J=LPCT1,NSTEP)
|
|---|
| 1597 | 912 FORMAT(/' AVERAGE LONGITUDINAL DISTRIBUTION IN STEPS OF ',
|
|---|
| 1598 | * F5.0,' G/CM**2 '/' ',115('=')/
|
|---|
| 1599 | * ' DEPTH',8X,A9,16X,A10,16X,A9,21X,A9 //
|
|---|
| 1600 | * (' ',F5.0,1X,F11.1,'+-',F11.1,1X,F12.0,'+-',F12.0,
|
|---|
| 1601 | * 2X,F10.1,'+-',F10.1,1X,1P,E16.6,'+-',E16.6,0P))
|
|---|
| 1602 | ENDIF
|
|---|
| 1603 |
|
|---|
| 1604 | IF ( FLGFIT ) THEN
|
|---|
| 1605 | C PERFORM FIT TO THE LONGITUDINAL DISTRIBUTION OF ALL CHARGED PARTICLES
|
|---|
| 1606 | C IF EGS IS SELECTED THIS IS THE DISTRIBUTION WHICH IS TO BE TAKEN
|
|---|
| 1607 | IF ( FEGS ) THEN
|
|---|
| 1608 | DO 730 J=0,NSTEP-LPCT1
|
|---|
| 1609 | DEP(J+1) = (J+LPCT1)*THSTEP
|
|---|
| 1610 | CHAPAR(J+1) = APLONG(J+LPCT1,7)
|
|---|
| 1611 | 730 CONTINUE
|
|---|
| 1612 | NSTP = NSTEP + 1 - LPCT1
|
|---|
| 1613 | WRITE(MONIOU,8229) 'AVERAGE ALL CHARGED PARTICLES'
|
|---|
| 1614 | C IF NKG IS SELECTED ONLY THE ELECTRON DISTRIBUTION IS AVAILABLE
|
|---|
| 1615 | ELSEIF ( FNKG ) THEN
|
|---|
| 1616 | DEP(1) = 0.D0
|
|---|
| 1617 | CHAPAR(1) = 0.D0
|
|---|
| 1618 | DO 731 J = 1,IALT(1)
|
|---|
| 1619 | DEP(J+1) = TLEV(J)
|
|---|
| 1620 | CHAPAR(J+1) = SEL(J)/ISHW
|
|---|
| 1621 | 731 CONTINUE
|
|---|
| 1622 | NSTP = IALT(1) + 1
|
|---|
| 1623 | WRITE(MONIOU,8229) 'AVERAGE NKG ELECTRONS'
|
|---|
| 1624 | C IF NONE IS SELECTED IT DOES NOT REALLY MAKE SENSE TO FIT
|
|---|
| 1625 | C BUT LET'S TAKE THEN ALL CHARGED WHICH ARE MUONS AND HADRONS
|
|---|
| 1626 | ELSE
|
|---|
| 1627 | DO 732 J=0,NSTEP-LPCT1
|
|---|
| 1628 | DEP(J+1) = (J+LPCT1)*THSTEP
|
|---|
| 1629 | CHAPAR(J+1) = APLONG(J+LPCT1,7)
|
|---|
| 1630 | 732 CONTINUE
|
|---|
| 1631 | NSTP = NSTEP + 1 - LPCT1
|
|---|
| 1632 | WRITE(MONIOU,8229) 'AVERAGE MUONS AND CHARGED HADRONS'
|
|---|
| 1633 | ENDIF
|
|---|
| 1634 | IF ( NSTP .GT. 6 ) THEN
|
|---|
| 1635 | C THERE ARE MORE THAN 6 STEP VALUES, A FIT SHOULD BE POSSIBLE.
|
|---|
| 1636 | C DO THE FIT: NPAR AND FPARAM GIVE THE NUMBER OF PARAMETERS USED
|
|---|
| 1637 | C AND THE FINAL VALUES FOR THE PARAMETERS. CHISQ GIVES THE CHI**2/DOF
|
|---|
| 1638 | C FOR THE FIT.
|
|---|
| 1639 | CALL LONGFT(FPARAM,CHI2)
|
|---|
| 1640 | WRITE(MONIOU,8230) FPARAM,CHI2,CHI2/SQRT(FPARAM(1))*100.D0
|
|---|
| 1641 | ELSE
|
|---|
| 1642 | WRITE(MONIOU,*) 'NO LONGI. FIT POSSIBLE, ',
|
|---|
| 1643 | * ' NSTP = ',NSTP,' TOO SMALL.'
|
|---|
| 1644 | ENDIF
|
|---|
| 1645 | ENDIF
|
|---|
| 1646 |
|
|---|
| 1647 |
|
|---|
| 1648 | C CONTROL PRINT OUTPUT OF CONSTANTS
|
|---|
| 1649 | IF ( DEBUG ) THEN
|
|---|
| 1650 | CALL STAEND
|
|---|
| 1651 | WRITE(MDEBUG,*) 'MAIN : STAEND CALLED'
|
|---|
| 1652 | ENDIF
|
|---|
| 1653 |
|
|---|
| 1654 | C
|
|---|
| 1655 | CBC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 1656 | C
|
|---|
| 1657 | C Modified by C. Bigongiari 2001 Jan 16
|
|---|
| 1658 | C
|
|---|
| 1659 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1660 | C call jcenddata(runh,rune)
|
|---|
| 1661 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1662 | C
|
|---|
| 1663 | C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 1664 |
|
|---|
| 1665 | WRITE(MONIOU,*)' '
|
|---|
| 1666 | CALL PRTIME(TTIME)
|
|---|
| 1667 | WRITE(MONIOU,101)
|
|---|
| 1668 | 101 FORMAT (/' ',10('='),' END OF RUN ',67('='))
|
|---|
| 1669 |
|
|---|
| 1670 | C CLOSE ALL OPEN UNITS
|
|---|
| 1671 | IF ( MONIOU .NE. 6 ) CLOSE( MONIOU )
|
|---|
| 1672 | IF ( MDEBUG .NE. 6 ) CLOSE( MDEBUG )
|
|---|
| 1673 | CLOSE( EXST )
|
|---|
| 1674 |
|
|---|
| 1675 | C
|
|---|
| 1676 | CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 1677 | C
|
|---|
| 1678 | C Modified by C. Bigongiari 2001 Jan 16
|
|---|
| 1679 | C
|
|---|
| 1680 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1681 | Cc CLOSE( PATAPE )
|
|---|
| 1682 | Cc IF ( LCERFI ) CLOSE( CETAPE )
|
|---|
| 1683 | Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 1684 | C
|
|---|
| 1685 | C- close files
|
|---|
| 1686 | C
|
|---|
| 1687 |
|
|---|
| 1688 | CLOSE( PATAPE )
|
|---|
| 1689 | IF ( LCERFI ) CLOSE( CETAPE )
|
|---|
| 1690 | C CALL JCENDRUN( RUNE )
|
|---|
| 1691 |
|
|---|
| 1692 | C
|
|---|
| 1693 | CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|---|
| 1694 | C
|
|---|
| 1695 |
|
|---|
| 1696 | STOP
|
|---|
| 1697 | END
|
|---|
