source: trunk/MagicSoft/Simulation/Corsika/Mmcs/lepacx.f@ 9029

Last change on this file since 9029 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: 12.7 KB
Line 
1 SUBROUTINE LEPACX( ECMCE,SDMLOG,LEPART,IPART )
2
3C-----------------------------------------------------------------------
4C LE(ADING) PA(RTICLE) C(HARGE) (E)X(CHANGE)
5C
6C CONSIDERS CHARGE EXCHANGE POSSIBILITY OF (ANTI)LEADING PARTICLE
7C CONSIDERS RESONANCE EXCITATION WITHOUT/WITH CHARGE EXCHANGE
8C LASTPI INCREASED: CREATE ONE CHARGED PION FOR CHARGE CONSERVATION
9C LASTPI UNCHANGED: NO CHARGE EXCHANGE
10C LASTPI DECREASED: CANCEL ONE CHARGED PION FOR CHARGE CONSERVATION
11C NRESPC INCRESAED BY 1, IF PI(+-) WILL BE GENERATED BY RESON. DECAY
12C NRESPN INCRESAED BY 1, IF PI(0) WILL BE GENERATED BY RESON. DECAY
13C NCPLUS INCREASED BY 1, IF POSITIVE CHARGE IS CREATED
14C NCPLUS DECREASED BY 1, IF NEGATIVE CHARGE IS CREATED
15C THIS SUBROUTINE IS CALLED FROM HDPM
16C ARGUMENTS:
17C ECMCE = ENERGY FOR CHARGE EXCHANGE (ECMDPM OR ECMDIF)
18C SDMLOG = ELABLG FOR NSD, DMLOG FOR DIFFRACTION
19C LEPART = PARTICLE CODE OF (ANTI)LEADER EXCHANGING CHARGE
20C IPART = PARTICLE NUMBER IN ARRAY OF SECONDARY PARTICLES
21C = 1 FOR LEADER, = 2 FOR ANTI-LEADER
22C-----------------------------------------------------------------------
23
24 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
25*KEEP,CONST.
26 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
27 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
28*KEEP,LEPAR.
29 COMMON /LEPAR/ LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS
30 INTEGER LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS
31*KEEP,RANDPA.
32 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
33 DOUBLE PRECISION FAC,U1,U2
34 REAL RD(3000)
35 INTEGER ISEED(103,10),NSEQ
36 LOGICAL KNOR
37*KEEP,RESON.
38 COMMON /RESON/ RDRES,RESRAN,IRESPAR
39 REAL RDRES(2),RESRAN(1000)
40 INTEGER IRESPAR
41
42*KEEP,RUNPAR.
43 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
44 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
45 * MONIOU,MDEBUG,NUCNUC,
46 * CETAPE,
47 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
48 * N1STTR,MDBASE,
49 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
50 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
51 * ,GHEISH,GHESIG
52 COMMON /RUNPAC/ DSN,HOST,USER
53 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
54 REAL STEPFC
55 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
56 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
57 * N1STTR,MDBASE
58 INTEGER CETAPE
59 CHARACTER*79 DSN
60 CHARACTER*20 HOST,USER
61
62 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
63 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
64 * ,GHEISH,GHESIG
65*KEND.
66
67C-----------------------------------------------------------------------
68
69 IF ( DEBUG ) WRITE(MDEBUG,*) 'LEPACX: LEPART=',LEPART
70
71C SET PROBABILITIES FOR RESONANCE PRODUCTION (PRESPR) AND FOR
72C CHARGE EXCHANGE OR RESONANCE PRODUCTION (PCEXRS)
73 IF ( ECMCE .LE. 19.4D0 ) THEN
74 PCEXRS = 0.45D0
75 PRESPR = 0.35D0
76 ELSEIF ( ECMCE .LT. 968.5D0 ) THEN
77 PCEXRS = 0.45D0 + 0.034509D0 * (SDMLOG - 5.29832D0)
78 PRESPR = 0.0881897D0 * (SDMLOG - 5.29832D0)
79 ELSE
80 PCEXRS = 0.72D0
81 PRESPR = 0.69D0
82 ENDIF
83 PRESPR = MAX( 0.35D0, PRESPR )
84 IF ( LEPART .EQ. 7 ) THEN
85C ASSUME 50% CHARGE EXCHANGE FOR GAMMA INITIATED INTERACTION
86 PCEXRS = 0.5D0
87 PRESPR = 0.D0
88 ENDIF
89
90C THROW RANDOM NUMBER TO LOOK FOR RES. PRODUCTION OR CHARGE EXCHANGE
91 CALL RMMAR( RD,2,1 )
92
93C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94C RESONANCE IS FORMED. IF ADDITIONAL CHARGE EXCHANGE, THEN SET LASTPI
95 IF ( RD(1) .LE. PRESPR ) THEN
96
97C FIRST FOR NUCLEONS (AS MOST FREQUENT)
98 IF ( LEPART .EQ. 13 ) THEN
99 IF ( RD(2) .LE. 0.5 ) THEN
100C NEUTRON ----> DELTA(-)
101 LEPART = 57
102 NRESPC = NRESPC + 1
103 NCPLUS = NCPLUS - 1
104 ELSEIF ( RD(2) .GT. TB3 ) THEN
105C NEUTRON ----> DELTA(0)
106 LEPART = 56
107 CALL RMMAR( RDRES(IPART),1,1 )
108 IF ( RDRES(IPART) .LE. TB3 ) THEN
109 NRESPN = NRESPN + 1
110 ELSE
111 NRESPC = NRESPC + 1
112 LASTPI = LASTPI - 1
113 ENDIF
114 ELSE
115C NEUTRON ----> DELTA(+)
116 LEPART = 55
117 CALL RMMAR( RDRES(IPART),1,1 )
118 IF ( RDRES(IPART) .LE. TB3 ) THEN
119 NRESPN = NRESPN + 1
120 LASTPI = LASTPI - 1
121 ELSE
122 NRESPC = NRESPC + 1
123 ENDIF
124 NCPLUS = NCPLUS + 1
125 ENDIF
126 ELSEIF ( LEPART .EQ. 14 ) THEN
127 IF ( RD(2) .LE. 0.5 ) THEN
128C PROTON ----> DELTA(++)
129 LEPART = 54
130 NRESPC = NRESPC + 1
131 NCPLUS = NCPLUS + 1
132 ELSEIF ( RD(2) .GT. TB3 ) THEN
133C PROTON ----> DELTA(+)
134 LEPART = 55
135 CALL RMMAR( RDRES(IPART),1,1 )
136 IF ( RDRES(IPART) .LE. TB3 ) THEN
137 NRESPN = NRESPN + 1
138 ELSE
139 NRESPC = NRESPC + 1
140 LASTPI = LASTPI + 1
141 ENDIF
142 ELSE
143C PROTON ----> DELTA(0)
144 LEPART = 56
145 CALL RMMAR( RDRES(IPART),1,1 )
146 IF ( RDRES(IPART) .LE. TB3 ) THEN
147 NRESPN = NRESPN + 1
148 LASTPI = LASTPI + 1
149 ELSE
150 NRESPC = NRESPC + 1
151 ENDIF
152 NCPLUS = NCPLUS - 1
153 ENDIF
154
155C NOW FOR PIONS
156 ELSEIF ( LEPART .EQ. 8 .OR. LEPART .EQ. 9 ) THEN
157 IF ( RD(2) .LE. 0.5 ) THEN
158C PI(+-) ----> RHO(+-)
159 LEPART = LEPART + 44
160 NRESPN = NRESPN + 1
161 ELSE
162C PI(+-) ----> RHO(0) ( ----> PI(+) + PI(-) )
163 NCPLUS = NCPLUS + 2 * LEPART - 17
164 LEPART = 51
165 NRESPC = NRESPC + 1
166 ENDIF
167
168C NOW FOR KAONS
169 ELSEIF ( LEPART .EQ. 11 .OR. LEPART .EQ. 12 ) THEN
170 IF ( RD(2) .LE. 0.5 ) THEN
171C K(+-) ----> K*(+-)
172 LEPART = LEPART + 52
173 CALL RMMAR( RDRES(IPART),1,1 )
174 IF ( RDRES(IPART) .LE. TB3 ) THEN
175 NRESPN = NRESPN + 1
176 ELSE
177 NRESPC = NRESPC + 1
178 LASTPI = LASTPI + 1
179 ENDIF
180 ELSE
181C K(+) ----> K*(0)
182C K(-) ----> ANTI-K*(0)
183 CALL RMMAR( RDRES(IPART),1,1 )
184 NCPLUS = NCPLUS + 2 * LEPART - 23
185 IF ( RDRES(IPART) .LE. TB3 ) THEN
186 NRESPC = NRESPC + 1
187 ELSE
188 NRESPN = NRESPN + 1
189 LASTPI = LASTPI + 1
190 ENDIF
191 LEPART = 3*LEPART + 29
192 ENDIF
193 ELSEIF ( LEPART .EQ. 10 .OR. LEPART .EQ. 16 ) THEN
194 IF ( RD(2) .LE. 0.5 ) THEN
195C K(0) ----> (ANTI) K*(0)
196 CALL RMMAR( RD,1,1 )
197 IF ( RD(1) .LE. 0.5 ) THEN
198 LEPART = 62
199 ELSE
200 LEPART = 65
201 ENDIF
202 CALL RMMAR( RDRES(IPART),1,1 )
203 IF ( RDRES(IPART) .LE. TB3 ) THEN
204 NRESPC = NRESPC + 1
205 LASTPI = LASTPI - 1
206 ELSE
207 NRESPN = NRESPN + 1
208 ENDIF
209 ELSE
210C K(0) ----> K*(+-)
211 CALL RMMAR( RD,1,1 )
212 IF ( RD(1) .LE. 0.5 ) THEN
213 LEPART = 63
214 NCPLUS = NCPLUS + 1
215 ELSE
216 LEPART = 64
217 NCPLUS = NCPLUS - 1
218 ENDIF
219 CALL RMMAR( RDRES(IPART),1,1 )
220 IF ( RDRES(IPART) .LE. TB3 ) THEN
221 NRESPN = NRESPN + 1
222 LASTPI = LASTPI - 1
223 ELSE
224 NRESPC = NRESPC + 1
225 ENDIF
226 ENDIF
227
228C NOW FOR ANTINUCLEONS
229 ELSEIF ( LEPART .EQ. 25 ) THEN
230 IF ( RD(2) .LE. 0.5 ) THEN
231C ANTINEUTRON ----> ANTI-DELTA(0)
232 LEPART = 60
233 CALL RMMAR( RDRES(IPART),1,1 )
234 IF ( RDRES(IPART) .LE. TB3 ) THEN
235 NRESPN = NRESPN + 1
236 ELSE
237 NRESPC = NRESPC + 1
238 LASTPI = LASTPI - 1
239 ENDIF
240 ELSEIF ( RD(2) .GT. TB3 ) THEN
241C ANTINEUTRON ----> ANTI-DELTA(+)
242 LEPART = 61
243 NRESPC = NRESPC + 1
244 NCPLUS = NCPLUS + 1
245 ELSE
246C ANTINEUTRON ----> ANTI-DELTA(-)
247 LEPART = 59
248 CALL RMMAR( RDRES(IPART),1,1 )
249 IF ( RDRES(IPART) .LE. TB3 ) THEN
250 NRESPN = NRESPN + 1
251 LASTPI = LASTPI - 1
252 ELSE
253 NRESPC = NRESPC + 1
254 ENDIF
255 NCPLUS = NCPLUS - 1
256 ENDIF
257 ELSEIF ( LEPART .EQ. 15 ) THEN
258 IF ( RD(2) .LE. 0.5 ) THEN
259C ANTIPROTON ----> ANTI-DELTA(--)
260 LEPART = 58
261 NRESPC = NRESPC + 1
262 NCPLUS = NCPLUS - 1
263 ELSEIF ( RD(2) .GT. TB3 ) THEN
264C ANTIPROTON ----> ANTI-DELTA(-)
265 LEPART = 59
266 CALL RMMAR( RDRES(IPART),1,1 )
267 IF ( RDRES(IPART) .LE. TB3 ) THEN
268 NRESPN = NRESPN + 1
269 ELSE
270 NRESPC = NRESPC + 1
271 LASTPI = LASTPI + 1
272 ENDIF
273 ELSE
274C ANTIPROTON ----> ANTI-DELTA(0)
275 LEPART = 60
276 CALL RMMAR( RDRES(IPART),1,1 )
277 IF ( RDRES(IPART) .LE. TB3 ) THEN
278 NRESPN = NRESPN + 1
279 LASTPI = LASTPI + 1
280 ELSE
281 NRESPC = NRESPC + 1
282 ENDIF
283 NCPLUS = NCPLUS + 1
284 ENDIF
285
286 ELSEIF ( LEPART .EQ. 7 ) THEN
287C NO RESONANCE FORMATION FOR INDUCING GAMMA RADIATION
288 IF ( DEBUG ) WRITE(MDEBUG,*) 'LEPACX: NO EXCHANGE'
289
290 ELSEIF ( (LEPART .GE. 18 .AND. LEPART .LE. 24) .OR.
291 * (LEPART .GE. 26 .AND. LEPART .LE. 32) ) THEN
292C NO RESONANCE FORMATION FOR STRANGE BARYONS
293 IF ( DEBUG ) WRITE(MDEBUG,*) 'LEPACX: NO EXCHANGE'
294
295 ELSE
296 WRITE(MONIOU,100) LEPART
297 100 FORMAT(1H ,'LEPACX: UNIDENTIFIED PARTICLE CODE= ',I4,
298 * ' FOR RESONANCE FORMATION')
299 ENDIF
300 IF ( DEBUG ) WRITE(MDEBUG,102)
301 * LEPART,LASTPI,NRESPC,NRESPN,NCPLUS
302 102 FORMAT(' LEPACX: LEPART,LASTPI,NRESPC,NRESPN,NCPLUS=',5I5)
303
304C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
305C CHARGE EXCHANGE WITHOUT RESONANCE FORMATION
306 ELSEIF ( RD(1) .LE. PCEXRS ) THEN
307
308C FIRST FOR NUCLEONS (AS MOST FREQUENT)
309 IF ( LEPART .EQ. 13 ) THEN
310C NEUTRON ( + PI(+) ) ----> PROTON ( + PI(0) )
311 LEPART = 14
312 LASTPI = LASTPI - 1
313 NCPLUS = NCPLUS + 1
314 ELSEIF ( LEPART .EQ. 14 ) THEN
315C PROTON ( + PI(0) ) ----> NEUTRON ( + PI(+) )
316 LEPART = 13
317 LASTPI = LASTPI + 1
318 NCPLUS = NCPLUS - 1
319
320C NOW FOR PIONS
321 ELSEIF ( LEPART .EQ. 8 .OR. LEPART .EQ. 9 ) THEN
322C PI(+-) ----> PI(0)
323 NCPLUS = NCPLUS + 2 * LEPART - 17
324 LEPART = 7
325 LASTPI = LASTPI + 1
326
327C NOW FOR KAONS
328 ELSEIF ( LEPART .EQ. 11 .OR. LEPART .EQ. 12 ) THEN
329C K(+-) ----> K(0) (S OR L)
330 NCPLUS = NCPLUS + 2 * LEPART - 23
331 IF ( RD(2) .LE. 0.5 ) THEN
332 LEPART = 10
333 ELSE
334 LEPART = 16
335 ENDIF
336 LASTPI = LASTPI + 1
337 ELSEIF ( LEPART .EQ. 10 .OR. LEPART .EQ. 16 ) THEN
338C K(0) ----> K(+-)
339 IF ( RD(2) .LE. 0.5 ) THEN
340 LEPART = 11
341 NCPLUS = NCPLUS + 1
342 ELSE
343 LEPART = 12
344 NCPLUS = NCPLUS - 1
345 ENDIF
346 LASTPI = LASTPI - 1
347
348C NOW FOR ANTINUCLEONS
349 ELSEIF ( LEPART .EQ. 25 ) THEN
350C ANTINEUTRON ( + PI(-) ) ----> ANTIPROTON ( + PI(0) )
351 LEPART = 15
352 LASTPI = LASTPI - 1
353 NCPLUS = NCPLUS - 1
354 ELSEIF ( LEPART .EQ. 15 ) THEN
355C ANTIPROTON ( + PI(0) ) ----> ANTINEUTRON ( + PI(-) )
356 LEPART = 25
357 LASTPI = LASTPI + 1
358 NCPLUS = NCPLUS + 1
359
360C NOW FOR GAMMA INDUCED REACTIONS (ITYPE=7)
361 ELSEIF ( LEPART .EQ. 7 ) THEN
362C TEST IF CHARGE EXCHANGE REACTION FOR PI(0)
363C PI(0) ----> PI(+-)
364 IF ( RD(2) .LE. 0.5 ) THEN
365 LEPART = 8
366 NCPLUS = NCPLUS + 1
367 ELSE
368 LEPART = 9
369 NCPLUS = NCPLUS - 1
370 ENDIF
371 LASTPI = LASTPI - 1
372
373 ELSEIF ( (LEPART .GE. 18 .AND. LEPART .LE. 24) .OR.
374 * (LEPART .GE. 26 .AND. LEPART .LE. 32) ) THEN
375C NO CHARGE EXCHANGE FOR STRANGE BARYONS
376 IF ( DEBUG ) WRITE(MDEBUG,*) 'LEPACX: NO EXCHANGE'
377
378 ELSE
379 WRITE(MONIOU,101) LEPART
380 101 FORMAT(1H ,'LEPACX: UNIDENTIFIED PARTICLE CODE= ',I4,
381 * ' FOR CHARGE EXCHANGE')
382 ENDIF
383 IF ( DEBUG ) WRITE(MDEBUG,*) 'LEPACX: LEPART,LASTPI,NCPLUS=',
384 * LEPART,LASTPI,NCPLUS
385 ELSE
386 IF ( DEBUG ) WRITE(MDEBUG,*) 'LEPACX: NO EXCHANGE'
387 ENDIF
388
389 RETURN
390 END
Note: See TracBrowser for help on using the repository browser.