source: trunk/MagicSoft/Simulation/Corsika/Mmcs/box2.f

Last change on this file 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: 23.3 KB
Line 
1 SUBROUTINE BOX2
2
3C-----------------------------------------------------------------------
4C
5C DETERMINES POINT OF INTERACTION OR DECAY FOR ANY PARTICLE
6C HEAVY PRIMARIES AND STRANGE BARYONS INCLUDED
7C ANNIHILATION CROSS SECTION INCLUDED
8C PRECISE MEAN FREE PATH FOR DECAYING PARTICLES
9C HAS INTERACTION LENGTH STATISTICS INCLUDED
10C THIS SUBROUTINE IS CALLED FROM MAIN
11C-----------------------------------------------------------------------
12
13 IMPLICIT NONE
14*KEEP,AIR.
15 COMMON /AIR/ COMPOS,PROBTA,AVERAW,AVOGAD
16 DOUBLE PRECISION COMPOS(3),PROBTA(3),AVERAW,AVOGAD
17*KEEP,CHISTA.
18 COMMON /CHISTA/ IHYCHI,IKACHI,IMUCHI,INNCHI,INUCHI,IPICHI
19 INTEGER IHYCHI(124),IKACHI(124),IMUCHI(124),
20 * INNCHI(124),INUCHI(124),IPICHI(124)
21*KEEP,CONST.
22 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
23 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
24*KEEP,KAONS.
25 COMMON /KAONS/ CKA
26 DOUBLE PRECISION CKA(80)
27*KEEP,MUPART.
28 COMMON /MUPART/ AMUPAR,BCUT,CMUON,FMUBRM,FMUORG
29 DOUBLE PRECISION AMUPAR(14),BCUT,CMUON(11)
30 LOGICAL FMUBRM,FMUORG
31*KEEP,NCSNCS.
32 COMMON /NCSNCS/ SIGN30,SIGN45,SIGN60,SIGO30,SIGO45,SIGO60,
33 * SIGA30,SIGA45,SIGA60,PNOA30,PNOA45,PNOA60,
34 * SIG30A,SIG45A,SIG60A
35 DOUBLE PRECISION SIGN30(56),SIGN45(56),SIGN60(56),
36 * SIGO30(56),SIGO45(56),SIGO60(56),
37 * SIGA30(56),SIGA45(56),SIGA60(56),
38 * PNOA30(1540,3),PNOA45(1540,3),PNOA60(1540,3),
39 * SIG30A(56),SIG45A(56),SIG60A(56)
40*KEEP,OBSPAR.
41 COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
42 * THETPR,PHIPR,NOBSLV
43 DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
44 * THETAP,THETPR(2),PHIP,PHIPR(2)
45 INTEGER NOBSLV
46*KEEP,PAM.
47 COMMON /PAM/ PAMA,SIGNUM
48 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
49*KEEP,PARPAR.
50 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
51 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
52 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
53 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
54 INTEGER ITYPE,LEVL
55*KEEP,PARPAE.
56 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
57 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
58 * (CURPAR(4), PHI ), (CURPAR(5), H ),
59 * (CURPAR(6), T ), (CURPAR(7), X ),
60 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
61 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
62 * (CURPAR(12),ECM )
63*KEEP,RANDPA.
64 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
65 DOUBLE PRECISION FAC,U1,U2
66 REAL RD(3000)
67 INTEGER ISEED(103,10),NSEQ
68 LOGICAL KNOR
69*KEEP,REST.
70 COMMON /REST/ CONTNE,TAR,LT
71 DOUBLE PRECISION CONTNE(3),TAR
72 INTEGER LT
73*KEEP,RUNPAR.
74 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
75 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
76 * MONIOU,MDEBUG,NUCNUC,
77 * CETAPE,
78 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
79 * N1STTR,MDBASE,
80 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
81 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
82 * ,GHEISH,GHESIG
83 COMMON /RUNPAC/ DSN,HOST,USER
84 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
85 REAL STEPFC
86 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
87 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
88 * N1STTR,MDBASE
89 INTEGER CETAPE
90 CHARACTER*79 DSN
91 CHARACTER*20 HOST,USER
92
93 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
94 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
95 * ,GHEISH,GHESIG
96*KEEP,SIGM.
97 COMMON /SIGM/ SIGMA,SIGANN,SIGAIR,FRACTN,FRCTNO
98 DOUBLE PRECISION SIGMA,SIGANN,SIGAIR,FRACTN,FRCTNO
99*KEEP,STRBAR.
100 COMMON /STRBAR/ CSTRBA
101 DOUBLE PRECISION CSTRBA(11)
102*KEEP,VENUS.
103 COMMON /VENUS/ ISH0,IVERVN,MTAR99,FVENUS,FVENSG
104 INTEGER ISH0,IVERVN,MTAR99
105 LOGICAL FVENUS,FVENSG
106*KEND.
107
108 DOUBLE PRECISION CHIBRM,CHIPRM,CHIINT,CHI1,CHI2,CKA2,COR1,DH,
109 * ELAB,ELABLG,ELABT,FRAPTN,FRPTNO,
110 * HDEC,HEIGH,PLAB,PLABLG,SIGBRM,SIGPRM,
111 * SIG45,S45SQ,S4530,THICK
112 REAL EKIN,GBRSGM,GPRSGM
113 INTEGER I,IA,IHY,IP,KA,MU,NI,NU
114 EXTERNAL HEIGH,THICK,GBRSGM,GPRSGM
115 DOUBLE PRECISION SIGGHE,CGHSIG
116 EXTERNAL CGHSIG
117C-----------------------------------------------------------------------
118
119 IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,8)
120 444 FORMAT(' BOX2 : CURPAR=',1P,8E10.3)
121
122 ITYPE = CURPAR(1)
123
124 BETA = SQRT( GAMMA**2 - 1.D0 ) / GAMMA
125
126C-----------------------------------------------------------------------
127C PHOTONS, ELECTRONS,PI(0), AND ETA ARE TREATED SEPARATELY (SEE BOX3)
128 IF ( ITYPE .LE. 3 .OR. ITYPE .EQ. 7 .OR. ITYPE .EQ. 17 .OR.
129 * (ITYPE .GE. 71 .AND. ITYPE .LE. 74) ) THEN
130 CHI = 0.D0
131 RETURN
132 ENDIF
133
134C-----------------------------------------------------------------------
135C RESONANCES ARE TREATED SEPARATELY (SEE BOX3)
136 IF ( ITYPE .GT. 50 .AND. ITYPE .LE. 65 ) THEN
137 CHI = 0.D0
138 RETURN
139 ENDIF
140
141 THICKH = THICK(H)
142 ELAB = PAMA(ITYPE) * GAMMA
143
144C-----------------------------------------------------------------------
145C MU + , MU - DECAYS AFTER ITS LIFE TIME
146C MUON INTERACTS BY BREMSSTRAHLUNG OR PAIR PRODUCTION
147 IF ( ITYPE .EQ. 5 .OR. ITYPE .EQ. 6 ) THEN
148 CALL RMMAR( RD,3,1 )
149 COR1 = -LOG(RD(1)) * C(25) * C(19)
150 CALL PRANGE(COR1)
151 DH = H - HEIGH( THICKH + CHI*COSTHE )
152 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : ITYPE,RD(1),CHIDEC=',
153 * ITYPE,RD(1),SNGL(CHI)
154 IF ( GAMMA .LE. 200.D0 ) THEN
155 FDECAY = .TRUE.
156
157 ELSE
158C AT HIGHER ENERGIES CHECK FOR MUON BREMSSTRAHLUNG AND PAIR PRODUCTION
159 ELABLG = LOG(ELAB)
160C CALCULATE MUON BREMSSTRAHLUNG CROSS SECTION FOR AIR
161 IF ( ELAB .LE. 1.D5 ) THEN
162 FRACTN = COMPOS(1)*GBRSGM( 7.,SNGL(ELAB))
163 FRCTNO = FRACTN + COMPOS(2)*GBRSGM( 8.,SNGL(ELAB))
164 SIGBRM = FRCTNO + COMPOS(3)*GBRSGM(18.,SNGL(ELAB))
165 ELSE
166C PRELIMINARY PARAMETRIZED FOR ULTRAHIGH ENERGIES
167 SIGBRM = EXP( ELABLG * 0.04437D0 - 1.4805D0 )
168 FRACTN = SIGBRM * 0.78D0
169 FRCTNO = SIGBRM * 0.99D0
170 ENDIF
171 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : SIGBRM=',SNGL(SIGBRM)
172C CALCULATE MEAN FREE PATH FOR BREMSSTRAHLUNG
173 CHIBRM = -LOG(RD(2)) * AVERAW / (AVOGAD * SIGBRM)
174 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : ITYPE,RD(2),CHIBRM=',
175 * ITYPE,RD(2),SNGL(CHIBRM)
176 CHI1 = MIN( CHIBRM, CHI )
177
178 IF ( ELAB .LE. 1.D6 ) THEN
179C CALCULATE MUON PAIR PRODUCTION CROSS SECTION FOR AIR
180 FRAPTN = COMPOS(1)*GPRSGM( 7.,SNGL(ELAB))
181 FRPTNO = FRAPTN + COMPOS(2)*GPRSGM( 8.,SNGL(ELAB))
182 SIGPRM = FRPTNO + COMPOS(3)*GPRSGM(18.,SNGL(ELAB))
183 ELSE
184C PRELIMINARY PARAMETRIZED FOR ULTRAHIGH ENERGIES
185 SIGPRM = EXP( ELABLG * 0.2067D0 + 0.9169D0 )
186 FRACTN = SIGPRM * 0.78D0
187 FRCTNO = SIGPRM * 0.99D0
188 ENDIF
189 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : SIGPRM=',SNGL(SIGPRM)
190C CALCULATE MEAN FREE PATH FOR PAIR PRODUCTION
191 CHIPRM = -LOG(RD(3)) * AVERAW / (AVOGAD * SIGPRM)
192 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : ITYPE,RD(3),CHIPRM=',
193 * ITYPE,RD(3),SNGL(CHIPRM)
194 CHI2 = MIN( CHIPRM, CHI1 )
195 IF ( CHI2 .EQ. CHI ) THEN
196 FDECAY = .TRUE.
197 ELSEIF ( CHI2 .EQ. CHIBRM ) THEN
198 FDECAY = .FALSE.
199 FMUBRM = .TRUE.
200C TARGET IS CHOSEN AT RANDOM FOR MUON BREMSSTRAHLUNG
201 CALL RMMAR( RD,1,1 )
202 IF ( RD(1)*SIGBRM .LE. FRACTN ) THEN
203C BREMSSTRAHLUNG WITH NITROGEN
204 LT = 1
205 TAR = 14.D0
206 ELSEIF ( RD(1)*SIGBRM .LE. FRCTNO ) THEN
207C BREMSSTRAHLUNG WITH OXYGEN
208 LT = 2
209 TAR = 16.D0
210 ELSE
211C BREMSSTRAHLUNG WITH ARGON
212 LT = 3
213 TAR = 40.D0
214 ENDIF
215 ELSEIF ( CHI2 .EQ. CHIPRM ) THEN
216 FDECAY = .FALSE.
217 FMUBRM = .FALSE.
218C TARGET IS CHOSEN AT RANDOM FOR MUON PAIR PRODUCTION
219 CALL RMMAR( RD,1,1 )
220 IF ( RD(1)*SIGPRM .LE. FRAPTN ) THEN
221C PAIR PRODUCTION WITH NITROGEN
222 LT = 1
223 TAR = 14.D0
224 ELSEIF ( RD(1)*SIGPRM .LE. FRPTNO ) THEN
225C PAIR PRODUCTION WITH OXYGEN
226 LT = 2
227 TAR = 16.D0
228 ELSE
229C PAIR PRODUCTION WITH ARGON
230 LT = 3
231 TAR = 40.D0
232 ENDIF
233 ENDIF
234 CHI = CHI2
235 ENDIF
236
237C DECAY LENGTH STATISTICS
238 MU = 1.D0 + DH * 1.D-5 / COSTHE
239 MU = MIN( MU, 123 )
240 IMUCHI( MU) = IMUCHI( MU) + 1
241 IMUCHI(124) = IMUCHI(124) + 1
242
243C-----------------------------------------------------------------------
244C CHARGED PIONS
245 ELSEIF ( ITYPE .EQ. 8 .OR. ITYPE .EQ. 9 ) THEN
246 PLAB = ELAB * BETA
247C CALCULATION OF CROSS SECTION IN THE GHEISHA ROUTINES
248 IF ( GHEISH .AND. (ELAB .LE. HILOELB) ) THEN
249 EKIN = ELAB - PAMA(ITYPE)
250 SIGAIR = CGHSIG(SNGL(PLAB),EKIN,ITYPE)
251 GHESIG = .TRUE.
252 ELSE
253 GHESIG = .FALSE.
254C SIGMA IS ENERGY DEPENDENT INELASTIC PION-NUCLEON CROSS SECTION
255 IF ( PLAB .LE. 5.D0 ) THEN
256 SIGMA = 20.64D0
257 ELSEIF ( PLAB .LT. 1.D3 ) THEN
258 PLABLG = LOG(PLAB)
259C INELASTIC CROSS SECTIONS FROM PARTICLE DATA GROUP
260C (A.BALDINI ET AL., LANDOLT-BOERNSTEIN NEW SERIES I/12A (1987) 193)
261 SIGMA = 24.3D0 - 12.3D0 * PLAB**(-1.91D0)
262 * + 0.324D0 * PLABLG**2 - 2.44D0 * PLABLG
263 ELSE
264C FACTOR 0.6667 GIVES RATIO BETWEEN PION AND NUCLEON CROSS SECTION
265 SIGMA = 22.01D0 * ELAB**.0642D0 * 0.6667D0
266 ENDIF
267C AUXIL. QUANTITIES FOR INTERPOLATION
268 SIG45 = SIGMA - 45.D0
269 S45SQ = SIG45**2 / 450.D0
270 S4530 = SIG45 / 30.D0
271C INELASTIC CROSS SECTIONS OF AIR FOR PROJECTILE WITH MASS NUMBER 1
272 SIGAIR = (1.D0 - 2.D0 * S45SQ) * SIG45A(1)
273 * +(S45SQ - S4530) * SIG30A(1)
274 * +(S45SQ + S4530) * SIG60A(1)
275 ENDIF
276 IF ( DEBUG ) WRITE(MDEBUG,*)'BOX2 : SIGMA,SIGAIR,GHESIG=',
277 * SNGL(SIGMA),SNGL(SIGAIR),GHESIG
278
279 CALL RMMAR( RD,2,1 )
280C MEAN FREE PATH FOR INTERACTION (CHIINT) OR DECAY (CHI)
281 CHIINT = -LOG(RD(1)) * AVERAW / (AVOGAD * SIGAIR)
282 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : ITYPE,RD(1),CHIINT=',
283 * ITYPE,RD(1),SNGL(CHIINT)
284 COR1 = -LOG(RD(2)) * C(25) * C(18)
285 CALL PRANGE(COR1)
286 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : ITYPE,RD(2),CHIDEC=',
287 * ITYPE,RD(2),SNGL(CHI)
288 CHI = MIN( CHIINT, CHI )
289 IF ( CHI .LT. CHIINT ) THEN
290 FDECAY = .TRUE.
291 ELSE
292 FDECAY = .FALSE.
293 ENDIF
294
295C INTERACTION LENGTH STATISTICS
296 IP = 1.D0 + CHI * 0.1D0
297 IP = MIN( IP, 123 )
298 IPICHI( IP) = IPICHI( IP) + 1
299 IPICHI(124) = IPICHI(124) + 1
300
301C-----------------------------------------------------------------------
302C NUCLEONS AND ANTINUCLEONS
303 ELSEIF ( ITYPE .EQ. 13 .OR. ITYPE .EQ. 14 .OR.
304 * ITYPE .EQ. 15 .OR. ITYPE .EQ. 25 ) THEN
305 PLAB = ELAB * BETA
306C CALCULATION OF CROSS SECTION IN THE GHEISHA ROUTINES
307 IF ( GHEISH .AND. (ELAB .LE. HILOELB) ) THEN
308 EKIN = ELAB - PAMA(ITYPE)
309 SIGAIR = CGHSIG(SNGL(PLAB),EKIN,ITYPE)
310 GHESIG = .TRUE.
311 ELSE
312 GHESIG = .FALSE.
313C SIGMA IS ENERGY DEPENDENT INELASTIC NUCLEON-NUCLEON CROSS SECTION
314 IF ( PLAB .LT. 1.D1 ) THEN
315 SIGMA = 29.9D0
316 ELSEIF ( PLAB .LT. 1.D3 ) THEN
317 PLABLG = LOG(PLAB)
318C INELASTIC CROSS SECTIONS FROM PARTICLE DATA GROUP
319C (A.BALDINI ET AL., LANDOLT-BOERNSTEIN NEW SERIES I/12B (1987) 150)
320 SIGMA = 30.9D0 - 28.9D0 * PLAB**(-2.46D0)
321 * + 0.192D0 * PLABLG**2 - 0.835D0 * PLABLG
322 ELSE
323 SIGMA = 22.01D0 * ELAB**.0642D0
324 ENDIF
325
326C ADD ANNIHILATION CROSS SECTION FOR ANTI-NUCLEONS
327 IF ( ITYPE .EQ. 15 .OR. ITYPE .EQ. 25 ) THEN
328C ANNIHILATION CROSS SECTIONS FROM PARTICLE DATA GROUP
329C (A.BALDINI ET AL., LANDOLT-BOERNSTEIN NEW SERIES I/12B (1987) 286)
330 SIGANN = 0.532D0 + 0.634D2 * PLAB**(-0.71D0)
331 SIGMA = MIN( 120.D0, SIGMA + SIGANN )
332 ENDIF
333C AUXIL. QUANTITIES FOR INTERPOLATION
334 SIG45 = SIGMA - 45.D0
335 S45SQ = SIG45**2 / 450.D0
336 S4530 = SIG45 / 30.D0
337C INELASTIC CROSS SECTIONS OF AIR FOR PROJECTILE WITH MASS NUMBER 1
338 SIGAIR = (1.D0 - 2.D0 * S45SQ) * SIG45A(1)
339 * +(S45SQ - S4530) * SIG30A(1)
340 * +(S45SQ + S4530) * SIG60A(1)
341 ENDIF
342 IF ( DEBUG ) WRITE(MDEBUG,*)'BOX2 : SIGMA,SIGAIR,GHESIG=',
343 * SNGL(SIGMA),SNGL(SIGAIR),GHESIG
344
345C MEAN FREE PATH FROM MOLECULAR WEIGHT, AVOGADRO'S CONSTANT AND SIGMA
346 CALL RMMAR( RD,1,1 )
347 CHI = -LOG(RD(1)) * AVERAW / (AVOGAD * SIGAIR)
348 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : ITYPE,RD(1),CHI=',
349 * ITYPE,RD(1),SNGL(CHI)
350
351C INTERACTION LENGTH STATISTICS
352 NU = 1.D0 + CHI * 0.1D0
353 NU = MIN( NU, 123 )
354 INUCHI( NU) = INUCHI( NU) + 1
355 INUCHI(124) = INUCHI(124) + 1
356
357C-----------------------------------------------------------------------
358C KAONS (PARTICLE TYPES 10,11,12,16)
359 ELSEIF ( ITYPE .EQ. 10 .OR. ITYPE .EQ. 11 .OR.
360 * ITYPE .EQ. 12 .OR. ITYPE .EQ. 16 ) THEN
361 PLAB = ELAB * BETA
362C CALCULATION OF CROSS SECTION IN THE GHEISHA ROUTINES
363 IF ( GHEISH .AND. (ELAB .LE. HILOELB) ) THEN
364 EKIN = ELAB - PAMA(ITYPE)
365 SIGAIR = CGHSIG(SNGL(PLAB),EKIN,ITYPE)
366 GHESIG = .TRUE.
367 ELSE
368 GHESIG = .FALSE.
369C SIGMA IS ENERGY DEPENDENT INELASTIC KAON-NUCLEON CROSS SECTION
370 IF ( PLAB .LE. 1.D1 ) THEN
371 SIGMA = 14.11D0
372 ELSEIF ( PLAB .LT. 1.D3 ) THEN
373 PLABLG = LOG(PLAB)
374C INELASTIC CROSS SECTIONS FROM PARTICLE DATA GROUP
375C (A.BALDINI ET AL., LANDOLT-BOERNSTEIN NEW SERIES I/12B (1987) 56)
376 SIGMA = 12.3D0 - 7.77D0 * PLAB**(-2.12D0)
377 * + 0.0326D0 * PLABLG**2 + 0.738D0 * PLABLG
378 ELSE
379C FACTOR 0.5541 GIVES RATIO BETWEEN KAON AND NUCLEON CROSS SECTION
380 SIGMA = 22.01D0 * ELAB**.0642D0 * 0.5541D0
381 ENDIF
382C AUXIL. QUANTITIES FOR INTERPOLATION
383 SIG45 = SIGMA - 45.D0
384 S45SQ = SIG45**2 / 450.D0
385 S4530 = SIG45 / 30.D0
386C INELASTIC CROSS SECTIONS OF AIR FOR PROJECTILE WITH MASS NUMBER 1
387 SIGAIR = (1.D0 - 2.D0 * S45SQ) * SIG45A(1)
388 * +(S45SQ - S4530) * SIG30A(1)
389 * +(S45SQ + S4530) * SIG60A(1)
390 ENDIF
391 IF ( DEBUG ) WRITE(MDEBUG,*)'BOX2 : SIGMA,SIGAIR,GHESIG=',
392 * SNGL(SIGMA),SNGL(SIGAIR),GHESIG
393
394 CALL RMMAR( RD,2,1 )
395C MEAN FREE PATH FOR INTERACTION (CHIINT) OR DECAY (CHI)
396 CHIINT = -LOG(RD(1)) * AVERAW / (AVOGAD * SIGAIR)
397 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : ITYPE,RD(1),CHIINT=',
398 * ITYPE,RD(1),SNGL(CHIINT)
399
400 IF ( ITYPE .EQ. 16 ) THEN
401 CKA2 = CKA(5)
402 ELSEIF ( ITYPE .EQ. 10 ) THEN
403 CKA2 = CKA(6)
404 ELSE
405 CKA2 = CKA(3)
406 ENDIF
407 COR1 = -LOG(RD(2)) * C(25) * CKA2
408 IF ( SIGNUM(ITYPE) .EQ. 0.D0 ) THEN
409C NEUTRAL KAONS
410 DH = BETA * GAMMA * COSTHE * COR1
411 HDEC = MAX( H - DH, -1.D5 )
412 CHI = ( THICK(HDEC) - THICKH ) / COSTHE
413 ELSE
414C CHARGED KAONS
415 CALL PRANGE(COR1)
416 ENDIF
417
418 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : ITYPE,RD(2),CHIDEC=',
419 * ITYPE,RD(2),SNGL(CHI)
420 CHI = MIN( CHIINT, CHI )
421 IF ( CHI .LT. CHIINT ) THEN
422 FDECAY = .TRUE.
423 ELSE
424 FDECAY = .FALSE.
425 ENDIF
426
427C INTERACTION LENGTH STATISTICS
428 KA = 1.D0 + CHI * 0.1D0
429 KA = MIN( KA, 123 )
430 IKACHI( KA) = IKACHI( KA) + 1
431 IKACHI(124) = IKACHI(124) + 1
432
433C-----------------------------------------------------------------------
434C STRANGE BARYONS ( LAMBDA, SIGMA(+,0,-),XI(0,-), OMEGA- )
435 ELSEIF ( (ITYPE .GE. 18 .AND. ITYPE .LE. 24) .OR.
436 * (ITYPE .GE. 26 .AND. ITYPE .LE. 32) ) THEN
437 PLAB = ELAB * BETA
438C CALCULATION OF CROSS SECTION IN THE GHEISHA ROUTINES
439 IF ( GHEISH .AND. (ELAB .LE. HILOELB) ) THEN
440 EKIN = ELAB - PAMA(ITYPE)
441 SIGAIR = CGHSIG(SNGL(PLAB),EKIN,ITYPE)
442C SET CROSS SECTION VALUE TO A SMALL NUMBER FOR SIGMA0 AND ANTI SIGMA0
443 IF ( ITYPE .EQ. 20 .OR. ITYPE .EQ. 28 ) THEN
444 SIGAIR = 1.D-3
445 ENDIF
446 GHESIG = .TRUE.
447 ELSE
448 GHESIG = .FALSE.
449C CROSS SECTION FOR BARYONS IS ASSUMED TO BE THE SAME AS FOR NUCLEONS
450C SIGMA IS ENERGY DEPENDENT INELASTIC NUCLEON-NUCLEON CROSS SECTION
451 IF ( PLAB .LT. 1.D1 ) THEN
452 SIGMA = 29.9D0
453 ELSEIF ( PLAB .LT. 1.D3 ) THEN
454 PLABLG = LOG(PLAB)
455C INELASTIC CROSS SECTIONS FROM PARTICLE DATA GROUP
456C (A.BALDINI ET AL., LANDOLT-BOERNSTEIN NEW SERIES I/12B (1987) 150)
457 SIGMA = 30.9D0 - 28.9D0 * PLAB**(-2.46D0)
458 * + 0.192D0 * PLABLG**2 - 0.835D0 * PLABLG
459 ELSE
460 SIGMA = 22.01D0 * ELAB**.0642D0
461 ENDIF
462C AUXIL. QUANTITIES FOR INTERPOLATION
463 SIG45 = SIGMA - 45.D0
464 S45SQ = SIG45**2 / 450.D0
465 S4530 = SIG45 / 30.D0
466C INELASTIC CROSS SECTIONS OF AIR FOR PROJECTILE WITH MASS NUMBER 1
467 SIGAIR = (1.D0 - 2.D0 * S45SQ) * SIG45A(1)
468 * +(S45SQ - S4530) * SIG30A(1)
469 * +(S45SQ + S4530) * SIG60A(1)
470 ENDIF
471 IF ( DEBUG ) WRITE(MDEBUG,*)'BOX2 : SIGMA,SIGAIR,GHESIG=',
472 * SNGL(SIGMA),SNGL(SIGAIR),GHESIG
473
474 CALL RMMAR( RD,2,1 )
475C MEAN FREE PATH FOR INTERACTION (CHIINT) OR DECAY (CHI)
476 IF ( ITYPE .GE. 18 .AND. ITYPE .LE. 21 ) THEN
477 COR1 = -LOG(RD(2)) * C(25) * CSTRBA(ITYPE-17)
478 ELSEIF ( ITYPE .GE. 22 .AND. ITYPE .LE. 24 ) THEN
479 COR1 = -LOG(RD(2)) * C(25) * CSTRBA(ITYPE-15)
480 ELSEIF ( ITYPE .GE. 26 .AND. ITYPE .LE. 29 ) THEN
481 COR1 = -LOG(RD(2)) * C(25) * CSTRBA(ITYPE-25)
482 ELSEIF ( ITYPE .GE. 30 .AND. ITYPE .LE. 32 ) THEN
483 COR1 = -LOG(RD(2)) * C(25) * CSTRBA(ITYPE-23)
484 ENDIF
485 IF ( SIGNUM(ITYPE) .EQ. 0.D0 ) THEN
486C NEUTRAL STRANGE BARYONS
487 DH = BETA * GAMMA * COSTHE * COR1
488 HDEC = MAX( H - DH, -1.D5 )
489 CHI = ( THICK(HDEC) - THICKH ) / COSTHE
490 ELSE
491C CHARGED STRANGE BARYONS
492 CALL PRANGE(COR1)
493 ENDIF
494 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : ITYPE,RD(2),CHIDEC=',
495 * ITYPE,RD(2),SNGL(CHI)
496 CHIINT = -LOG(RD(1)) * AVERAW / (AVOGAD * SIGAIR)
497 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : ITYPE,RD(1),CHIINT=',
498 * ITYPE,RD(1),SNGL(CHIINT)
499 CHI = MIN( CHIINT, CHI )
500 IF ( CHI .LT. CHIINT ) THEN
501 FDECAY = .TRUE.
502 ELSE
503 FDECAY = .FALSE.
504 ENDIF
505C GHEISHA CANNOT TREAT SIGMA0 AND ANTI-SIGMA0, LET THEM DECAY
506 IF (GHESIG .AND. (ITYPE.EQ.20 .OR. ITYPE.EQ.28))FDECAY = .TRUE.
507
508C INTERACTION LENGTH STATISTICS
509 IHY = 1.D0 + CHI * 0.1D0
510 IHY = MIN( IHY, 123 )
511 IHYCHI(IHY) = IHYCHI(IHY) + 1
512 IHYCHI(124) = IHYCHI(124) + 1
513
514C-----------------------------------------------------------------------
515C HEAVY PRIMARIES ( ITYPE = 100 * A + Z , FE -> ITYPE = 5626 )
516C ( APPEARING AT FIRST INTERACTION AND AS REMANENTS OF THE PRIMARY )
517 ELSEIF ( ITYPE .GT. 100 ) THEN
518 IA = ITYPE / 100
519 IF ( IA .GT. 56 ) THEN
520 WRITE(MONIOU,*) 'BOX2 : UNEXPECTED PARTICLE TYPE=',ITYPE
521 STOP
522 ENDIF
523C MEAN FREE PATH OF THE HEAVY PRIMARY IS DEDUCED FROM THAT OF A NUCLEON
524C ONLY INELASTIC SCATTERING AT INTERACTIONS WITH HEAVY PRIMARY/FRAGMENT
525 ELAB = (PAMA(13) + PAMA(14)) * 0.5D0 * GAMMA
526 PLAB = ELAB * BETA
527C CALCULATION OF CROSS SECTION IN THE GHEISHA ROUTINES
528 ELABT = ELAB * IA
529
530c> *** modified by fs (22/09/98) *******************************
531
532c IF ( GHEISH .AND. (ELAB .LE. HILOELB) .AND.
533c * (ITYPE.EQ.402 .OR. ITYPE.EQ.201 .OR. ITYPE.EQ.301) ) THEN
534 IF ( GHEISH .AND. (ELAB .LE. HILOELB) .AND.
535 * (ITYPE.LE.101) ) THEN
536
537c> *** end of modification *************************************
538
539 EKIN = ELABT - PAMA(ITYPE)
540 SIGGHE = CGHSIG(SNGL(PLAB),EKIN,ITYPE)
541 IF ( SIGGHE .LE. 0. ) THEN
542 GHESIG = .FALSE.
543 ELSE
544 GHESIG = .TRUE.
545 SIGAIR = SIGGHE
546 ENDIF
547 ELSE
548 GHESIG = .FALSE.
549 ENDIF
550 IF ( .NOT. GHESIG ) THEN
551C SIGMA IS ENERGY DEPENDENT INELASTIC NUCLEON-NUCLEON CROSS SECTION
552 IF ( PLAB .LT. 1.D1 ) THEN
553 SIGMA = 29.9D0
554 ELSEIF ( PLAB .LT. 1.D3 ) THEN
555 PLABLG = LOG(PLAB)
556C INELASTIC CROSS SECTIONS FROM PARTICLE DATA GROUP
557C (A.BALDINI ET AL., LANDOLT-BOERNSTEIN NEW SERIES I/12B (1987) 150)
558 SIGMA = 30.9D0 - 28.9D0 * PLAB**(-2.46D0)
559 * + 0.192D0 * PLABLG**2 - 0.835D0 * PLABLG
560 ELSE
561 SIGMA = 22.01D0 * ELAB**.0642D0
562 ENDIF
563C AUXIL. QUANTITIES FOR INTERPOLATION
564 SIG45 = SIGMA - 45.D0
565 S45SQ = SIG45**2 / 450.D0
566 S4530 = SIG45 / 30.D0
567C INELASTIC CROSS SECTIONS OF AIR FOR PROJECTILE WITH MASS NUMBER IA
568 SIGAIR = (1.D0 - 2.D0 * S45SQ) * SIG45A(IA)
569 * +(S45SQ - S4530) * SIG30A(IA)
570 * +(S45SQ + S4530) * SIG60A(IA)
571 ENDIF
572 IF ( DEBUG ) WRITE(MDEBUG,*)'BOX2 : SIGMA,SIGAIR,GHESIG=',
573 * SNGL(SIGMA),SNGL(SIGAIR),GHESIG
574
575C MEAN FREE PATH FROM MOLECULAR WEIGHT, AVOGADRO'S CONSTANT AND SIGMA
576 IF ( SIGAIR .EQ. 0.D0 ) WRITE(MONIOU,*)
577 * 'BOX2: SIGAIR=0.D0, PROGRAM STOPPED',
578 * 'CHECK SELECTED CROSS SECTIONS AND PRIMARIES'
579 CALL RMMAR( RD,1,1 )
580 CHI = -LOG(RD(1)) * AVERAW / (AVOGAD * SIGAIR)
581 IF(DEBUG)WRITE(MDEBUG,*)'BOX2 : ITYPE,RD(1),CHI=',
582 * ITYPE,RD(1),SNGL(CHI)
583
584C INTERACTION LENGTH STATISTICS
585 NI = 1.D0 + CHI * 0.1D0
586 NI = MIN( NI, 123 )
587 INNCHI( NI) = INNCHI( NI) + 1
588 INNCHI(124) = INNCHI(124) + 1
589
590C-----------------------------------------------------------------------
591C ERROR IN PARTICLE CODE
592 ELSE
593 WRITE(MONIOU,*) 'BOX2 : UNEXPECTED PARTICLE TYPE=',ITYPE
594 STOP
595 ENDIF
596
597 RETURN
598 END
Note: See TracBrowser for help on using the repository browser.