source: trunk/MagicSoft/Simulation/Corsika/Mmcs/inprm.f@ 18477

Last change on this file since 18477 was 286, checked in by harald, 25 years ago
This is the start point for further developments of the Magic Monte Carlo Simulation written by Jose Carlos Gonzales. Now it is under control of one CVS repository for the whole collaboration. Everyone should use this CVS repository for further developments.
File size: 30.2 KB
Line 
1C=======================================================================
2
3 SUBROUTINE INPRM
4
5C-----------------------------------------------------------------------
6C IN(PUT) PR(I)M(ARY)
7C
8C TAKES INPUT PRIMARY ENERGY FROM SPECIFIED SPECTRUM
9C CHECKS INPUT VARIABLES FOR CONSISTENCY AND LIMITATIONS
10C WRITES DATA BASE FILE
11C THIS SUBROUTINE IS CALLED FROM MAIN
12C-----------------------------------------------------------------------
13
14 IMPLICIT NONE
15*KEEP,ATMOS.
16 COMMON /ATMOS/ AATM,BATM,CATM,DATM
17 DOUBLE PRECISION AATM(5),BATM(5),CATM(5),DATM(5)
18*KEEP,BUFFS.
19 COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
20 INTEGER MAXBUF,MAXLEN
21 PARAMETER (MAXBUF=39*7)
22 PARAMETER (MAXLEN=12)
23 REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
24 * RUNE(MAXBUF),DATAB(MAXBUF)
25 INTEGER LH
26 CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
27 EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
28 EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
29*KEEP,CONST.
30 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
31 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
32*KEEP,DPMFLG.
33 COMMON /DPMFLG/ NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM
34 INTEGER NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM
35*KEEP,ELABCT.
36 COMMON /ELABCT/ ELCUT
37 DOUBLE PRECISION ELCUT(4)
38*KEEP,ETHMAP.
39 COMMON /ETHMAP/ ECTMAP,ELEFT
40 DOUBLE PRECISION ECTMAP,ELEFT
41*KEEP,LONGI.
42 COMMON /LONGI/ APLONG,HLONG,PLONG,SPLONG,THSTEP,THSTPI,
43 * NSTEP,LLONGI,FLGFIT
44 DOUBLE PRECISION APLONG(0:1040,9),HLONG(0:1024),PLONG(0:1040,9),
45 * SPLONG(0:1040,9),THSTEP,THSTPI
46 INTEGER NSTEP
47 LOGICAL LLONGI,FLGFIT
48*KEEP,MAGANG.
49 COMMON /MAGANG/ ARRANG,ARRANR,COSANG,SINANG
50 DOUBLE PRECISION ARRANG,ARRANR,COSANG,SINANG
51*KEEP,MAGNET.
52 COMMON /MAGNET/ BX,BZ,BVAL,BNORMC,BNORM,COSB,SINB,BLIMIT
53 DOUBLE PRECISION BX,BZ,BVAL,BNORMC
54 REAL BNORM,COSB,SINB,BLIMIT
55*KEEP,NKGI.
56 COMMON /NKGI/ SEL,SELLG,STH,ZEL,ZELLG,ZSL,DIST,
57 * DISX,DISY,DISXY,DISYX,DLAX,DLAY,DLAXY,DLAYX,
58 * OBSATI,RADNKG,RMOL,TLEV,TLEVCM,IALT
59 DOUBLE PRECISION SEL(10),SELLG(10),STH(10),ZEL(10),ZELLG(10),
60 * ZSL(10),DIST(10),
61 * DISX(-10:10),DISY(-10:10),
62 * DISXY(-10:10,2),DISYX(-10:10,2),
63 * DLAX (-10:10,2),DLAY (-10:10,2),
64 * DLAXY(-10:10,2),DLAYX(-10:10,2),
65 * OBSATI(2),RADNKG,RMOL(2),TLEV(10),TLEVCM(10)
66 INTEGER IALT(2)
67*KEEP,OBSPAR.
68 COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
69 * THETPR,PHIPR,NOBSLV
70 DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
71 * THETAP,THETPR(2),PHIP,PHIPR(2)
72 INTEGER NOBSLV
73*KEEP,PARPAR.
74 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
75 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
76 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
77 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
78 INTEGER ITYPE,LEVL
79*KEEP,PARPAE.
80 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
81 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
82 * (CURPAR(4), PHI ), (CURPAR(5), H ),
83 * (CURPAR(6), T ), (CURPAR(7), X ),
84 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
85 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
86 * (CURPAR(12),ECM )
87*KEEP,PRIMSP.
88 COMMON /PRIMSP/ PSLOPE,LLIMIT,ULIMIT,LL,UL,SLEX,ISPEC
89 DOUBLE PRECISION PSLOPE,LLIMIT,ULIMIT,LL,UL,SLEX
90 INTEGER ISPEC
91*KEEP,RANDPA.
92 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
93 DOUBLE PRECISION FAC,U1,U2
94 REAL RD(3000)
95 INTEGER ISEED(103,10),NSEQ
96 LOGICAL KNOR
97*KEEP,REJECT.
98 COMMON /REJECT/ AVNREJ,
99 * ALTMIN,ANEXP,THICKA,THICKD,CUTLN,EONCUT,
100 * FNPRIM
101 DOUBLE PRECISION AVNREJ(10)
102 REAL ALTMIN(10),ANEXP(10),THICKA(10),THICKD(10),
103 * CUTLN,EONCUT
104 LOGICAL FNPRIM
105*KEEP,RUNPAR.
106 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
107 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
108 * MONIOU,MDEBUG,NUCNUC,
109 * CETAPE,
110 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
111 * N1STTR,MDBASE,
112 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
113 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
114 * ,GHEISH,GHESIG
115 COMMON /RUNPAC/ DSN,HOST,USER
116 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
117 REAL STEPFC
118 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
119 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
120 * N1STTR,MDBASE
121 INTEGER CETAPE
122 CHARACTER*79 DSN
123 CHARACTER*20 HOST,USER
124
125 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
126 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
127 * ,GHEISH,GHESIG
128*KEEP,STACKF.
129 COMMON /STACKF/ STACK,STACKP,EXST,NSHIFT,NOUREC,ICOUNT,NTO,NFROM
130 INTEGER MAXSTK
131 PARAMETER (MAXSTK = 12*340*2)
132 DOUBLE PRECISION STACK(MAXSTK)
133 INTEGER STACKP,EXST,NSHIFT,NOUREC,ICOUNT,NTO,NFROM
134*KEEP,VERS.
135 COMMON /VERS/ VERNUM,MVDATE,VERDAT
136 DOUBLE PRECISION VERNUM
137 INTEGER MVDATE
138 CHARACTER*18 VERDAT
139*KEEP,VENUS.
140 COMMON /VENUS/ ISH0,IVERVN,MTAR99,FVENUS,FVENSG
141 INTEGER ISH0,IVERVN,MTAR99
142 LOGICAL FVENUS,FVENSG
143*KEEP,CEREN1.
144 COMMON /CEREN1/ CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD,
145 * CERSIZ,LCERFI
146 DOUBLE PRECISION CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD
147 REAL CERSIZ
148 LOGICAL LCERFI
149*KEEP,CEREN2.
150 COMMON /CEREN2/ PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
151 * DCERX,DCERY,ACERX,ACERY,
152 * XCMAX,YCMAX,EPSX,EPSY,
153 * DCERXI,DCERYI,FCERX,FCERY,
154 * XSCATT,YSCATT,CERXOS,CERYOS,
155 * NCERX,NCERY,ICERML
156 REAL PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
157 * DCERX,DCERY,ACERX,ACERY,
158 * XCMAX,YCMAX,EPSX,EPSY,
159 * DCERXI,DCERYI,FCERX,FCERY,
160 * XSCATT,YSCATT,CERXOS(20),CERYOS(20)
161 INTEGER NCERX,NCERY,ICERML
162*KEND.
163
164 DOUBLE PRECISION EFRAC,HEIGH,H0,OOO,THICK
165 REAL VERVEN
166 INTEGER I,IBL,IDPM,ILONG,ISO,J,L
167 LOGICAL LTHIN
168 EXTERNAL HEIGH,THICK
169 CHARACTER*1 MARK
170 CHARACTER*9 LSTDSN
171C-----------------------------------------------------------------------
172
173 WRITE(MONIOU,504)
174 504 FORMAT(//' ',10('='),' SHOWER PARAMETERS ', 50('=') )
175
176C WRITE ENERGY SPECTRUM TO HEADER
177 RUNH(16) = PSLOPE
178 RUNH(17) = LLIMIT
179 RUNH(18) = ULIMIT
180
181 EVTH(58) = PSLOPE
182 EVTH(59) = LLIMIT
183 EVTH(60) = ULIMIT
184
185 IF ( PRMPAR(1) .GE. 6000.D0 .OR. PRMPAR(1) .LE. 0.D0 ) THEN
186 WRITE(MONIOU,*)'INCORRECT SELECTION OF PRIMARY PARTICLE TYPE = '
187 * ,INT(PRMPAR(1))
188 WRITE(MONIOU,*)'PLEASE READ THE MANUALS'
189 STOP
190 ENDIF
191C CHECK WETHER NUCLEUS IS A SINGLE NUCLEON
192 IF (PRMPAR(1) .EQ. 100.D0 ) PRMPAR(1) = 13.D0
193 IF (PRMPAR(1) .EQ. 101.D0 ) PRMPAR(1) = 14.D0
194 WRITE(MONIOU,*)'PRIMARY PARTICLE IDENTIFICATION IS ',
195 * NINT(PRMPAR(1))
196C CHECK RECOMMENDED ENERGY RANGE
197 IF ( FVENUS .AND.
198 * ULIMIT.GT.2.D7 .AND. PRMPAR(1).GE.8.D0 ) THEN
199 WRITE(MONIOU,502) ULIMIT
200 502 FORMAT(' INTERACTION MODEL DOUBTFUL FOR THE SELECTED PRIMARY ',
201 * 'ENERGY OF ',E10.3,' GEV'/' PLEASE READ THE MANUALS')
202 STOP
203 ENDIF
204
205
206
207
208c> *** modified by fs (22/09/98) *******************************
209
210
211c IF ( PRMPAR(1) .GT. 101.D0 ) THEN
212c IF ( GHEISH ) THEN
213cC GHEISHA CAN TREAT ONLY DEUTERONS, TRITONS, AND ALPHA PARTICLES
214c IF ( PRMPAR(1) .NE. 201.D0 .AND. PRMPAR(1) .NE. 301.D0
215c * .AND. PRMPAR(1) .NE. 402.D0 ) THEN
216c IF ( LLIMIT .LT. HILOELB * INT(PRMPAR(1)/100.D0) ) THEN
217c WRITE(MONIOU,503) INT(PRMPAR(1)/100.D0),LLIMIT
218c STOP
219c ENDIF
220c ENDIF
221c ELSE
222c IF ( LLIMIT .LT. HILOELB * INT(PRMPAR(1)/100.D0) ) THEN
223c WRITE(MONIOU,503) INT(PRMPAR(1)/100.D0),LLIMIT
224c 503 FORMAT(' NUCLEUS WITH A =',I2,' AND PRIMARY ENERGY =',
225c * 1PE10.3,' GEV TOO LOW FOR HIGH ENERGY INTERACTION MODEL'/
226c * ' AND CANNOT BE TREATED BY LOW ENERGY INTERACTION MODEL'/
227c * ' PLEASE READ THE MANUALS')
228c STOP
229c ENDIF
230c ENDIF
231c ENDIF
232
233
234c> *** end of modification ****************************************
235
236C DEFINE ENERGY RANGE AND ENERGY SPECTRUM OF PRIMARY
237 IF ( LLIMIT .EQ. ULIMIT ) THEN
238 ISPEC = 0
239 WRITE(MONIOU,506) LLIMIT
240 506 FORMAT(' PRIMARY ENERGY IS FIXED AT ',1PE10.3,
241 * ' GEV' )
242 ELSE
243 ISPEC = 1
244 WRITE(MONIOU,505) PSLOPE,LLIMIT,ULIMIT
245 505 FORMAT(' PRIMARY ENERGY IS TAKEN FROM SPECTRUM VIA MONTE CARLO'/
246 * 5X,' SLOPE OF PRIMARY SPECTRUM = ',1P,E10.3/
247 * 5X,' LOWER LIMIT CUT-OFF FOR PRIMARY SPECTRUM = ',E10.3,' GEV'/
248 * 5X,' UPPER LIMIT CUT-OFF FOR PRIMARY SPECTRUM = ',E10.3,' GEV'/)
249 IF ( PSLOPE .NE. -1.D0 ) THEN
250 LL = LLIMIT ** (PSLOPE + 1.D0)
251 UL = ULIMIT ** (PSLOPE + 1.D0)
252 SLEX = 1.D0 / (PSLOPE + 1.D0)
253 ELSE
254 LL = ULIMIT / LLIMIT
255 ENDIF
256 ENDIF
257
258C FIRST INTERACTION TARGET FIXED ?
259 IF ( N1STTR .EQ. 1 ) THEN
260 WRITE(MONIOU,508) 'NITROGEN'
261 508 FORMAT(' TARGET OF FIRST INTERACTION IS FIXED TO ',A8)
262 ELSEIF ( N1STTR .EQ. 2 ) THEN
263 WRITE(MONIOU,508) 'OXYGEN '
264 ELSEIF ( N1STTR .EQ. 3 ) THEN
265 WRITE(MONIOU,508) 'ARGON '
266 ELSE
267 N1STTR = 0
268 WRITE(MONIOU,*)'TARGET OF FIRST INTERACTION IS CHOSEN RANDOMLY'
269 ENDIF
270
271C CHECK ANGULAR SETTINGS
272 IF ( THETPR(1) .LT. 0.D0 ) THEN
273 WRITE(MONIOU,*)'UNALLOWED CHOICE OF THETPR = ',SNGL(THETPR(1)),
274 * ' DEGREES'
275 WRITE(MONIOU,*)'PLEASE READ THE MANUALS'
276 STOP
277 ENDIF
278c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
279c> IF ( THETPR(2) .GT. 70.D0 ) THEN
280c> WRITE(MONIOU,*)'UNALLOWED CHOICE OF THETPR = ',SNGL(THETPR(2)),
281c> * ' DEGREES'
282c> WRITE(MONIOU,*)'PLEASE READ THE MANUALS'
283c> STOP
284c> ELSEIF ( THETPR(2) .GT. 45.D0 ) THEN
285c> WRITE(MONIOU,*)'UNALLOWED CHOICE OF THETPR = ',SNGL(THETPR(2)),
286c> * ' DEGREES'
287c> WRITE(MONIOU,*)'#########################################'
288c> WRITE(MONIOU,*)'# IN DOUBTFUL CASES CONTACT THE AUTHORS #'
289c> WRITE(MONIOU,*)'#########################################'
290c> STOP
291c> ENDIF
292c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
293C INCIDENCE ANGLE FIXED ?
294 IF ( THETPR(1) .EQ. THETPR(2) .AND. PHIPR(1) .EQ. PHIPR(2) ) THEN
295 FIXINC = .TRUE.
296 WRITE(MONIOU,517) THETPR(1),PHIPR(1)
297 517 FORMAT(' THETA OF INCIDENCE IS FIXED TO ',F10.2,' DEGREES'/
298 * ' PHI OF INCIDENCE IS FIXED TO ',F10.2,' DEGREES')
299 ELSE
300 FIXINC = .FALSE.
301 WRITE(MONIOU,527) THETPR,PHIPR
302 527 FORMAT(' THETA OF INCIDENCE CHOSEN FROM ',F10.2,'...',F10.2,
303 * ' DEGREES'/
304 * ' PHI OF INCIDENCE CHOSEN FROM ',F10.2,'...',F10.2,
305 * ' DEGREES')
306 ENDIF
307 EVTH(81) = THETPR(1)
308 EVTH(82) = THETPR(2)
309 EVTH(83) = PHIPR(1)
310 EVTH(84) = PHIPR(2)
311 THETPR(1) = THETPR(1)*PI/180.D0
312 THETPR(2) = THETPR(2)*PI/180.D0
313 PHIPR(1) = PHIPR(1) *PI/180.D0
314 PHIPR(2) = PHIPR(2) *PI/180.D0
315
316C-----------------------------------------------------------------------
317C PRMPAR, OBSLEV, NOBSLV
318 PRMPAR(2) = 0.D0
319 PRMPAR(6) = 0.D0
320 PRMPAR(7) = 0.D0
321 PRMPAR(8) = 0.D0
322
323C ORDERING OF OBSERVATION LEVELS FROM TOP TO BOTTOM
324 IF ( NOBSLV .GT. 1 ) THEN
325 215 CONTINUE
326 DO 11 I = 2,NOBSLV
327 IF ( OBSLEV(I) .GT. OBSLEV(I-1) ) THEN
328 OOO = OBSLEV(I)
329 OBSLEV(I) = OBSLEV(I-1)
330 OBSLEV(I-1) = OOO
331 GOTO 215
332 ENDIF
333 11 CONTINUE
334 ENDIF
335C CHECK WETHER OBSERVATION LEVELS ARE IN ALLOWED RANGE
336 DO 12 I = 1,NOBSLV
337 IF ( OBSLEV(I) .GE. HEIGH(0.D0) ) THEN
338 WRITE(MONIOU,120)I,OBSLEV(I),HEIGH(0.D0)
339 120 FORMAT(' UNALLOWED CHOICE OF OBSLEV '/' OBSERVATION LEVEL ',
340 * I2,' IS AT ',F12.3,' CM, WHICH IS ABOVE ',
341 * F12.3,' CM'/' PLEASE READ THE MANUALS')
342 STOP
343 ENDIF
344 IF ( OBSLEV(I) .LE. -1.D5 ) THEN
345 WRITE(MONIOU,121)I,OBSLEV(I)
346 121 FORMAT(' UNALLOWED CHOICE OF OBSLEV '/' OBSERVATION LEVEL ',
347 * I2,' IS AT ',F12.3,' CM, WHICH IS BELOW ',
348 * '-1.D5 CM'/' PLEASE READ THE MANUALS')
349 STOP
350 ENDIF
351 THCKOB(I) = THICK(OBSLEV(I))
352 12 CONTINUE
353
354C WRITE OBSERVATION LEVELS TO HEADER (IN CM)
355 RUNH(5) = REAL(NOBSLV)
356 EVTH(47) = REAL(NOBSLV)
357 DO 114 I = 1,NOBSLV
358 RUNH(5+I) = OBSLEV(I)
359 EVTH(47+I) = OBSLEV(I)
360 114 CONTINUE
361
362C FIRST INTERACTION HEIGHT FIXED ?
363 IF ( FIX1I ) THEN
364 IF ( FIXHEI .GE. HEIGH(0.D0) ) THEN
365 WRITE(MONIOU,122)FIXHEI,HEIGH(0.D0)
366 122 FORMAT(' UNALLOWED CHOICE OF FIXHEI '/' FIRST INTERACTION ',
367 * 'IS FIXED AT ',F12.3,' CM, WHICH IS ABOVE ',
368 * F12.3,' CM'/' PLEASE READ THE MANUALS')
369 STOP
370 ENDIF
371 IF ( FIXHEI .LE. OBSLEV(NOBSLV) ) THEN
372 WRITE(MONIOU,123)FIXHEI,OBSLEV(NOBSLV)
373 123 FORMAT(' UNALLOWED CHOICE OF FIXHEI '/' FIRST INTERACTION ',
374 * 'IS FIXED AT ',F12.3,' CM, '/' WHICH IS BELOW ',
375 * 'LOWEST OBSERVATION LEVEL AT ',F12.3,' CM'
376 * /' PLEASE READ THE MANUALS')
377 STOP
378 ENDIF
379 WRITE(MONIOU,507) FIXHEI
380 507 FORMAT(' HEIGHT OF FIRST INTERACTION IS FIXED TO ',1P,E10.2,
381 * ' CM')
382 IF ( N1STTR .GE. 1 .AND. N1STTR .LE. 3 ) THEN
383 IF ( PRMPAR(1) .LE. 3.D0 ) THEN
384 WRITE(MONIOU,516) INT(PRMPAR(1))
385 516 FORMAT(' TARGET OF FIRST INTERACTION CANNOT BE FIXED FOR ',
386 * 'PRIMARY TYPE ',I5/' PLEASE READ THE MANUALS')
387 STOP
388 ELSEIF ( N1STTR .EQ. 1 ) THEN
389 WRITE(MONIOU,*)'TARGET OF FIRST INTERACTION IS NITROGEN'
390 ELSEIF ( N1STTR .EQ. 2 ) THEN
391 WRITE(MONIOU,*)'TARGET OF FIRST INTERACTION IS OXYGEN'
392 ELSEIF ( N1STTR .EQ. 3 ) THEN
393 WRITE(MONIOU,*)'TARGET OF FIRST INTERACTION IS ARGON'
394 ENDIF
395 ELSE
396 WRITE(MONIOU,*)
397 * 'TARGET OF FIRST INTERACTION IS CHOSEN AT RANDOM'
398 ENDIF
399 ELSE
400 FIXHEI = 0.D0
401 WRITE(MONIOU,*) 'HEIGHT OF FIRST INTERACTION IS CHOSEN RANDOMLY'
402 ENDIF
403
404C STARTING ALTITUDE WITHIN ATMOSPHERE?
405 IF ( THICK0 .LT. 0.D0 ) THEN
406 WRITE(MONIOU,130)THICK0
407 130 FORMAT(' UNALLOWED STARTING ALTITUDE WITH NEGATIVE MASS OVERLAY'
408 * ,E12.3/' PLEASE READ THE MANUALS')
409 STOP
410 ENDIF
411 IF ( THICK0 .GE. THCKOB(NOBSLV) ) THEN
412 WRITE(MONIOU,131) THICK0
413 131 FORMAT(' UNALLOWED STARTING ALTITUDE AT ',F12.3,' G/CM**2',
414 * ' WHICH IS BELOW LOWEST OBSERVATION LEVEL'/
415 * ' PLEASE READ THE MANUALS')
416 STOP
417 ENDIF
418 H0 = HEIGH(THICK0)
419 IF ( THICK0 .EQ. 0.D0 ) THEN
420 WRITE(MONIOU,518) H0, THICK0
421 WRITE(MONIOU,*)' WHICH IS AT TOP OF ATMOSPHERE'
422 ELSE
423 WRITE(MONIOU,518) H0, THICK0
424 ENDIF
425 518 FORMAT(' STARTING ALTITUDE AT ',1P,F13.2,' CM (=',
426 * E7.1,' G/CM**2)')
427 WRITE(MONIOU,203) (OBSLEV(I),THCKOB(I),I=1,NOBSLV)
428 203 FORMAT(/' OBSERVATION LEVELS IN CM AND IN G/CM**2 ',
429 * 1P /(5X, 2E20.8 /))
430
431C LONGITUDINAL SHOWER DEVELOPMENT
432 IF ( LLONGI ) THEN
433 THSTEP = NINT(THSTEP)
434 THSTEP = MAX(THSTEP,1.D0)
435 THSTEP = MIN(THSTEP,1040.D0)
436 THSTPI = 1.D0/THSTEP
437 NSTEP = INT(THCKOB(NOBSLV)*THSTPI)
438 IF ( NSTEP .GE. 1040 ) THEN
439 NSTEP = 1040
440 THSTEP = THCKOB(NOBSLV)/NSTEP
441 WRITE(MONIOU,*)'LONGITUDINAL SHOWER SAMPLING MODIFIED'
442 ENDIF
443 WRITE(MONIOU,925) NSTEP+1,THSTEP
444 925 FORMAT(/' LONGITUDINAL SHOWER DEVELOPMENT:'/
445 * ' SHOWER IS SAMPLED IN ',I4,
446 * ' STEPS OF ',F6.1,' G/CM**2')
447C GET HEIGHT VALUES IN CM FOR USE IN EGS
448 DO 478 J = 0,NSTEP
449 HLONG(J) = HEIGH(J*THSTEP)
450 478 CONTINUE
451 IF ( FLGFIT ) THEN
452 WRITE(MONIOU,*)
453 * ' FIT TO CHARGED PARTICLE LONG. DISTRIBUTION ENABLED'
454 ELSE
455 WRITE(MONIOU,*)
456 * ' FIT TO CHARGED PARTICLE LONG. DISTRIBUTION DISABLED'
457 ENDIF
458 WRITE(MONIOU,*)
459 ENDIF
460
461C-----------------------------------------------------------------------
462C CHECK INPUT OF ENERGY CUTS
463 IF ( GHEISH .AND. ELCUT(1) .LT. 0.05D0 ) THEN
464 WRITE(MONIOU,*)'ELCUT(1) SELECTED INCORRECT TO ',ELCUT(1),' GEV'
465 WRITE(MONIOU,*)'PLEASE READ THE MANUALS '
466 STOP
467 ELSEIF ( .NOT.GHEISH .AND. ELCUT(1) .LT. 0.3D0 ) THEN
468 WRITE(MONIOU,*)'ELCUT(1) SELECTED INCORRECT TO ',ELCUT(1),' GEV'
469 WRITE(MONIOU,*)'PLEASE READ THE MANUALS '
470 STOP
471 ENDIF
472 IF ( ELCUT(2) .LT. 0.05D0 ) THEN
473 WRITE(MONIOU,*)'ELCUT(2) SELECTED INCORRECT TO ',ELCUT(2),' GEV'
474 WRITE(MONIOU,*)'PLEASE READ THE MANUALS '
475 STOP
476 ENDIF
477 IF ( ELCUT(3) .LT. 0.003D0 ) THEN
478 WRITE(MONIOU,*)'ELCUT(3) SELECTED INCORRECT TO ',ELCUT(3),' GEV'
479 WRITE(MONIOU,*)'PLEASE READ THE MANUALS '
480 STOP
481 ENDIF
482 IF ( ELCUT(4) .LT. 0.003D0 ) THEN
483 WRITE(MONIOU,*)'ELCUT(4) SELECTED INCORRECT TO ',ELCUT(4),' GEV'
484 WRITE(MONIOU,*)'PLEASE READ THE MANUALS '
485 STOP
486 ENDIF
487 WRITE(MONIOU,703) ECTMAP,ELCUT
488 703 FORMAT (' PARTICLES WITH LORENTZ FACTOR LARGER THAN',1P,E15.4,
489 * ' ARE PRINTED OUT'/' SHOWER PARTICLES ENERGY CUT :'/
490 * ' FOR HADRONS : ',E15.4,' GEV'/
491 * ' FOR MUONS : ',E15.4,' GEV'/
492 * ' FOR ELECTRONS : ',E15.4,' GEV'/
493 * ' FOR PHOTONS : ',E15.4,' GEV'//)
494
495 DO 774 I = 1,4
496 RUNH(20+I) = ELCUT(I)
497 EVTH(60+I) = ELCUT(I)
498 774 CONTINUE
499
500C-----------------------------------------------------------------------
501C PARAMETERS OF EARTH MAGNETIC FIELD OF MIDDLE EUROPE
502C +X DIRECTION IS NORTH, +Y DIRECTION IS EAST, +Z DIRECTION IS DOWN
503 BVAL = SQRT( BX**2 + BZ**2 )
504C BNORM HAS DIMENSIONS OF MEV/CM
505 BNORM = BVAL * C(25) * 1.D-16
506C BNORMC HAS DIMENSIONS OF GEV/CM
507 BNORMC = BNORM * 1.D-3
508 SINB = BZ / BVAL
509 COSB = BX / BVAL
510 WRITE(MONIOU,*)'EARTH MAGNETIC FIELD STRENGTH IS ',SNGL(BVAL),
511 * ' MICROTESLA'
512 WRITE(MONIOU,*)' WITH INCLINATION ANGLE ',
513 * SNGL(ASIN(SINB)*180./PI),' DEGREES'
514 IF ( BVAL .GE. 10000.D0 ) THEN
515 WRITE(MONIOU,*)'YOU WANT TO MAGNETIZE THE GALAXY ?'
516 WRITE(MONIOU,*)'PLEASE READ THE MANUALS !'
517 STOP
518 ENDIF
519C LIMITING FACTOR FOR STEP SIZE OF ELECTRON IN MAGNETIC FIELD
520 BLIMIT = 0.2 / BNORM
521 EVTH(71) = BX
522 EVTH(72) = BZ
523C ANGLE BETWEEN ARRAY X-DIRECTION AND MAGNETIC NORD
524C POSITIV, IF X-DIRECTION OF ARRAY POINTS TO EASTERN DIRECTION
525 ARRANR = ARRANG * PI / 180.D0
526 COSANG = COS(ARRANR)
527 SINANG = SIN(ARRANR)
528 EVTH(93) = ARRANR
529 IF ( ARRANG .NE. 0.D0 ) THEN
530 WRITE(MONIOU,*)
531 WRITE(MONIOU,*)'DETECTOR COORDINATE SYSTEM IS ROTATED AWAY ',
532 * 'FROM NORTH BY ',SNGL(ARRANG),' DEGREES'
533 ENDIF
534
535C-----------------------------------------------------------------------
536C DEFINE CERENKOV ARRAY
537 NCERX = MAX( NCERX, 1 )
538 NCERY = MAX( NCERY, 1 )
539 ACERX = ABS(ACERX)
540 ACERY = ABS(ACERY)
541 DCERX = MAX( ABS(DCERX), 0.001 )
542 DCERY = MAX( ABS(DCERY), 0.001 )
543 XCMAX = (ACERX + (NCERX-1) * DCERX) * 0.5
544 YCMAX = (ACERY + (NCERY-1) * DCERY) * 0.5
545 DCERXI = 1./DCERX
546 EPSX = ACERX * 0.5 * DCERXI
547 DCERYI = 1./DCERY
548 EPSY = ACERY * 0.5 * DCERYI
549 IF ( MOD(NCERX,2) .EQ. 0 ) THEN
550 FCERX = -0.5
551 ELSE
552 FCERX = 0.0
553 ENDIF
554 IF ( MOD(NCERY,2) .EQ. 0 ) THEN
555 FCERY = -0.5
556 ELSE
557 FCERY = 0.0
558 ENDIF
559
560 WRITE(MONIOU,472)
561 * ACERX*.01,ACERY*.01, DCERX*.01,DCERY*.01, NCERX,NCERY
562 472 FORMAT(/' CERENKOV ARRAY:'/
563 * 5X,' CERENKOV STATIONS ARE ',F6.2,' X ',F6.2,' M**2 LARGE'/
564 * 5X,' THE GRID SPACING IS ',F6.2,' AND ',F6.2,' M',/
565 * 5X,' THERE ARE ',I3,' X ',I3,' STATIONS IN X/Y DIRECTIONS'/
566 * 5X,' THE CERENKOV ARRAY IS CENTERED AROUND (0., 0.)'/)
567C CALCULATE CERENKOV YIELD FACTOR FROM WAVELENGTH BAND
568 IF ( WAVLGL .LT. 100.D0 .OR. WAVLGU .GT. 700.D0
569 * .OR. WAVLGL .GE. WAVLGU ) THEN
570 WRITE(MONIOU,*)'CERENKOV WAVELENGTH BAND FROM ',SNGL(WAVLGL),
571 * ' TO ',SNGL(WAVLGU),' NANOMETER'
572 WRITE(MONIOU,*)' IS OUT OF VALIDITY RANGE'
573 WRITE(MONIOU,*)'PLEASE READ THE MANUALS'
574 STOP
575 ENDIF
576 WRITE(MONIOU,*)'CERENKOV WAVELENGTH BAND FROM ',SNGL(WAVLGL),
577 * ' TO ',SNGL(WAVLGU),' NANOMETER'
578C WAVELENGTH IS CONVERTED FROM NM TO CM
579 CYIELD = (1.D0/WAVLGL - 1.D0/WAVLGU) * 2.D7 * PI / C(50)
580C CALCULATE FACTOR FOR ETA DENSITY NORML.(ETA AT SEA LEVEL = 0.283D-3)
581 ETADSN = 0.283D-3 * CATM(1) / BATM(1)
582
583 IF ( CERSIZ .GT. 0. ) THEN
584 WRITE(MONIOU,*)'CERENKOV BUNCH SIZE IS SET TO=',CERSIZ
585 ELSE
586 WRITE(MONIOU,*)'CERENKOV BUNCH SIZE IS CALCULATED FOR EACH ',
587 * 'SHOWER'
588 ENDIF
589
590 IF ( LCERFI ) THEN
591 WRITE(MONIOU,*)'CERENKOV PHOTONS ARE WRITTEN TO SEPARATE FILE'
592 ELSE
593 WRITE(MONIOU,*)'CERENKOV PHOTONS ARE WRITTEN TO PARTICLE ',
594 * 'OUTPUT FILE'
595 ENDIF
596
597c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
598c Next block of code has been modified, and is passed to MAIN
599c----------------------------------------------------------------------
600cC SCATTERING OF CENTER OF CHERENKOV ARRAY RELATIVE TO SHOWER AXIS
601c ICERML = MIN(MAX(ICERML,1),20)
602c XSCATT = ABS(XSCATT)
603c YSCATT = ABS(YSCATT)
604c IF ( ICERML .GE. 1 ) THEN
605c WRITE(MONIOU,5225)ICERML,XSCATT,YSCATT
606c 5225 FORMAT(' DEFINE MULTIPLE CERENKOV ARRAYS TO USE EACH',
607c * ' SHOWER SEVERAL TIMES'/ ' USE EACH EVENT ',I2,' TIMES'/
608c * ' THE EVENTS ARE SCATTERED QUASI RANDOMLY IN THE RANGE '/
609c * ' X = +- ',F10.0,' Y = +- ',F10.0)
610c DO 4438 I=1,ICERML
611c CALL SELCOR(CERXOS(I),CERYOS(I))
612c WRITE(MONIOU,4437) I,CERXOS(I),CERYOS(I)
613c 4437 FORMAT(' CORE OF EVENT ',I2,' AT ',2F12.2)
614c 4438 CONTINUE
615c XCMAX = XCMAX + XSCATT
616c YCMAX = YCMAX + YSCATT
617c ENDIF
618c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
619
620C STORE CERENKOV PARAMETERS IN EVENTHEADER
621 EVTH(86) = NCERX
622 EVTH(87) = NCERY
623 EVTH(88) = DCERX
624 EVTH(89) = DCERY
625 EVTH(90) = ACERX
626 EVTH(91) = ACERY
627 IF ( LCERFI ) THEN
628 EVTH(92) = 1.
629 ELSE
630 EVTH(92) = 0.
631 ENDIF
632 EVTH(96) = WAVLGL
633 EVTH(97) = WAVLGU
634
635c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
636c Next block of code has been passed to MAIN
637c----------------------------------------------------------------------
638c EVTH(98) = FLOAT(ICERML)
639c DO 480 I=1,20
640c EVTH( 98+I) = CERXOS(I)
641c EVTH(118+I) = CERYOS(I)
642c 480 CONTINUE
643c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
644
645C-----------------------------------------------------------------------
646C FLAG FOR ADDITIONAL MUON INFORMATION
647 IF ( FMUADD ) THEN
648 WRITE(MONIOU,*)
649 WRITE(MONIOU,*)'ADDITIONAL INFORMATION ON MUON ORIGIN IS',
650 * ' WRITTEN TO PARTICLE TAPE'
651 EVTH(94) = 1.
652 ELSE
653 EVTH(94) = 0.
654 ENDIF
655
656C PRINTOUT OF INFORMATIONS FOR DEBUGGING
657 IF ( DEBUG ) WRITE(MONIOU,484) MDEBUG
658 484 FORMAT(/' ATTENTION ! DEBUGGING IS ACTIVE'/
659 * ' ====> DEBUG INFORMATION WRITTEN TO UNIT ',I3//)
660
661c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
662c Next block of code is obsolete.
663c Now it's used "jcio" routines (C)
664cC-----------------------------------------------------------------------
665cC OPEN OUTPUT DATA SET FOR RUN
666c IBL = INDEX(DSN,' ')
667c DSN(IBL:73) = 'DAT000000'
668c WRITE(DSN(IBL+3:IBL+8),'(I6)') NRRUN
669c DO 274 L = IBL+3,IBL+8
670c IF ( DSN(L:L) .EQ. ' ' ) DSN(L:L) = '0'
671c 274 CONTINUE
672cC OPEN DATASET FOR PARTICLE OUTPUT
673c OPEN(UNIT=PATAPE,FILE=DSN,STATUS='NEW',
674c * FORM='UNFORMATTED',ACCESS='SEQUENTIAL')
675c WRITE(MONIOU,579) DSN
676c 579 FORMAT(/' PARTICLE OUTPUT TO DIRECTORY : ',A79)
677cC WRITE RUNHEADER TO OUTPUT BUFFER
678c CALL TOBUF( RUNH,0 )
679c
680cC OPEN OUTPUT DATA SET FOR CERENKOV PHOTONS
681c IF ( LCERFI ) THEN
682c DSN(IBL:73) = 'CER000000'
683c WRITE(DSN(IBL+3:IBL+8),'(I6)') NRRUN
684c DO 249 L = IBL+3,IBL+8
685c IF (DSN(L:L) .EQ. ' ' ) DSN(L:L) = '0'
686c 249 CONTINUE
687c OPEN(UNIT=CETAPE,FILE=DSN,STATUS='NEW',
688c * FORM='UNFORMATTED',ACCESS='SEQUENTIAL')
689c WRITE(MONIOU,580) DSN
690c 580 FORMAT(' CERENKOV OUTPUT TO DIRECTORY : ',A79)
691c CALL TOBUFC( RUNH,0 )
692c ELSE
693c WRITE(MONIOU,580) DSN
694c ENDIF
695cC RESET DSN
696c DSN(IBL:73) = ' '
697c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
698
699C OPEN ON EXTERNAL STACK
700C BLOCKS OF 32640 BYTES = 4080 REAL*8 = 340 PARTICLES
701 OPEN(UNIT=EXST,STATUS='SCRATCH',
702 * FORM='UNFORMATTED',ACCESS='DIRECT',RECL=MAXSTK)
703
704
705C-----------------------------------------------------------------------
706C WRITE DATA SET FOR INFORMATION BANK
707 IF (FDBASE ) THEN
708C OPEN OUTPUT DATA SET FOR RUN
709 IBL = INDEX(DSN,' ')
710 DSN(IBL:79) = 'DAT000000.dbase'
711 WRITE(DSN(IBL+3:IBL+8),'(I6)') NRRUN
712 DO 275 L = IBL+3,IBL+8
713 IF ( DSN(L:L) .EQ. ' ' ) DSN(L:L) = '0'
714 275 CONTINUE
715 OPEN(UNIT=MDBASE,FILE=DSN,STATUS='NEW')
716 WRITE(MONIOU,581) DSN
717 581 FORMAT(/' DBASE OUTPUT TO DIRECTORY : ',A79)
718C RESET DSN
719 DSN(IBL+9:IBL+14) = ' '
720
721 LSTDSN(1:3) = 'LST'
722 LSTDSN(4:9) = DSN(IBL+3:IBL+8)
723 VERVEN=FLOAT(IVERVN)/1000.
724 IF ( LLONGI ) THEN
725 ILONG = 1
726 ELSE
727 ILONG = 0
728 ENDIF
729 IF ( EVTH(75) .NE. 0. ) THEN
730 ISO = 0
731 ELSE
732 ISO = 1
733 ENDIF
734C SET DPMFLAG (0=VENUS, 1=HDPM, 2=SIBYLL, 3=QGSJET, 4=DPMJET)
735 IF ( EVTH( 76) .NE. 0. ) THEN
736 IDPM = 0
737 ELSEIF ( EVTH(139) .NE. 0. ) THEN
738 IDPM = 2
739 ELSEIF ( EVTH(141) .NE. 0. ) THEN
740 IDPM = 3
741 ELSEIF ( EVTH(143) .NE. 0. ) THEN
742 IDPM = 4
743 ELSE
744 IDPM = 1
745 ENDIF
746C INCREMENT DPMFLAG FOR VARIOUS CROSS SECTIONS
747C BY (0=HDPM-SIG, 10=VENUSSIG, 20=SIBYLLSIG, 30=QGSSIG, 4=DPMJETSIG)
748 IF ( EVTH(145) .NE. 0 ) THEN
749 IDPM = IDPM + 10
750 ELSEIF ( EVTH(140) .NE. 0 ) THEN
751 IDPM = IDPM + 20
752 ELSEIF ( EVTH(142) .NE. 0 ) THEN
753 IDPM = IDPM + 30
754 ELSEIF ( EVTH(144) .NE. 0 ) THEN
755 IDPM = IDPM + 40
756 ENDIF
757 MARK = '1'
758 LTHIN = .FALSE.
759 EFRAC = 0.D0
760
761 WRITE(MDBASE,666)VERNUM,MARK,MVDATE,VERVEN,
762 $INT(RUNH(3))+19000000,INT(EVTH(80)),INT(EVTH(79)),INT(EVTH(78)),
763 $INT(EVTH(77)),INT(RUNH(2)),INT(PRMPAR(1)),
764 $LLIMIT,ULIMIT,PSLOPE,INT(RUNH(20)),
765 $INT(RUNH(19)),INT(EVTH(76)),INT(EVTH(75)),ISO,IDPM,
766 $NFLAIN,NFLDIF,NFLPI0,NFLPIF,
767 $NFLCHE,NFRAGM,ILONG,THSTEP,
768 $BX,BZ,NOBSLV,
769 $OBSLEV(1),OBSLEV(2),OBSLEV(3),
770 $OBSLEV(4),OBSLEV(5),OBSLEV(6),
771 $OBSLEV(7),OBSLEV(8),OBSLEV(9),
772 $OBSLEV(10),ELCUT(1),ELCUT(2),ELCUT(3),
773 $ELCUT(4),EVTH(81),EVTH(82),EVTH(83),
774 $EVTH(84),FIXHEI,N1STTR,THICK0,
775 $STEPFC,ARRANG,INT(EVTH(94)),NSEQ,
776 $ISEED(1,1),ISEED(2,1),ISEED(3,1),ISEED(1,2),
777 $ISEED(2,2),ISEED(3,2),ISEED(1,3),
778 $ISEED(2,3),ISEED(3,3),0,DSN,
779 $LSTDSN,' JDD300.01',' JDD300.01',
780 $NSHOW,HOST,USER,LTHIN,EFRAC
781
782 666 FORMAT('#version#',F6.3,A1,'#versiondate#',I9,'#venusversion#',
783 $F6.3,'#rundate#',I9,/,'#computer#',I2,'#horizont#',I2,'#neutrino#'
784 $,I2,'#cerenkov#',I2,'#runnumber#',I7,/,'#primary#',I5,
785 $'#e_range_l#',E15.7,'#e_range_u#',E15.7,/,'#slope#',E15.7,'#nkg#',
786 $I2,'#egs#',I2,'#venus#',I2,'#gheisha#',I2,'#isobar#',I2,'#hdpm#',
787 $I2,/,'#hadflag1#',I2,'#hadflag2#',I2,'#hadflag3#',I2,'#hadflag4#',
788 $I2,'#hadflag5#',I2,'#hadflag6#',I2,/,'#longi#',I2,'#longistep#',
789 $E15.7,'#magnetx#',E15.7,/,'#magnetz#',E15.7,'#nobslev#',I3,/,
790 $'#obslev1#',E15.7,'#obslev2#',E15.7,'#obslev3#',E15.7,/,
791 $'#obslev4#',E15.7,'#obslev5#',E15.7,'#obslev6#',E15.7,/,
792 $'#obslev7#',E15.7,'#obslev8#',E15.7,'#obslev9#',E15.7,/,
793 $'#obslev10#',E15.7,'#hcut#',E15.7,'#mcut#',E15.7,/,'#ecut#',E15.7,
794 $'#gcut#',E15.7,'#thetal#',E15.7,/,'#thetau#',E15.7,'#phil#',E15.7,
795 $'#phiu#',E15.7,/,'#fixhei#',E15.7,'#n1sttr#',I3,'#fixchi#',E15.7,
796 $/,'#stepfc#',E15.7,'#arrang#',E15.7,'#muaddi#',I2,'#nseq#',I2,/,
797 $'#seq1seed1#',I9,'#seq1seed2#',I9,'#seq1seed3#',I9,/,'#seq2seed1#'
798 $,I9,'#seq2seed2#',I9,'#seq2seed3#',I9,/,'#seq3seed1#',I9,
799 $'#seq3seed2#',I9,'#seq3seed3#',I9,/,'#size#',I10,'#dsn_events#',
800 $A59,/,'#dsn_prtout# ',A9,'#tape_name#',A10,'#backup#',A10,/,
801 $'#howmanyshowers#',I10,'#host#',A20,'#user#',A20,/
802 $'#thinning#',L4,'#thinninglevel#',E15.7)
803
804C RESET DSN
805 DSN(IBL:IBL+14) = ' '
806 CLOSE(UNIT=MDBASE)
807 ENDIF
808
809 WRITE(MONIOU,*)'NUMBER OF SHOWERS TO GENERATE =',NSHOW
810 WRITE(MONIOU,*)
811 RETURN
812 END
Note: See TracBrowser for help on using the repository browser.