source: trunk/MagicSoft/Simulation/Corsika/Mmcs/nihila.f@ 12392

Last change on this file since 12392 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.8 KB
Line 
1 SUBROUTINE NIHILA
2
3C-----------------------------------------------------------------------
4C (AN)NIHILA(TION)
5C
6C TREATES ANNIHILATION OF ANTINUCLEONS WITH FREE NUCLEONS
7C MOMENTA CONSERVED IN ALL DIRECTIONS
8C ENERGY CONSERVED BY MULTPLICATION OF ALL MOMENTA WITH A CORRECTION
9C FACTOR, CONSERVING MOMENTUM BALANCE
10C THIS SUBROUTINE IS CALLED FROM BOX60, BOX61, BOX62, AND BOX63
11C
12C DESIGN : D. HECK IK3 FZK KARLSRUHE
13C-----------------------------------------------------------------------
14
15 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
16*KEEP,ANNI.
17 COMMON /ANNI/ CAN,CANN
18 DOUBLE PRECISION CAN(50),CANN(50)
19*KEEP,CONST.
20 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
21 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
22*KEEP,MULT.
23 COMMON /MULT/ EKINL,MSMM,MULTMA,MULTOT
24 DOUBLE PRECISION EKINL
25 INTEGER MSMM,MULTMA(37,13),MULTOT(37,13)
26*KEEP,PAM.
27 COMMON /PAM/ PAMA,SIGNUM
28 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
29*KEEP,PARPAR.
30 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
31 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
32 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
33 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
34 INTEGER ITYPE,LEVL
35*KEEP,PARPAE.
36 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
37 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
38 * (CURPAR(4), PHI ), (CURPAR(5), H ),
39 * (CURPAR(6), T ), (CURPAR(7), X ),
40 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
41 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
42 * (CURPAR(12),ECM )
43*KEEP,RANDPA.
44 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
45 DOUBLE PRECISION FAC,U1,U2
46 REAL RD(3000)
47 INTEGER ISEED(103,10),NSEQ
48 LOGICAL KNOR
49*KEEP,RUNPAR.
50 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
51 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
52 * MONIOU,MDEBUG,NUCNUC,
53 * CETAPE,
54 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
55 * N1STTR,MDBASE,
56 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
57 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
58 * ,GHEISH,GHESIG
59 COMMON /RUNPAC/ DSN,HOST,USER
60 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
61 REAL STEPFC
62 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
63 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
64 * N1STTR,MDBASE
65 INTEGER CETAPE
66 CHARACTER*79 DSN
67 CHARACTER*20 HOST,USER
68
69 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
70 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
71 * ,GHEISH,GHESIG
72*KEEP,VKIN.
73 COMMON /VKIN/ BETACM
74 DOUBLE PRECISION BETACM
75*KEND.
76
77 DOUBLE PRECISION E(10),PHIPAR(10),PL(10),PTR(10),PTSQ(10)
78 DOUBLE PRECISION PX(10),PY(10)
79 DIMENSION ISEQ(10),NTYP(10)
80C-----------------------------------------------------------------------
81
82 IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9)
83 444 FORMAT(' NIHILA: CURPAR=',1P,9E10.3)
84
85 IREP = 0
86
87C RANDOM DECISION FOR ANNIHILATION WITH PROTON OR NEUTRON
88
89 20 CONTINUE
90 CALL RMMAR( RD,2,1 )
91 IREP = IREP + 1
92C AFTER THE 5TH TRY, QUIT THE ANNIHILATION WITHOUT ANY PION GENERATED
93 IF ( IREP .GT. 5 ) GOTO 999
94 IF ( RD(1) .LE. 0.5 ) THEN
95
96C-----------------------------------------------------------------------
97C ANTIPROTON - PROTON AND ANTINEUTRON - NEUTRON ANNIHILATION
98 IF ( RD(2) .LE. CANN(1) ) THEN
99C ANNIHILATION INTO PI+, PI-
100 NPIPOS = 1
101 NPINEG = 1
102 NPIZ = 0
103 ELSEIF ( RD(2) .LE. CANN(2) ) THEN
104C ANNIHILATION INTO PI+, PI-, PI0
105 NPIPOS = 1
106 NPINEG = 1
107 NPIZ = 1
108 ELSEIF ( RD(2) .LE. CANN(3) ) THEN
109C ANNIHILATION INTO PI+, PI-, 2PI0
110 NPIPOS = 1
111 NPINEG = 1
112 NPIZ = 2
113 ELSEIF ( RD(2) .LE. CANN(4) ) THEN
114C ANNIHILATION INTO PI+, PI-, 3PI0
115 NPIPOS = 1
116 NPINEG = 1
117 NPIZ = 3
118 ELSEIF ( RD(2) .LE. CANN(5) ) THEN
119C ANNIHILATION INTO PI+, PI-, 4PI0
120 NPIPOS = 1
121 NPINEG = 1
122 NPIZ = 4
123 ELSEIF ( RD(2) .LE. CANN(6) ) THEN
124C ANNIHILATION INTO 2PI+, 2PI-
125 NPIPOS = 2
126 NPINEG = 2
127 NPIZ = 0
128 ELSEIF ( RD(2) .LE. CANN(7) ) THEN
129C ANNIHILATION INTO 2PI+, 2PI-, PI0
130 NPIPOS = 2
131 NPINEG = 2
132 NPIZ = 1
133 ELSEIF ( RD(2) .LE. CANN(8) ) THEN
134C ANNIHILATION INTO 2PI+, 2PI-, 2PI0
135 NPIPOS = 2
136 NPINEG = 2
137 NPIZ = 2
138 ELSEIF ( RD(2) .LE. CANN(9) ) THEN
139C ANNIHILATION INTO 2PI+, 2PI-, 3PI0
140 NPIPOS = 2
141 NPINEG = 2
142 NPIZ = 3
143 ELSEIF ( RD(2) .LE. CANN(10) ) THEN
144C ANNIHILATION INTO 3PI+, 3PI-
145 NPIPOS = 3
146 NPINEG = 3
147 NPIZ = 0
148 ELSEIF ( RD(2) .LE. CANN(11) ) THEN
149C ANNIHILATION INTO 3PI+, 3PI-, PI0
150 NPIPOS = 3
151 NPINEG = 3
152 NPIZ = 1
153 ELSEIF ( RD(2) .LE. CANN(12) ) THEN
154C ANNIHILATION INTO 3PI+, 3PI-, 2PI0
155 NPIPOS = 3
156 NPINEG = 3
157 NPIZ = 2
158 ELSE
159C ANNIHILATION INTO 4PI0
160 NPIPOS = 0
161 NPINEG = 0
162 NPIZ = 4
163 ENDIF
164
165 ELSE
166C-----------------------------------------------------------------------
167C ANTIPROTON - NEUTRON (OR ANTINEUTRON - PROTON) ANNIHILATION
168 IF ( RD(2) .LE. CANN(13) ) THEN
169C ANNIHILATION INTO PI-, PI0
170 NPIPOS = 0
171 NPINEG = 1
172 NPIZ = 1
173 ELSEIF ( RD(2) .LE. CANN(14) ) THEN
174C ANNIHILATION INTO PI-, 2PI0
175 NPIPOS = 0
176 NPINEG = 1
177 NPIZ = 2
178 ELSEIF ( RD(2) .LE. CANN(15) ) THEN
179C ANNIHILATION INTO PI-, 3PI0
180 NPIPOS = 0
181 NPINEG = 1
182 NPIZ = 3
183 ELSEIF ( RD(2) .LE. CANN(16) ) THEN
184C ANNIHILATION INTO PI-, 4PI0
185 NPIPOS = 0
186 NPINEG = 1
187 NPIZ = 4
188 ELSEIF ( RD(2) .LE. CANN(17) ) THEN
189C ANNIHILATION INTO PI-, 5PI0
190 NPIPOS = 0
191 NPINEG = 1
192 NPIZ = 5
193 ELSEIF ( RD(2) .LE. CANN(18) ) THEN
194C ANNIHILATION INTO PI+, 2PI-
195 NPIPOS = 1
196 NPINEG = 2
197 NPIZ = 0
198 ELSEIF ( RD(2) .LE. CANN(19) ) THEN
199C ANNIHILATION INTO PI+, 2PI-, PI0
200 NPIPOS = 1
201 NPINEG = 2
202 NPIZ = 1
203 ELSEIF ( RD(2) .LE. CANN(20) ) THEN
204C ANNIHILATION INTO PI+, 2PI-, 2PI0
205 NPIPOS = 1
206 NPINEG = 2
207 NPIZ = 2
208 ELSEIF ( RD(2) .LE. CANN(21) ) THEN
209C ANNIHILATION INTO PI+, 2PI-, 3PI0
210 NPIPOS = 1
211 NPINEG = 2
212 NPIZ = 3
213 ELSEIF ( RD(2) .LE. CANN(22) ) THEN
214C ANNIHILATION INTO PI+, 2PI-, 4PI0
215 NPIPOS = 1
216 NPINEG = 2
217 NPIZ = 4
218 ELSEIF ( RD(2) .LE. CANN(23) ) THEN
219C ANNIHILATION INTO 2PI+, 3PI-
220 NPIPOS = 2
221 NPINEG = 3
222 NPIZ = 0
223 ELSEIF ( RD(2) .LE. CANN(24) ) THEN
224C ANNIHILATION INTO 2PI+, 3PI-, PI0
225 NPIPOS = 2
226 NPINEG = 3
227 NPIZ = 1
228 ELSEIF ( RD(2) .LE. CANN(25) ) THEN
229C ANNIHILATION INTO 2PI+, 3PI-, 2PI0
230 NPIPOS = 2
231 NPINEG = 3
232 NPIZ = 2
233 ELSEIF ( RD(2) .LE. CANN(26) ) THEN
234C ANNIHILATION INTO 2PI+, 3PI-, 3PI0
235 NPIPOS = 2
236 NPINEG = 3
237 NPIZ = 3
238 ELSE
239C ANNIHILATION INTO 3PI+, 4PI-
240 NPIPOS = 3
241 NPINEG = 4
242 NPIZ = 0
243 ENDIF
244
245C CHARGE INVERSION IF ANTINEUTRON ANNIHILATES WITH PROTON
246 IF ( ITYPE .EQ. 25 ) THEN
247 NPINEG = NPINEG - 1
248 NPIPOS = NPIPOS + 1
249 ENDIF
250
251 ENDIF
252
253 NPI = NPIPOS + NPINEG + NPIZ
254 FNPI = 1.D0 / NPI
255 GCMI = 1.D0 / GCM
256
257C-----------------------------------------------------------------------
258C CHARGE ASSIGNMENT
259
260 DO 26 I = 1,NPI
261 IF ( I .LE. NPIZ ) THEN
262C NEUTRAL PIONS
263 NTYP(I) = 7
264 ELSEIF ( I .LE. NPIZ+NPIPOS ) THEN
265C POSITIVE PIONS
266 NTYP(I) = 8
267 ELSE
268C NEGATIVE PIONS
269 NTYP(I) = 9
270 ENDIF
271 26 CONTINUE
272
273C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274C KINEMATIC CALCULATIONS
275
276 ISCALE = 0
277 27 CONTINUE
278 ISCALE = ISCALE + 1
279C AFTER THE 5TH TRY, TAKE A NEW SET OF PIONS
280 IF ( ISCALE .GT. 5 ) GOTO 20
281
282C DISTRIBUTUION OF TRANSVERSE MOMENTA PTR
283
284 CORECT = 5.3333333333D0 * FNPI**1.5D0 * ECM
285 DO 28 I = 1,NPI
286 PTR(I) = PTRANS(DUMMY) * (1.33333333D0 + CORECT)
287 28 CONTINUE
288
289 SUMPX = 0.D0
290 SUMPY = 0.D0
291 CALL RMMAR( RD,NPI,1 )
292 DO 29 I = 1,NPI
293C SELECT EMISSION ANGLE BY REDUCED RESIDUAL DIRECTION
294 IF ( SUMPX .NE. 0.D0 .OR. SUMPY .NE. 0.D0 ) THEN
295 PHISUM = ATAN2( SUMPY, SUMPX )
296 ELSE
297 PHISUM = 0.D0
298 ENDIF
299 PHIPAR(I) = PHISUM + PI + PI * (NPI+1-I) * (2.*RD(I)-1.) * FNPI
300 PX(I) = COS( PHIPAR(I) ) * PTR(I)
301 PY(I) = SIN( PHIPAR(I) ) * PTR(I)
302 SUMPX = SUMPX + PX(I)
303 SUMPY = SUMPY + PY(I)
304 29 CONTINUE
305
306C CORRECTION OF TRANSVERSE MOMENTA TO KEEP TRANSVERSE MOMENTUM BALANCE
307 SUMPT2 = ECM
308 DPX = SUMPX * FNPI
309 DPY = SUMPY * FNPI
310 DO 30 I = 1,NPI
311 PX(I) = PX(I) - DPX
312 PY(I) = PY(I) - DPY
313 IF ( PX(I) .NE. 0.D0 .OR. PY(I) .NE. 0.D0 ) THEN
314 PHIPAR(I) = ATAN2( PY(I), PX(I) )
315 ELSE
316 PHIPAR(I) = 0.D0
317 ENDIF
318 PTSQ(I) = PX(I)**2 + PY(I)**2
319 SUMPT2 = SUMPT2 - SQRT( PAMA(NTYP(I))**2 + PTSQ(I) )
320 30 CONTINUE
321
322C CHECK, IF C.M. ENERGY IS EXHAUSTED BY TRANSVERSE MOMENTA
323C IF SO, TRY ANOTHER SET OF TRANSVERSE MOMENTA
324 IF ( SUMPT2 .LE. 0.D0 ) GOTO 27
325
326C DISTRIBUTION OF LONGITUDINAL MOMENTA PL
327
328C SUM1PL IS SUM OF ABS. VALUES OF LONGITUDINAL MOMENTA
329 F = SUMPT2 * FNPI
330 SUM1PL = 0.D0
331 DO 31 I = 1,NPI
332 FWHM = F + 0.5D0 * SQRT( PTSQ(I) )
333 PL(I) = ABS( RANNOR(0.D0,FWHM) )
334 SUM1PL = SUM1PL + PL(I)
335C SET SEQUENCE COUNTER
336 ISEQ(I) = I
337 31 CONTINUE
338
339C SORT ISEQ IN DECREASING SIZE OF THE LONGITUDINAL MOMENTUM
340 DO 33 I = 1,NPI
341 DO 32 K = I+1,NPI
342 IF ( PL(ISEQ(I)) .LT. PL(ISEQ(K)) ) THEN
343 IHELP = ISEQ(I)
344 ISEQ(I) = ISEQ(K)
345 ISEQ(K) = IHELP
346 ENDIF
347 32 CONTINUE
348 33 CONTINUE
349
350C TRY TO BALANCE LONGITUDINAL MOMENTA (TO MINIMIZE CORRECTIONS)
351C START WITH LONG. MOMENTA IN FORWARD/BACKWARD DIRECTION BY RANDOM
352 CALL RMMAR( RD,1,1 )
353 IF ( RD(1) .LE. 0.5 ) PL(ISEQ(1)) = - PL(ISEQ(1))
354 SUMPL = PL(ISEQ(1))
355 DO 34 I = 2,NPI
356 SUM1PL = SUM1PL - PL(ISEQ(I))
357C IF THERE IS NOT ENOUGH MOMENTUM LEFT, SELECT FORWARD/BACKWARD TO
358C BALANCE MOMENTUM, ELSE CHOOSE DIRECTION BY RANDOM
359 IF ( PL(ISEQ(I))+ABS(SUMPL) .GT. SUM1PL ) THEN
360 IF ( PL(ISEQ(I))*SUMPL .GT. 0.D0 ) PL(ISEQ(I)) = -PL(ISEQ(I))
361 ELSE
362 CALL RMMAR( RD,1,1 )
363 IF ( RD(1) .LE. 0.5 ) PL(ISEQ(I)) = - PL(ISEQ(I))
364 ENDIF
365 SUMPL = SUMPL + PL(ISEQ(I))
366 34 CONTINUE
367
368C CORRECTION OF LONGITUDINAL MOMENTA TO KEEP MOMENTUM BALANCE
369 DPL = SUMPL * FNPI
370 DO 35 I = 1,NPI
371 PL(I) = PL(I) - DPL
372 35 CONTINUE
373
374C ITERATIVE CORRECTION OF ALL MOMENTA TO KEEP ENERGY BALANCE
375
376 IREPET = 0
377 36 CONTINUE
378 IREPET = IREPET + 1
379 IF ( IREPET .GT. 10 ) GOTO 27
380 ETOT = 0.D0
381C CHECK ENERGY CONSERVATION
382 DO 37 I = 1,NPI
383 PTSQ(I) = PX(I)**2 + PY(I)**2
384 E(I) = SQRT( PAMA(NTYP(I))**2 + PTSQ(I) + PL(I)**2 )
385 ETOT = ETOT + E(I)
386 37 CONTINUE
387
388 ECORR = ECM / ETOT - 1.D0
389
390C LOOK WHETHER ENERGY IS CONSERVED WITHIN 1 %
391 IF ( ABS(ECORR) .GT. .01D0 ) THEN
392C FACTOR IS MODIFIED WITH EMPIRICAL TERM 1/GCM FOR FASTER CONVERGENCE
393 FACT = (0.5D0+GCMI) * ECORR * 0.02D0 * NPI
394 DO 38 I = 1,NPI
395 PX(I) = PX(I) * ( FACT + 1.D0 )
396 PY(I) = PY(I) * ( FACT + 1.D0 )
397 PL(I) = PL(I) * ( FACT * 20.D0 + 1.D0 )
398 38 CONTINUE
399 GOTO 36
400 ENDIF
401
402C LORENTZ TRANSFORMATION FROM C.M. TO LAB. FRAME
403
404 DO 40 K = 5,8
405 SECPAR(K) = CURPAR(K)
406 40 CONTINUE
407 DO 41 I = 1,NPI
408 PLLAB = GCM * ( PL(I) + BETACM * E(I) )
409 IF ( PLLAB .LE. 0.D0 ) GOTO 41
410 CTHETA = PLLAB / SQRT( PTSQ(I) + PLLAB**2 )
411 IF ( CTHETA .LT. C(27) ) GOTO 41
412 CALL ADDANG( COSTHE,PHI, CTHETA,PHIPAR(I), SECPAR(3),SECPAR(4) )
413 IF ( SECPAR(3) .LT. C(29) ) GOTO 41
414 SECPAR(1) = NTYP(I)
415 SECPAR(2) = GCM / PAMA(NTYP(I)) * ( PL(I) * BETACM + E(I) )
416 CALL TSTACK
417 41 CONTINUE
418
419 MSMM = MSMM + NPI
420 999 CONTINUE
421
422 RETURN
423 END
Note: See TracBrowser for help on using the repository browser.