source: trunk/MagicSoft/Simulation/Corsika/Mmcs/nucint.f@ 10099

Last change on this file since 10099 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: 13.1 KB
Line 
1 SUBROUTINE NUCINT
2
3C-----------------------------------------------------------------------
4C NUC(LEAR) INT(ERACTION)
5C
6C SELECTS TYPE OF INTERACTION PROCESS FOR ISOBAR MODEL, ACCORDING TO ECM
7C ISOBAR MASSES INDEPENDENT OF RESPECTIVE ENERGY RANGES
8C HEAVY PRIMARIES AND STRANGE BARYONS INCLUDED
9C THIS SUBROUTINE IS CALLED FROM MAIN
10C-----------------------------------------------------------------------
11
12 IMPLICIT NONE
13*KEEP,AIR.
14 COMMON /AIR/ COMPOS,PROBTA,AVERAW,AVOGAD
15 DOUBLE PRECISION COMPOS(3),PROBTA(3),AVERAW,AVOGAD
16*KEEP,CONST.
17 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
18 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
19*KEEP,GENER.
20 COMMON /GENER/ GEN,ALEVEL
21 DOUBLE PRECISION GEN,ALEVEL
22*KEEP,IRET.
23 COMMON /IRET/ IRET1,IRET2
24 INTEGER IRET1,IRET2
25*KEEP,KAONS.
26 COMMON /KAONS/ CKA
27 DOUBLE PRECISION CKA(80)
28*KEEP,MULT.
29 COMMON /MULT/ EKINL,MSMM,MULTMA,MULTOT
30 DOUBLE PRECISION EKINL
31 INTEGER MSMM,MULTMA(37,13),MULTOT(37,13)
32*KEEP,PAM.
33 COMMON /PAM/ PAMA,SIGNUM
34 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
35*KEEP,PARPAR.
36 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
37 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
38 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
39 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
40 INTEGER ITYPE,LEVL
41*KEEP,PARPAE.
42 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
43 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
44 * (CURPAR(4), PHI ), (CURPAR(5), H ),
45 * (CURPAR(6), T ), (CURPAR(7), X ),
46 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
47 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
48 * (CURPAR(12),ECM )
49*KEEP,POLAR.
50 COMMON /POLAR/ POLART,POLARF
51 DOUBLE PRECISION POLART,POLARF
52*KEEP,RANDPA.
53 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
54 DOUBLE PRECISION FAC,U1,U2
55 REAL RD(3000)
56 INTEGER ISEED(103,10),NSEQ
57 LOGICAL KNOR
58*KEEP,RANGE.
59 COMMON /RANGE/ CC
60 DOUBLE PRECISION CC(20)
61*KEEP,RUNPAR.
62 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
63 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
64 * MONIOU,MDEBUG,NUCNUC,
65 * CETAPE,
66 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
67 * N1STTR,MDBASE,
68 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
69 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
70 * ,GHEISH,GHESIG
71 COMMON /RUNPAC/ DSN,HOST,USER
72 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
73 REAL STEPFC
74 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
75 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
76 * N1STTR,MDBASE
77 INTEGER CETAPE
78 CHARACTER*79 DSN
79 CHARACTER*20 HOST,USER
80
81 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
82 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
83 * ,GHEISH,GHESIG
84*KEEP,SIGM.
85 COMMON /SIGM/ SIGMA,SIGANN,SIGAIR,FRACTN,FRCTNO
86 DOUBLE PRECISION SIGMA,SIGANN,SIGAIR,FRACTN,FRCTNO
87*KEEP,STATI.
88 COMMON /STATI/ SABIN,SBBIN,INBIN,IPBIN,IKBIN,IHBIN
89 DOUBLE PRECISION SABIN(37),SBBIN(37)
90 INTEGER INBIN(37),IPBIN(37),IKBIN(37),IHBIN(37)
91*KEEP,VKIN.
92 COMMON /VKIN/ BETACM
93 DOUBLE PRECISION BETACM
94*KEND.
95
96 DOUBLE PRECISION BETA3,COSMU,COSTCM,COSTH3,GAMMA3,
97 * PHIMU,PHI3,WORK1,WORK2
98 INTEGER I,IGO,KJ
99C-----------------------------------------------------------------------
100
101 IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9)
102 444 FORMAT(' NUCINT: CURPAR=',1P,9E10.3)
103
104C SET GENERATION AND LEVEL OF LAST INTERACTION
105 SECPAR( 9) = GEN
106 SECPAR(10) = ALEVEL
107C RESET POLARIZATION, NOT USED FOR PARTICLES OTHER THAN MUONS YET
108 SECPAR(11) = 0.D0
109 SECPAR(12) = 0.D0
110
111C CALCULATE KIN. ENERGY BIN
112 EKINL = PAMA(ITYPE) * ( GAMMA - 1.D0 )
113 IF ( EKINL .GE. .1D0 ) THEN
114 KJ = INT( MIN( 37.D0, 4.D0 + 3.D0*LOG10(EKINL) ) )
115 ELSE
116 KJ = 1
117 ENDIF
118
119C-----------------------------------------------------------------------
120C PION INCIDENT
121 IF ( ITYPE .EQ. 8 .OR. ITYPE .EQ. 9 ) THEN
122 IF ( DEBUG ) WRITE(MDEBUG,*) 'NUCINT: PION EKINL=',SNGL(EKINL)
123 IPBIN(KJ) = IPBIN(KJ) + 1
124
125C DECAY OR INTERACTION FOR PIONS ?
126 IF ( FDECAY ) THEN
127 DO 10 I = 5,8
128 SECPAR(I) = CURPAR(I)
129 10 CONTINUE
130C DECAY PI(+,-) ----> MU(+,-) + (ANTI)-NEUTRINO(MU)
131 WORK1 = C(48) * GAMMA
132 WORK2 = C(49) * BETA * WORK1
133 CALL RMMAR( RD,2,1 )
134 COSTCM = 2.D0 * RD(1) - 1.D0
135 GAMMA3 = WORK1 + COSTCM * WORK2
136 BETA3 = SQRT( 1.D0 - 1.D0 / GAMMA3**2 )
137 COSTH3 = MIN( 1.D0, ( GAMMA * GAMMA3 - C(48) )
138 * /( BETA * GAMMA * BETA3 * GAMMA3 ) )
139 PHI3 = PI2 * RD(2)
140C MUON / NEUTRINO IS DROPPED
141 CALL ADDANG( COSTHE,PHI, COSTH3,PHI3, COSMU,PHIMU )
142 IF ( COSMU .GT. C(29) ) THEN
143C DIRECTION OF PION IN THE MUON CM SYSTEM (= DIRECTION OF POLARIZATION)
144C SEE: G. BARR ET AL., PHYS. REV. D39 (1989) 3532, EQ. 5
145C POLART IS COS OF ANGLE BETWEEN PION AND LABORATORY IN THE MU CM
146C POLARF IS ANGLE PHI AROUND THE LAB DIRECTION IN THE MU CM
147C POLART, POLARF ARE WITH RESPECT TO THE MU DIRECTION IN THE LAB SYSTEM
148 POLART = ( 2.D0*PAMA(8)*GAMMA*C(7)/(PAMA(5)*GAMMA3)
149 * - C(7) - 1.D0 ) / ( BETA3 * (1.D0 - C(7)) )
150 POLARF = PHI3 - PI
151C PION DIRECTION IS DIRECTION OF POLARIZATION FOR PI+, OPPOSITE FOR PI-
152 IF ( ITYPE .EQ. 9 ) THEN
153 POLART = -POLART
154 POLARF = POLARF + PI
155 ENDIF
156C GET THE POLARIZATION DIRECTION IN THE MU CM RELATIVE TO THE CORSIKA
157C COORDINATE SYSTEM
158 CALL ADDANG( COSMU,PHIMU, POLART,POLARF, POLART,POLARF )
159C MUON IS WRITTEN TO STACK
160 SECPAR( 1) = CURPAR(1) - 3.D0
161 SECPAR( 2) = GAMMA3
162 SECPAR( 3) = COSMU
163 SECPAR( 4) = PHIMU
164 SECPAR(11) = POLART
165 SECPAR(12) = POLARF
166 CALL TSTACK
167 SECPAR(11) = 0.D0
168 SECPAR(12) = 0.D0
169 ENDIF
170 IRET1 = 1
171 RETURN
172 ENDIF
173
174C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175C PION INTERACTS
176
177C CALCULATE GAMMA, BETA AND ENERGY IN CENTER OF MASS
178 ECM = SQRT( C(45) * GAMMA + C(46) )
179 GCM = (PAMA(ITYPE) * GAMMA + PAMA(14)) / ECM
180 BETACM = SQRT( 1.D0 - 1.D0 / GCM**2 )
181
182C LOW ENERGY HADRONIC INTERACTIONS
183C USE GHEISHA IF THE CROSS SECTION HAS BEEN CALCULATED FOR GHEISHA
184 IF ( GHEISH ) THEN
185 IF ( GHESIG ) THEN
186 CALL CGHEI
187 ELSE
188 CALL SDPM
189 ENDIF
190 ELSE
191C DETERMINE TYPE OF INTERACTION FOR PIONS
192 IF ( ECM .GT. CC(8) ) THEN
193C DUAL PARTON MODEL
194 CALL SDPM
195 ELSEIF ( ECM .GT. CC(7) ) THEN
196C HEAVY ISOBAR + HEAVY MESON
197 CALL BOX69
198 ELSEIF ( ECM .GT. CC(6) ) THEN
199 CALL RMMAR( RD,1,1 )
200 IF ( RD(1) .LE. 0.5 ) THEN
201C HEAVY ISOBAR + PION
202 CALL BOX68
203 ELSE
204C HEAVY MESON + NUCLEON
205 CALL BOX67
206 ENDIF
207 ELSEIF ( ECM .GT. CC(5) ) THEN
208C LIGHT ISOBAR + PION
209 CALL BOX66
210 ELSE
211C ELASTIC SCATTERING
212 CALL BOX65
213 ENDIF
214 ENDIF
215
216C-----------------------------------------------------------------------
217C NUCLEON OR ANTINUCLEON INCIDENT
218 ELSEIF ( ITYPE .EQ. 13 .OR. ITYPE .EQ. 14 .OR.
219 * ITYPE .EQ. 15 .OR. ITYPE .EQ. 25 ) THEN
220C CALCULATE GAMMA, BETA AND ENERGY IN CENTER OF MASS
221 GCM = SQRT( GAMMA * 0.5D0 + 0.5D0 )
222 ECM = PAMA(ITYPE) * GCM * 2.D0
223 BETACM = SQRT( 1.D0 - 1.D0 / GCM**2 )
224 IF ( DEBUG ) WRITE(MDEBUG,*) 'NUCINT: NUCL EKINL=',SNGL(EKINL)
225 INBIN(KJ) = INBIN(KJ) + 1
226
227C LOW ENERGY HADRONIC INTERACTIONS
228C USE GHEISHA IF THE CROSS SECTION HAS BEEN CALCULATED FOR GHEISHA
229 IF ( GHEISH ) THEN
230 IF ( GHESIG ) THEN
231 CALL CGHEI
232 ELSE
233 CALL SDPM
234 ENDIF
235 ELSE
236C DETERMINE TYPE OF INTERACTION FOR NUCLEONS AND ANTINUCLEONS
237 IF ( ECM .GT. CC(4) ) THEN
238C DUAL PARTON MODEL
239 CALL SDPM
240 ELSEIF ( ECM .GT. CC(3) ) THEN
241C USE THE INTERACTION ROUTINES OF PKF GRIEDER
242C 2 HEAVY ISOBARS AND ANNIHILATION
243 CALL BOX63
244 ELSEIF ( ECM .GT. CC(2) ) THEN
245C 1 HEAVY ISOBAR + NUCLEON AND ANNIHILATION
246 CALL BOX62
247 ELSEIF ( ECM .GT. CC(1) ) THEN
248C 1 LIGHT ISOBAR + NUCLEON AND ANNIHILATION
249 CALL BOX61
250 ELSE
251C ELASTIC SCATTERING AND ANNIHILATION
252 CALL BOX60
253 ENDIF
254 ENDIF
255
256C-----------------------------------------------------------------------
257C KAON INCIDENT
258 ELSEIF ( ITYPE .EQ. 11 .OR. ITYPE .EQ. 12 .OR.
259 * ITYPE .EQ. 10 .OR. ITYPE .EQ. 16 ) THEN
260 IF ( DEBUG ) WRITE(MDEBUG,*) 'NUCINT: KAON EKINL=',SNGL(EKINL)
261 IKBIN(KJ) = IKBIN(KJ) + 1
262
263C DECAY OR INTERACTION FOR KAONS ?
264 IF ( FDECAY ) THEN
265C KAON DECAYS. DETERMINE DECAY MODE FOR KAONS AND SET LIFE TIME
266 IF ( ITYPE .EQ. 10 ) THEN
267C K(0,L)-MESON
268 IGO = 4
269 ELSEIF ( ITYPE .EQ. 11 ) THEN
270C K(+)-MESON
271 IGO = 1
272 ELSEIF ( ITYPE .EQ. 12 ) THEN
273C K(-)-MESON
274 IGO = 2
275 ELSE
276C K(0,S)-MESON
277 IGO = 3
278 ENDIF
279 CALL KDECAY( IGO )
280 RETURN
281
282 ELSE
283C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
284C KAON INTERACTS
285C CALCULATE GAMMA, BETA AND ENERGY IN CENTER OF MASS
286 ECM = SQRT( CKA(13) * GAMMA + CKA(14) )
287 GCM = ( PAMA(ITYPE) * GAMMA + PAMA(14) ) / ECM
288 BETACM = SQRT( 1.D0 - 1.D0 / GCM**2 )
289C LOW ENERGY HADRONIC INTERACTIONS
290C USE GHEISHA IF THE CROSS SECTION HAS BEEN CALCULATED FOR GHEISHA
291 IF ( GHEISH ) THEN
292 IF ( GHESIG ) THEN
293 CALL CGHEI
294 ELSE
295 CALL SDPM
296 ENDIF
297 ELSE
298C KAON INTERACTS. DETERMINE TYPE OF INTERACTION FOR KAONS
299 IF ( ECM .GT. CC(12) ) THEN
300C DUAL PARTON MODEL
301 CALL SDPM
302C USE THE INTERACTION ROUTINES OF PKF GRIEDER
303 ELSEIF ( ECM .GT. CC(11) ) THEN
304C HEAVY ISOBAR + STRANGE MESON
305 CALL BOX74
306 ELSEIF ( ECM .GT. CC(10) ) THEN
307 CALL RMMAR( RD,1,1 )
308 IF ( RD(1) .GT. CKA(21) ) THEN
309C HEAVY ISOBAR + KAON
310 CALL BOX73
311 ELSE
312C STRANGE MESON + NUCLEON
313 CALL BOX72
314 ENDIF
315 ELSEIF ( ECM .GT. CC(9) ) THEN
316C LIGHT ISOBAR + KAON
317 CALL BOX71
318 ELSE
319C ELASTIC SCATTERING
320 CALL BOX70
321 ENDIF
322 ENDIF
323 ENDIF
324
325C-----------------------------------------------------------------------
326C STRANGE BARYON (LAMDA, SIGMA) INCIDENT
327 ELSEIF ( (ITYPE .GE. 18 .AND. ITYPE .LE. 24) .OR.
328 * (ITYPE .GE. 26 .AND. ITYPE .LE. 32) ) THEN
329 IF ( DEBUG ) WRITE(MDEBUG,*) 'NUCINT: SBAR EKINL=',SNGL(EKINL)
330 IHBIN(KJ) = IHBIN(KJ) + 1
331C DECAY OR INTERACTION FOR STRANGE BARYONS?
332 IF ( FDECAY ) THEN
333 CALL STRDEC
334 RETURN
335 ENDIF
336C CALCULATE GAMMA, BETA AND ENERGY IN CENTER OF MASS
337 ECM = SQRT( 2.D0 * PAMA(ITYPE) * PAMA(14) * GAMMA
338 * + PAMA(ITYPE)**2 + PAMA(14)**2 )
339 GCM = ( PAMA(ITYPE) * GAMMA + PAMA(14)) / ECM
340 BETACM = SQRT( 1.D0 - 1.D0 / GCM**2 )
341C LOW ENERGY HADRONIC INTERACTIONS
342C USE GHEISHA IF THE CROSS SECTION HAS BEEN CALCULATED FOR GHEISHA
343 IF ( GHEISH ) THEN
344 IF ( GHESIG ) THEN
345 CALL CGHEI
346 ELSE
347C VENUS TREATS STRANGE BARYONS
348 CALL SDPM
349 ENDIF
350 ELSE
351 IF ( ECM .GT. CC(4) ) THEN
352 CALL SDPM
353 ELSE
354C USE THE INTERACTION ROUTINES OF PKF GRIEDER
355C ELASTIC SCATTERING
356 CALL BOX60
357 ENDIF
358 ENDIF
359
360C-----------------------------------------------------------------------
361C HEAVY PRIMARY INCIDENT
362 ELSEIF ( ITYPE .GT. 100 ) THEN
363 IF (DEBUG) WRITE(MDEBUG,*) 'NUCINT: HEAVY PRIMARY EKINL=',
364 * SNGL(EKINL)
365C USE GHEISHA IF THE CROSS SECTION HAS BEEN CALCULATED FOR GHEISHA
366C (THIS MIGHT BE THE CASE FOR DEUTERONS, TRITONS AND ALPHAS)
367 IF ( GHEISH ) THEN
368 IF ( GHESIG ) THEN
369 CALL CGHEI
370 ELSE
371 CALL SDPM
372 ENDIF
373 ELSE
374C TREAT HEAVY PRIMARY IN SDPM
375 CALL SDPM
376 ENDIF
377
378C-----------------------------------------------------------------------
379C ILLEGAL PARTICLE
380 ELSE
381 WRITE(MONIOU,*) 'NUCINT: ILLEGAL PARTICLE = ',ITYPE
382 STOP
383 ENDIF
384
385C-----------------------------------------------------------------------
386C KILL PARTICLE
387 IRET1 = 1
388
389 RETURN
390 END
Note: See TracBrowser for help on using the repository browser.