source: branches/start/MagicSoft/Simulation/Corsika/Mmcs/parnum.f@ 18569

Last change on this file since 18569 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: 10.0 KB
Line 
1 SUBROUTINE PARNUM( INUMFL )
2
3C-----------------------------------------------------------------------
4C PART(ICLE TYPE) NUM(BERS)
5C
6C DETERMINES THE NUMBERS OF SECONDARY PARTICLES FOR EACH TYPE
7C THIS SUBROUTINE IS CALLED FROM HDPM
8C ARGUMENT
9C INUMFL = 0 CORRECT DETERMINATION OF PARTICLE NUMBERS
10C = 1 SOMETHING WENT WRONG WITH NEUTRAL PARTICLE NUMBERS
11C-----------------------------------------------------------------------
12
13 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
14*KEEP,CONST.
15 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
16 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
17*KEEP,EDECAY.
18 COMMON /EDECAY/ CETA
19 DOUBLE PRECISION CETA(5)
20*KEEP,INDICE.
21 COMMON /INDICE/ NNUCN,NKA0,NHYPN,NETA,NETAS,NPIZER,
22 * NNC,NKC,NHC,NPC,NCH,NNN,NKN,NHN,NET,NPN
23 INTEGER NNUCN(2:3),NKA0(2:3),NHYPN(2:3),NETA(2:3,1:4),
24 * NETAS(2:3),NPIZER(2:3),
25 * NNC,NKC,NHC,NPC,NCH,NNN,NKN,NHN,NET,NPN
26*KEEP,INTER.
27 COMMON /INTER/ AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB,
28 * FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3,
29 * RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG,
30 * WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN,
31 * IDIF,ITAR
32 DOUBLE PRECISION AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB,
33 * FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3,
34 * RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG,
35 * WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN
36 INTEGER IDIF,ITAR
37*KEEP,LEPAR.
38 COMMON /LEPAR/ LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS
39 INTEGER LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS
40*KEEP,RANDPA.
41 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
42 DOUBLE PRECISION FAC,U1,U2
43 REAL RD(3000)
44 INTEGER ISEED(103,10),NSEQ
45 LOGICAL KNOR
46*KEEP,RATIOS.
47 COMMON /RATIOS/ RPI0R,RPIER,RPEKR,RPEKNR,PPICH,PPINCH,PPNKCH,
48 * ISEL,NEUTOT,NTOTEM
49 DOUBLE PRECISION RPI0R,RPIER,RPEKR,RPEKNR,PPICH,PPINCH,PPNKCH
50 INTEGER ISEL,NEUTOT,NTOTEM
51*KEEP,RUNPAR.
52 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
53 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
54 * MONIOU,MDEBUG,NUCNUC,
55 * CETAPE,
56 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
57 * N1STTR,MDBASE,
58 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
59 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
60 * ,GHEISH,GHESIG
61 COMMON /RUNPAC/ DSN,HOST,USER
62 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
63 REAL STEPFC
64 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
65 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
66 * N1STTR,MDBASE
67 INTEGER CETAPE
68 CHARACTER*79 DSN
69 CHARACTER*20 HOST,USER
70
71 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
72 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
73 * ,GHEISH,GHESIG
74*KEND.
75
76 REAL RDETA
77C-----------------------------------------------------------------------
78
79 IF ( DEBUG ) WRITE(MDEBUG,*) 'PARNUM: NCH,NEUTOT,NTOTEM=',
80 * NCH,NEUTOT,NTOTEM
81
82 INUMFL = 0
83C RESET PARTICLE NUMBERS
84 NNC = 0
85 NKC = 0
86 NHC = 0
87 NPC = 0
88C ISEL IS 1 MEANS VERY LOW MULTIPLICITY
89C CREATE ONLY PIONS (TO RISKY TO CREATE OTHER PARTICLES)
90 IF ( ISEL .EQ. 1 ) THEN
91 NNN = 0
92 NKN = 0
93 NET = 0
94 NHN = 0
95 NPN = 0
96 NETAS(2) = 0
97 NETAS(3) = 0
98C CREATE RANDOM NUMBERS
99 CALL RMMAR( RD,NTOTEM,1 )
100 DO 1000 I = 1,NTOTEM
101 IF ( RD(I) .LE. TB3 ) THEN
102 NPC = NPC + 1
103 ELSE
104 NPN = NPN + 1
105 ENDIF
106 1000 CONTINUE
107C NO NEUTRAL PARTICLES FOR THE 3RD STRING EXCEPT EVENTUALLY PI(0)
108 NNUCN(3) = 0
109 NKA0(3) = 0
110 NHYPN(3) = 0
111 NETAS(3) = 0
112 NPIZER(3) = MAX( 0, NINT(RC3TO2/(1.D0+RC3TO2)*DBLE(NPN)) )
113 IF ( DEBUG ) WRITE(MDEBUG,*) ' ISEL=1, NTOTEM=',NTOTEM
114
115 ELSE
116
117C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118C NOW THE CASE OF HAVING ENOUGH PARTICLES TO BE ABLE TO CREATE
119C KAONS, NUCLEONS, AND HYPERONS TOO.
120
121C ...FOR NEUTRALS
122 NCOUNT = 0
123C BEGIN OF REJECT LOOP
124 1002 K = 1
125 CALL RMMAR( RD,NEUTOT+3,1 )
126C DETERMINE NUMBER OF PI(0), ETA, K0S/K0 PAIRS, NEUTRON/ANTINEUTRON
127C PAIRS, AND NEUTRAL HYPERON PAIRS AND SUM UP THE GAMMAS
128C FOR 1ST + 2ND STRING: J IS 2; FOR 3RD STRING: J IS 3
129 SGAMMA = 0.D0
130 DO 1010 J = 2,3
131 NNUCN(J) = 0
132 NKA0(J) = 0
133 NHYPN(J) = 0
134 NETA(J,1) = 0
135 NETA(J,2) = 0
136 NETA(J,3) = 0
137 NETA(J,4) = 0
138 NPIZER(J) = 0
139 IF ( J .EQ. 2 ) THEN
140C SET BOUNDARY FOR GAMMA SUM
141 GABOU = SEUGF
142 NNTOT = INT(FNEUT2)
143C CALCULATE BOUNDARY NNTOT OF PARTICLE LOOP RATHER AT RANDOM THAN BY
144C ROUNDING OF FNEUT2 TO AVOID DIGITIZING EFFECTS ON THE NEUTRAL
145C PARTICLE COMPOSITION AT COLLISIONS WITH LOW MULTIPLICITY
146 IF ( NNTOT+RD(NEUTOT+2) .GE. FNEUT2 ) NNTOT = NNTOT+1
147 ELSE
148 IF ( RC3TO2 .LE. 0.D0 ) GOTO 1010
149 GABOU = GABOU + SEUGF* RC3TO2
150 NNTOT = INT(FNEUT)
151 IF ( NNTOT+RD(NEUTOT+3) .GE. FNEUT ) NNTOT = NNTOT+1
152 ENDIF
153 IF ( DEBUG ) WRITE(MDEBUG,*) ' J,NNTOT=',J,NNTOT
154C START NEUTRAL PARTICLE PRODUCTION LOOP
155 1003 CONTINUE
156 IF ( K .LT. NNTOT ) THEN
157 RNDM = RD(K)
158 ELSEIF ( K .EQ. NNTOT ) THEN
159C RENORMALIZE THE RANDOM NUMBER, THAT ONLY PI(0) OR ETA IS PRODUCED
160C BUT PAIR PRODUCTION BECOMES IMPOSSIBLE
161 RNDM = RD(K) * RPIER
162 ELSEIF ( K .GT. NNTOT ) THEN
163 GOTO 1010
164 ENDIF
165 IF ( RNDM .LE. RPI0R ) THEN
166C PI(0)
167 SGAMMA = SGAMMA + 2.D0
168 NPIZER(J) = NPIZER(J) + 1
169 K = K + 1
170
171 ELSEIF ( RNDM .LE. RPIER ) THEN
172C ETA
173 CALL RMMAR( RDETA,1,1 )
174 IF ( RDETA .LE. CETA(1) ) THEN
175 SGAMMA = SGAMMA + 2.D0
176 NETA(J,1) = NETA(J,1) + 1
177 ELSEIF ( RDETA .LE. CETA(2) ) THEN
178 SGAMMA = SGAMMA + 6.D0
179 NETA(J,2) = NETA(J,2) + 1
180 ELSEIF ( RDETA .LE. CETA(3) ) THEN
181 SGAMMA = SGAMMA + 2.D0
182 NETA(J,3) = NETA(J,3) + 1
183 ELSE
184 SGAMMA = SGAMMA + 1.D0
185 NETA(J,4) = NETA(J,4) + 1
186 ENDIF
187 K = K + 1
188
189 ELSEIF ( RNDM .LE. RPEKR ) THEN
190C K0S/K0L PAIR; RPEKR IS NORMALIZED FOR K0 PAIR FORMATION
191C THE UA5 GAMMA YIELD DOES NOT INCLUDE GAMMAS FROM K DECAY !!!
192C SEE: ANSORGE ET AL., Z. PHYS. C43 (1989) 75
193 NKA0(J) = NKA0(J) + 1
194 K = K + 2
195 ELSEIF ( RNDM .LE. RPEKNR ) THEN
196C NEUTRON-ANTINEUTRON PAIR
197 NNUCN(J) = NNUCN(J) + 1
198 K = K + 2
199 ELSE
200C HYPERON-ANTIHYPERON PAIR
201C AVERAGE NEUTRAL HYPERON PAIR L0 --> .357*2 GAMMAS = 0.714 GAMMAS
202C S0 --> L0 + 1 GAMMA = 1.714 GAMMAS
203C THEY ARE INCLUDED IN UA5 GAMMA MULTIPLICITIES, THEREFORE COUNT
204 SGAMMA = SGAMMA + 2.428D0
205 NHYPN(J) = NHYPN(J) + 1
206 K = K + 2
207 ENDIF
208 GOTO 1003
209 1010 CONTINUE
210 IF ( DEBUG ) WRITE(MDEBUG,1020) ( 2*NNUCN(J),2*NKA0(J),
211 * 2*NHYPN(J),NETA(J,1),NETA(J,2),NETA(J,3),NETA(J,4),
212 * NPIZER(J),J=2,3 ), NNTOT,GABOU,SGAMMA,SGAMMA/GABOU
213 1020 FORMAT(' PARNUM: NEUTRALS (1.,2.STRING)=',8I5,/
214 * ' NEUTRALS (3. STRING) =',8I5,/
215 * ' NNTOT,SEUGF2+3,SGAMMA,RATIO=',I6,3(2X,F10.5))
216C REJECT ALL NEUTRALS, IF SUM OF GAMMAS DEVIATES BY MORE THAN SIGMA
217 IF ( (SGAMMA - GABOU)**2 .GT. GABOU ) THEN
218 NCOUNT = NCOUNT + 1
219C AFTER 20 TRIES SET FLAG INUMFL TO 1 AND RETURN
220 IF ( NCOUNT .LE. 20 ) GOTO 1002
221 INUMFL = 1
222 RETURN
223 ENDIF
224C ALL NEUTRALS
225 NNN = NNUCN(2) + NNUCN(3)
226 NKN = NKA0(2) + NKA0(3)
227 NHN = NHYPN(2) + NHYPN(3)
228 NETAS(2) = NETA(2,1) + NETA(2,2) + NETA(2,3) + NETA(2,4)
229 NETAS(3) = NETA(3,1) + NETA(3,2) + NETA(3,3) + NETA(3,4)
230 NET = NETAS(2) + NETAS(3)
231 NPN = NPIZER(2) + NPIZER(3)
232
233C ...FOR CHARGED
234 I = 1
235 CALL RMMAR( RD,NCH-1,1 )
236C START CHARGED PARTICLE PRODUCTION LOOP
237 1101 CONTINUE
238 RNDM = RD(I)
239 IF ( RNDM .LT. PPICH ) THEN
240C PI(+-)
241 NPC = NPC + 1
242 I = I + 1
243 ELSEIF ( RNDM .LT. PPINCH ) THEN
244C PROTON/ANTIPROTON PAIR
245 NNC = NNC + 1
246 I = I + 2
247 ELSEIF ( RNDM .LT. PPNKCH ) THEN
248C KAON(+,-) PAIR
249 NKC = NKC + 1
250 I = I + 2
251 ELSE
252C CHARGED HYPERON/ANTIHYPERON PAIR
253 NHC = NHC + 1
254 I = I + 2
255 ENDIF
256 IF ( I .LT. NCH ) THEN
257 GOTO 1101
258 ELSEIF ( I .EQ. NCH ) THEN
259C ONLY 1 CHARGED PARTICLE TO BE PRODUCED WHICH IS PI(+-)
260 NPC = NPC + 1
261 ENDIF
262C CORRECT CHARGED PION NUMBER FOR DECAY OF ETA'S
263 NCORR = 2 * ( NETA(2,3) + NETA(2,4) + NETA(3,3) + NETA(3,4) )
264 NPC = MAX( 0, NPC - NCORR )
265 IF ( DEBUG ) WRITE(MDEBUG,*) ' NPC,NPN,NCORR,LASTPI=',
266 * NPC,NPN,NCORR,LASTPI
267 ENDIF
268C CORRECT NUMBER OF CHARGED AND NEUTRAL PIONS FOR RESONANCE DECAY
269C (NRESPC, NRESPN)
270 NPC = MAX( 0, NPC - NRESPC + LASTPI )
271C INCREASE NPN ADDITIONALLY BY 1 TO MEET UA5 DATA, WHICH REPRODUCE ON
272C AVERAGE ONE EXCHANGED CHARGE (LASTPI = +1).
273 NPN = MAX( 0, NPN - NRESPN - LASTPI + 1 )
274C TOTAL NUMBER OF CHARGED PARTICLES
275 NCH = (NNC + NKC + NHC) * 2 + NPC
276C NOW ALL PARTICLES ARE DETERMINED
277 IF ( DEBUG ) WRITE(MDEBUG,*)
278 * 'PARNUM: TOT.CHARGED=',2*NNC,2*NKC,2*NHC,NPC,
279 * 'PARNUM: TOT.NEUTRAL=',2*NNN,2*NKN,2*NHN,NET,NPN
280
281 RETURN
282 END
Note: See TracBrowser for help on using the repository browser.