source: trunk/MagicSoft/Simulation/Corsika/Mmcs/hdpm.f@ 6495

Last change on this file since 6495 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: 32.5 KB
Line 
1 SUBROUTINE HDPM
2
3C-----------------------------------------------------------------------
4C H(ADRONIC) D(UAL) P(ARTON) M(ODEL)
5C
6C GENERATOR OF HADRONIC COLLISION INSPIRED BY DUAL PARTON MODEL
7C THIS SUBROUTINE IS CALLED FROM SDPM
8C-----------------------------------------------------------------------
9
10 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
11*KEEP,CONST.
12 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
13 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
14*KEEP,DPMFLG.
15 COMMON /DPMFLG/ NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM
16 INTEGER NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM
17*KEEP,ELADPM.
18 COMMON /ELADPM/ ELMEAN,ELMEAA,IELDPM,IELDPA
19 DOUBLE PRECISION ELMEAN(37),ELMEAA(37)
20 INTEGER IELDPM(37,13),IELDPA(37,13)
21*KEEP,ELASTY.
22 COMMON /ELASTY/ ELAST,IELIS,IELHM,IELNU,IELPI
23 DOUBLE PRECISION ELAST
24 INTEGER IELIS(20),IELHM(20),IELNU(20),IELPI(20)
25*KEEP,INDICE.
26 COMMON /INDICE/ NNUCN,NKA0,NHYPN,NETA,NETAS,NPIZER,
27 * NNC,NKC,NHC,NPC,NCH,NNN,NKN,NHN,NET,NPN
28 INTEGER NNUCN(2:3),NKA0(2:3),NHYPN(2:3),NETA(2:3,1:4),
29 * NETAS(2:3),NPIZER(2:3),
30 * NNC,NKC,NHC,NPC,NCH,NNN,NKN,NHN,NET,NPN
31*KEEP,INTER.
32 COMMON /INTER/ 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 * IDIF,ITAR
37 DOUBLE PRECISION AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB,
38 * FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3,
39 * RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG,
40 * WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN
41 INTEGER IDIF,ITAR
42*KEEP,ISTA.
43 COMMON /ISTA/ IFINET,IFINNU,IFINKA,IFINPI,IFINHY
44 INTEGER IFINET,IFINNU,IFINKA,IFINPI,IFINHY
45*KEEP,LEPAR.
46 COMMON /LEPAR/ LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS
47 INTEGER LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS
48*KEEP,MULT.
49 COMMON /MULT/ EKINL,MSMM,MULTMA,MULTOT
50 DOUBLE PRECISION EKINL
51 INTEGER MSMM,MULTMA(37,13),MULTOT(37,13)
52*KEEP,NEWPAR.
53 COMMON /NEWPAR/ EA,PT2,PX,PY,TMAS,YR,ITYP,
54 * IA1,IA2,IB1,IB2,IC1,IC2,ID1,ID2,IE1,IE2,IF1,IF2,
55 * IG1,IG2,IH1,IH2,II1,II2,IJ1,NTOT
56 DOUBLE PRECISION EA(3000),PT2(3000),PX(3000),PY(3000),TMAS(3000),
57 * YR(3000)
58 INTEGER ITYP(3000),
59 * IA1,IA2,IB1,IB2,IC1,IC2,ID1,ID2,IE1,IE2,IF1,IF2,
60 * IG1,IG2,IH1,IH2,II1,II2,IJ1,NTOT
61*KEEP,PAM.
62 COMMON /PAM/ PAMA,SIGNUM
63 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
64*KEEP,PARPAR.
65 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
66 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
67 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
68 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
69 INTEGER ITYPE,LEVL
70*KEEP,PARPAE.
71 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
72 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
73 * (CURPAR(4), PHI ), (CURPAR(5), H ),
74 * (CURPAR(6), T ), (CURPAR(7), X ),
75 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
76 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
77 * (CURPAR(12),ECM )
78*KEEP,RANDPA.
79 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
80 DOUBLE PRECISION FAC,U1,U2
81 REAL RD(3000)
82 INTEGER ISEED(103,10),NSEQ
83 LOGICAL KNOR
84*KEEP,RATIOS.
85 COMMON /RATIOS/ RPI0R,RPIER,RPEKR,RPEKNR,PPICH,PPINCH,PPNKCH,
86 * ISEL,NEUTOT,NTOTEM
87 DOUBLE PRECISION RPI0R,RPIER,RPEKR,RPEKNR,PPICH,PPINCH,PPNKCH
88 INTEGER ISEL,NEUTOT,NTOTEM
89*KEEP,RESON.
90 COMMON /RESON/ RDRES,RESRAN,IRESPAR
91 REAL RDRES(2),RESRAN(1000)
92 INTEGER IRESPAR
93
94*KEEP,REST.
95 COMMON /REST/ CONTNE,TAR,LT
96 DOUBLE PRECISION CONTNE(3),TAR
97 INTEGER LT
98*KEEP,RUNPAR.
99 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
100 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
101 * MONIOU,MDEBUG,NUCNUC,
102 * CETAPE,
103 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
104 * N1STTR,MDBASE,
105 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
106 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
107 * ,GHEISH,GHESIG
108 COMMON /RUNPAC/ DSN,HOST,USER
109 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
110 REAL STEPFC
111 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
112 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
113 * N1STTR,MDBASE
114 INTEGER CETAPE
115 CHARACTER*79 DSN
116 CHARACTER*20 HOST,USER
117
118 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
119 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
120 * ,GHEISH,GHESIG
121*KEND.
122
123C-----------------------------------------------------------------------
124
125 IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9)
126 444 FORMAT(' HDPM : CURPAR=',1P,9E10.3)
127
128C SET ANTI-LEADER TO PROTON OR NEUTRON; TARGET IS ALWAYS NUCLEON
129 CALL RMMAR( RD,1,1 )
130 IF ( RD(1) .LE. CONTNE(LT) ) THEN
131 ITAR = 13
132 ELSE
133 ITAR = 14
134 ENDIF
135
136
137C CALCULATE LAB AND CM ENERGY
138 IF ( ITYPE .NE. 1 ) THEN
139 ELAB = PAMA(ITYPE) * GAMMA
140 PLAB = ELAB * BETA
141 S = PAMA(ITYPE)**2 + PAMA(ITAR)**2 + 2.D0*PAMA(ITAR)*ELAB
142 ELSE
143C FOR GAMMA-INDUCED REACTION TAKE PI(0) AS LEADING PARTICLE
144 ITYPE = 7
145 ELAB = GAMMA
146 PLAB = ELAB
147 S = PAMA(ITAR)**2 + 2.D0*PAMA(ITAR)*ELAB
148 ENDIF
149
150 ECMDPM = SQRT(S)
151 IF ( DEBUG ) WRITE(MDEBUG,*) 'HDPM : ITYPE,ELAB,PLAB,S,ECMDPM=',
152 * ITYPE,SNGL(ELAB),SNGL(PLAB),SNGL(S),SNGL(ECMDPM)
153
154C LN(S), LN(S)**2 AND RAPIDITY OF C. M. SYSTEM IN LAB
155 SLOG = LOG(S)
156 SLOGSQ = SLOG**2
157 SMLOG = LOG( 2.D0 * PAMA(ITAR) * ELAB )
158 ELABLG = LOG(ELAB)
159 EPLUSP = ELAB + PLAB
160* YCM = 0.5D0 * LOG( (ELAB+PAMA(ITAR)+PLAB)/(ELAB+PAMA(ITAR)-PLAB) )
161 YCM = 0.5D0 * LOG( (EPLUSP**2 +PAMA(ITAR)*EPLUSP)/
162 * (PAMA(ITYPE)**2+PAMA(ITAR)*EPLUSP) )
163 IF ( DEBUG ) WRITE(MDEBUG,*)'HDPM : SLOG,SLOGSQ,YCM=',
164 * SNGL(SLOG),SNGL(SLOGSQ),SNGL(YCM)
165
166C-----------------------------------------------------------------------
167C RETURN POINT IF CALCULATION OF PARTICLES GOES WRONG
168 1 CONTINUE
169
170 IF ( ITYPE .NE. 7 ) THEN
171C CHOOSE NUMBER OF INTERACTIONS IN TARGET
172 CALL TARINT
173 ELSE
174C FOR GAMMA-INDUCED REACTIONS TAKE ALWAYS ONE COLLISION
175 GNU = 1.D0
176 ENDIF
177
178C-----------------------------------------------------------------------
179C NO DIFFRACTION IF
180C OR THE NUMBER OF INTERACTIONS IN TARGET IS CHOSEN RANDOMLY
181C AND MORE THAN ONE INTERACTION TAKES PLACE
182C OR PRIMARY PARTICLE IS GAMMA (PI0)
183C NOW NFLDIF DECIDES WHETHER DIFFRACTIVE PROCESS POSSIBLE OR NOT
184 IF ( ( NFLAIN.EQ.0 .AND. GNU.GT.1.D0 .AND. NFLDIF.EQ.0 )
185 * .OR. ( ITYPE .EQ. 7 ) ) THEN
186 IDIF = 0
187 ELSE
188C SET DIFFRACTION FLAG IF RANDOM NUMBER < PROBABILITY
189 CALL RMMAR( RD,1,1 )
190C IDIF IS 0 : NO DIFFRACTION ; IDIF IS 1 : DIFFRACTION
191C DIFFRACTION RISES WITH ENERGY AND SATURATES AT 10000 GEV
192C ### DAS TUT ES ABER NICHT: ES IST KONSTANT 0.15 (SIEHE DPFUNC) !!!!
193 IF ( RD(1) .GT. DPFUNC(ECMDPM) ) THEN
194 IDIF = 0
195 ELSE
196 IDIF = 1
197 ENDIF
198 ENDIF
199
200
201C PRINTOUT FOR DEBUG
202 IF ( DEBUG ) THEN
203 WRITE(MDEBUG,*) ' DIFFRACTIVE INTERACTION (0/1) = ',IDIF
204 ENDIF
205
206C SET COUNTER FOR REPEAT TO 0
207 NREPRD = 0
208
209C GENERATION OF INTERACTION
210 1919 CONTINUE
211
212C FLAG TO CHECK NUMBER OF SECONDARIES;
213C IS CHANGED TO 1 IF SECONDARY MULTIPLICITY IS LOW
214 ISEL = 0
215C SET LEADING PARTICLE TO INCOMING PARTICLE AND ANTI-LEADER TO NUCLEON
216C (AS IT COMES FROM TARGET NUCLEUS) BOTH MAY BE CHANGED BY LEPACX
217 LEPAR1 = ITYPE
218 LEPAR2 = ITAR
219
220 IF ( IDIF .EQ. 0 ) THEN
221C-----------------------------------------------------------------------
222C NON SINGLE DIFFRACTIVE PROCESS STARTS HERE
223
224 CALL NSD
225C CHARGE EXCHANGE ENABLED? EXCHANGE LEADER AND ANTI-LEADER
226 LASTPI = 0
227 NRESPC = 0
228 NRESPN = 0
229 NCPLUS = 0
230 IF ( NFLCHE .EQ. 0 ) THEN
231 CALL LEPACX( ECMDPM,ELABLG,LEPAR1,1 )
232 CALL LEPACX( ECMDPM,ELABLG,LEPAR2,2 )
233 ENDIF
234 1921 CONTINUE
235 CALL RNEGBI( NCH,AVCH,ECMDPM )
236C NCH IS # OF ALL CHARGED PARTICLES INCLUDING EXCESS FROM TARGET
237 IF ( NCH .LT. 1 ) THEN
238 IF ( LEPAR1 .LT. 50 .OR. LEPAR2 .LT. 50 ) THEN
239 NREPRD = NREPRD + 1
240 IF ( NREPRD .GT. 10 ) GOTO 1
241 GOTO 1921
242 ELSE
243C INTERACTION IS ONLY RESONANCE PRODUCTION
244 ISEL = 1
245 ENDIF
246 ENDIF
247C WIDTH PLATEAU FOR CLUSTERS AND FOR CALCULATION OF CENTR.RAP.DENSITY
248 DELRAP = 0.6722D0 * (2.95D0 + 0.0302D0 * SLOG)
249C SET RSLOG FOR CALCULATION OF PARTICLE RATIOS
250 RSLOG = SLOG
251C AVERAGE TRANSVERSE MOMENTUM
252 CALL AVEPT( ECMDPM,SLOG )
253
254 ELSE
255C-----------------------------------------------------------------------
256C SINGLE DIFFRACTIVE PROCESS STARTS HERE
257
258 1920 CONTINUE
259 CALL DIFRAC( NRETDF )
260 IF ( NRETDF .EQ. 1 ) GOTO 1
261C CHARGE EXCHANGE ENABLED? EXCHANGE CHARGE OF DIFFRACTING PARTICLE
262 LASTPI = 0
263 NRESPC = 0
264 NRESPN = 0
265 NCPLUS = 0
266 IF ( NFLCHE .EQ. 0 ) THEN
267 IF ( YY0 .GT. 0.D0 ) THEN
268C PROJECTILE DIFFRACTION
269 CALL LEPACX( ECMDIF,DMLOG,LEPAR1,1 )
270 ELSE
271C TARGET DIFFRACTION
272 CALL LEPACX( ECMDIF,DMLOG,LEPAR2,2 )
273 ENDIF
274 ENDIF
275C FLUCTUATION OF MULTIPLICITY ACCORDING TO NEG.BIN. DISTRIBUTION
276 CALL RNEGBI( NCH,AVCH,ECMDIF )
277C REPEAT CALCULATION AS SOMETHING WENT WRONG
278 IF ( NCH .LT. 1 ) THEN
279 IF ( (YY0 .GT. 0.D0 .AND. LEPAR1 .LT. 50) .OR.
280 * (YY0 .LT. 0.D0 .AND. LEPAR2 .LT. 50) ) THEN
281 NREPRD = NREPRD + 1
282 IF ( NREPRD .GT. 10 ) GOTO 1
283 GOTO 1920
284 ELSE
285C DIFFRACTIVE INTERACTION IS ONLY RESONANCE PRODUCTION
286 ISEL = 1
287 ENDIF
288 ENDIF
289C SET RSLOG FOR CALCULATION OF PARTICLE RATIOS
290 RSLOG = DLOG
291C HERE WE USE ECMDPM, BECAUSE THE MOMENTUM TRANSFER IS DEPENDENT
292C ON THE ENERGY OF THE TOTAL SYSTEM AND NOT ON THE DIFFRACTING MASS
293 CALL AVEPT( ECMDPM,SLOG )
294
295 ENDIF
296
297C-----------------------------------------------------------------------
298C NOW FOR NSD AND DIFFRACTIVE PROCESSES
299
300C IN CASE OF LOW MULTIPLICITY SET FLAG ISEL
301 IF ( NCH .LE. 2 ) ISEL=1
302C FNCH IS FLUCTUATING TOT.NUMBER OF CHARGED PARTICLES FOR ALL 3 STRINGS
303 FNCH = DBLE(NCH)
304C RATIO ALL CHARGED PARTICLES WITH FLUCTUATION/WITHOUT FLUCTUATION
305 XZ = FNCH / AVCH
306C FNCH3 IS FLUCTUATING NUMBER OF CHARGED PARTICLES FOR 3RD STRING
307 FNCH3 = XZ * AVCH3
308C FNCH2 IS FLUCTUATING NUMBER OF CHARGED PARTICLES 1ST AND 2ND STRING
309 FNCH2 = FNCH - FNCH3
310C RC3TO2 IS RATIO (CHARGED 3RD STRING)/(CHARGED 1ST AND 2ND STRING)
311 IF ( FNCH2 .NE. 0.D0 ) THEN
312 RC3TO2 = FNCH3 / FNCH2
313 ELSE
314 RC3TO2 = 0.D0
315 ENDIF
316 IF ( DEBUG ) WRITE(MDEBUG,*) ' FNCH,FNCH2,FNCH3,RC3TO2=',
317 * SNGL(FNCH),SNGL(FNCH2),SNGL(FNCH3),SNGL(RC3TO2)
318
319C IS NUMBER OF NEUTRALS FLUCTUATING AS NUMBER OF CHARGED ?
320 IF ( NFLPIF .EQ. 0 .OR. IDIF .EQ. 1 .OR. ECMDPM .LT. 60.D0 ) THEN
321C SET NUMBER OF GAMMAS ACCORDING TO NEG. BIN. VARIABLE XZ
322C AS NUMBER OF NEUTRALS FLUCTUATES AS CHARGED.
323 SEUGF = SEUGP * XZ
324 ZG = XZ
325 ELSE
326C NFLPIF IS 1 MEANS: # OF PI(0) FLUCTUATES AS MEASURED AT COLLIDER
327 IF ( ECMDPM .LT. 200.D0 ) THEN
328 SEUGF = SEUGP * XZ
329* SEUGF = (0.0786D0*SLOG-0.010D0)*FNCH2 + (0.391D0*SLOG+0.305D0)
330 ELSE
331C DETERMINE NEW NUMBER OF GAMMAS WITH FLUCTUATION AROUND SEUGP*XZ
332 AGR = EXP(-XZ)
333 DGR = SEUGP * XZ * (0.9823D0 - 0.3756D0 * AGR)
334 SGS = DGR * (0.14718D0 + 2.53213D0 * AGR)
335 723 CONTINUE
336 SEUGF = 0.88D0 * RANNOR(DGR,SGS)
337 IF ( SEUGF .LT. 1.D0 ) GOTO 723
338 ENDIF
339C SET NEGATIVE BINOMIAL VARIABLE ZG FOR GAMMAS
340 ZG = SEUGF / SEUGP
341 ENDIF
342 SEUGF = MAX( 1.D0, SEUGF )
343 IF ( DEBUG ) WRITE(MDEBUG,*) 'HDPM :XZ,ZG,SEUGF=',
344 * SNGL(XZ),SNGL(ZG),SNGL(SEUGF)
345
346C-----------------------------------------------------------------------
347C RATIO ALL-NUCLEON/ALL-CHARGED
348C PARAMETRISATION FROM UA5, NUCL. PHYS. B291 (1987) 445 EQ.(2.4)
349 RNUCCH = MAX( 0.D0, -0.008D0 + 0.00865D0 * RSLOG )
350C NUMBER FOR DIRECT NEUTRON/ANTINEUTRON PRODUCTION 1ST AND 2ND STRING
351C MULTIPLY BY 0.5 BECAUSE RATIO RNUCCH GIVES (ALL-NUCL)/(ALL-CHARGED)
352C AND HERE ONLY THE NEUTRON-ANTINEUTRONS ARE COUNTED
353 FNUCN = 0.5D0 * RNUCCH * FNCH2
354C RATIO (ALL CHARGED SIGMAS)/(ALL CHARGED) IS 1/3 OF ALL STRANGE BARYON
355C PARAMETRISATION FORM UA5, NUCL. PHYS. B291 (1987) 445 EQ.(2.5)
356 RHYPCH = MAX( 0.D0, (-0.007D0 + 0.0028D0 * RSLOG) * OB3 )
357C NEUTRAL STRANGE BARYONS ARE DOUBLE OF CHARGED STRANGE BARYONS
358 FHYPN = 2.D0 * RHYPCH * FNCH2
359C CORRECT NUMBER OF GAMMAS FROM NEUTRAL HYPERON DECAY S0-->L+GAMMA
360 SEUGFC = MAX( 0.D0, SEUGF - 0.5D0 * FHYPN )
361C RATIO CHARGED-KAON/CHARGED PIONS
362C PARAMETRISATION FROM UA5, NUCL. PHYS. B291 (1987) 445 EQ.(2.7)
363 RKPI = MAX (0.D0, 0.024D0 + 0.0062D0 * RSLOG )
364C RKCH IS RATIO (CHARGED-KAON)/(ALL-CHARGED) DERIVED FROM RKPI;
365C THE FACTOR 0.5 IN FRONT OF RNUCCH IS BECAUSE ONLY HALF OF NUCLEONS
366C ARE P/PBAR. THE 1.17 IS AN APROXIMATE CONVERSION FACTOR FROM
367C (ALL-NUCL)/(ALL-CHARGED) TO (ALL-NUCL)/(CHARGED-PI), WHICH IS A BIT
368C ENERGY DEPENDENT (1.14 ...1.21) SEE GEICH-GIMBEL TABLE 7.1
369 RKCH = RKPI / (1.D0 + RKPI + (0.5D0*RNUCCH+RHYPCH) * 1.17D0)
370C K0/K0-BAR FOR 1ST AND 2ND STRING
371C NEUTRAL KAONS ARE PRODUCED WITH THE SAME RATE AS CHARGED KAONS
372 FKA0 = RKCH * FNCH2
373C RATIO ETA/PI(0) IS ASSUMED TO BE INDEPENDENT OF ENERGY = 0.19
374C SEE: ANSORGE ET AL. (UA5-COLLABORATION) Z.PHYS.C43(1989)75
375* RETPI0 = 0.19D0
376C RATIO ETA/PI(0) IS ASSUMED TO BE DEPENDENT ON ENERGY
377C SEE: GEICH-GIMBEL,INT.J.MOD.PHYS.A4(1989)1527 TAB.7.1
378C HECK'S FIT: RETPI0 IS 0.06 + 0.006*RSLOG + 0.0011*RSLOG**2
379 RETPI0 = 0.06D0 + 0.006D0 * RSLOG + 0.0011D0 * RSLOG**2
380C AUXIL1 IS FRACTION OF PI(0)/(PI(0)+ETA)
381 AUXIL1 = 1.D0 / (1.D0 + RETPI0)
382C NUMBER OF GAMMAS FROM PI(0) IS 2, FROM ETA IS 3.216 IN AVERAGE;
383C AUXIL2 IS NUMBER OF GAMMA-PRODUCING-PARTICLES: PI(0) AND ETA
384 AUXIL2 = SEUGFC / ( AUXIL1 * 2.D0 + (1.D0 - AUXIL1) * 3.216D0 )
385 FETA = (1.D0 - AUXIL1) * AUXIL2
386 FPI0 = AUXIL1 * AUXIL2
387C CORRECT FPI0 BY DECAYS OF STRANGE BARYONS; NEUTRAL: FHYPN*0.357
388C CHARGED: 0.5*FNCH2*RHYPCH*0.5157; IT YIELDS FHYPN*(0.357+0.12893)
389 FPI0 = MAX( 0.D0, FPI0 - FHYPN * 0.486D0 )
390C SUMMED NEUTRAL PARTICLES FOR 1ST AND 2ND STRING
391 FNEUT2 = FNUCN + FKA0 + FHYPN + FETA + FPI0
392C NEUTRAL PARTICLES FROM 3RD STRING
393 FNEUT3 = RC3TO2 * FNEUT2
394C TOTAL NUMBER OF NEUTRALS
395 FNEUT = FNEUT2 + FNEUT3
396 NEUTOT = NINT( FNEUT )
397C CALCULATE TOTAL NUMBER OF PARTICLES TO BE CREATED
398 NTOTEM = NCH + NEUTOT
399 IF ( DEBUG ) WRITE(MDEBUG,*)
400 * ' FNUCN,FKA0,FHYPN,FETA,FPI0,FNEUT2,FNEUT3,NTOTEM=',
401 * SNGL(FNUCN),SNGL(FKA0),SNGL(FHYPN),SNGL(FETA),SNGL(FPI0),
402 * SNGL(FNEUT2),SNGL(FNEUT3),NTOTEM
403C LIMIT OF SECONDARIES PRODUCED (GIVEN BY SIZE OF ARRAY : 3000)
404C LIMIT IS ARRAY SIZE - SIZE OF LARGEST TARGET NUCLEUS(40)
405 IF ( NTOTEM .GE. 2956 ) THEN
406 WRITE(MONIOU,*) 'HDPM : REJECT EVENT WITH ',NTOTEM,
407 * ' SECONDARIES'
408 GOTO 1
409 ENDIF
410C SPECIAL TREATMENT IF MULTIPLICITY IS TOO LOW
411 IF ( NTOTEM .LE. 3 ) ISEL = 1
412
413C FRACTION OF THE VARIOUS NEUTRAL PARTICLES (NN, K(0), L+S0 AS PAIRS)
414C NORMALIZE WITH THE SUM OF ALL NEUTRAL PARTICLES
415 FNORML = 1.D0 / ( 0.5D0 * (FNUCN+FKA0+FHYPN) + FETA + FPI0 )
416 RNUCNR = FNUCN * FNORML * 0.5D0
417 RKA0R = FKA0 * FNORML * 0.5D0
418 RHYPNR = FHYPN * FNORML * 0.5D0
419 RETAR = FETA * FNORML
420 RPI0R = FPI0 * FNORML
421C CUMULATED RATIOS (NN, K(0), LAMBDA+SIGMA0 AS PAIRS)
422 RPIER = RPI0R + RETAR
423 RPEKR = RPIER + RKA0R
424 RPEKNR = RPEKR + RNUCNR
425C THEN THE REMAINDER (1-RPEKNR) MUST BE NEUTRAL HYPERON PAIRS
426 IF ( DEBUG ) WRITE(MDEBUG,*)
427 * ' RPI0R,RETAR,RKA0R,RNUCNR,RHYPNR,FNORML=',
428 * SNGL(RPI0R),SNGL(RETAR),SNGL(RKA0R),SNGL(RNUCNR),SNGL(RHYPNR),
429 * SNGL(FNORML)
430
431C PROBABILITY TO PRODUCE CHARGED PIONS IS PROBABILITY NOT TO PRODUCE
432C CHARGED KAONS OR PROTONS OR CHARGED HYPERONS, WHERE PROTON/ANTIPROTON
433C IS HALF OF (ALL-NUCL)/(ALL-CHARGED)
434 AUXIL = RKCH + 0.5D0 * RNUCCH + RHYPCH
435 AUXIL3 = 1.D0 - AUXIL
436C RENORMALIZATION AS P/P_BAR, K+-, AND HYPERONS ARE PRODUCED IN PAIRS
437C AUXIL2 IS INVERSE OF NORMALISATION
438 AUXIL2 = 1.D0 / (1.D0 - 0.5D0 * AUXIL)
439C CUMULATED PROBABILITIES (PP, K+-, SIGMA+- AS PAIRS)
440 PPICH = AUXIL3 * AUXIL2
441 PPINCH = PPICH + 0.25D0 * RNUCCH * AUXIL2
442 PPNKCH = PPINCH + 0.5D0 * RKCH * AUXIL2
443C THEN THE REMAINDER (1-PPNKCH) MUST BE CHARGED HYPERON PAIRS
444 IF ( DEBUG ) WRITE(MDEBUG,*) ' PPICH,PPINCH,PPNKCH=',
445 * SNGL(PPICH),SNGL(PPINCH),SNGL(PPNKCH)
446
447C NOW SELECT HOW MANY PARTICLES OF EACH TYPE ARE PRODUCED
448 CALL PARNUM( INUMFL )
449 IF ( INUMFL .NE. 0 ) GOTO 1919
450
451C DEFINE PARTICLE NUMBERS WHERE SPECIAL RAPIDITY IS CALCULATED
452C FOR PARTICLES FROM TARGET (THIRD STRING)
453 PPP = RC3TO2 / (1.D0+RC3TO2)
454C NUMBER OF PARTICLES IN PROTON ANTIPROTON PAIRS FROM TARGET
455 ITA = NINT(PPP * 2.D0 * NNC)
456C NUMBER OF PARTICLES IN K+ K- PAIRS FROM TARGET
457 ITB = NINT(PPP * 2.D0 * NKC)
458C NUMBER OF PARTICLES IN SIGMA+ SIGMA- PAIRS FROM TARGET
459 ITC = NINT(PPP * 2.D0 * NHC)
460C NUMBER OF PI+ PI- FROM TARGET
461 ITD = NINT(PPP * NPC )
462C CALCULATE BOUNDARIES
463 IA1 = 2
464 IA2 = IA1 + ITA
465 IB1 = IA1 + 2 * NNC
466 IB2 = IB1 + ITB
467 IC1 = IB1 + 2 * NKC
468 IC2 = IC1 + ITC
469 ID1 = IC1 + 2 * NHC
470 ID2 = ID1 + ITD
471 IE1 = ID1 + NPC
472C NUMBER OF PARTICLES IN NEUTRON ANTINEUTRON PAIRS FROM TARGET
473 IE2 = IE1 + 2 * NNUCN(3)
474 IF1 = IE1 + 2 * NNN
475C NUMBER OF PARTICLES IN K0S K0L PAIRS FROM TARGET
476 IF2 = IF1 + 2 * NKA0(3)
477 IG1 = IF1 + 2 * NKN
478C NUMBER OF PARTICLES IN NEUTRAL HYPERON PAIRS FROM TARGET
479 IG2 = IG1 + 2 * NHYPN(3)
480 IH1 = IG1 + 2 * NHN
481C NUMBER OF ETA FROM TARGET
482 IH2 = IH1 + NETAS(3)
483 II1 = IH1 + NET
484C NUMBER OF PI(0) FROM TARGET
485 II2 = II1 + NPIZER(3)
486 IJ1 = II1 + NPN
487 IF ( DEBUG ) THEN
488 WRITE(MDEBUG,*) ' CHARGED FROM TARGET:',ITA,ITB,ITC,ITD
489 WRITE(MDEBUG,*) ' NEUTRAL FROM TARGET:',
490 * 2*NNUCN(3),2*NKA0(3),2*NHYPN(3),NETAS(3),NPIZER(3)
491 WRITE(MDEBUG,*) ' NTOTEM,IJ1=',NTOTEM,IJ1
492 ENDIF
493C REDEFINE TOTAL NUMBER OF SECONDARY PARTICLES : NTOTEM
494C BY CHARGE EXCHANGE AND RESONANCE FORMATION THIS NUMBER MAY BE ALTERED
495 NTOTEM = IJ1 - 2
496
497C-----------------------------------------------------------------------
498C RATIO OF RAPIDITY DENSITY TO MEAN PSEUDORAPIDITY IN CENTER
499C PARAMETRISATION SEE CAPDEVIELLE, J.PHYS.G:NUCL.PHYS.15(1989)909,EQ.6
500 IF ( XZ .LT. 1.5D0 ) THEN
501 RDS = (0.24396D0 + 0.70150424D0 * XZ)**2
502 ELSE
503 RDS = (0.55685D0 + 0.48664753D0 * XZ)**2
504 ENDIF
505C CALCULATE NOW: DN/DY AT Y = 0; DC0 IS AVERAGE PSEUDORAPIDITY DENSITY
506C TRAP IS RATIO (RAPID.DENS.)/(PSEUDORAP.DENS.) IN CENTER OF RAPIDITY
507 TRAP = 1.25D0
508 IF ( IDIF .EQ. 0 .AND. ECMDPM .GT. 19.4D0 )
509 * TRAP = MAX( 1.D0, 1.28852D0 - 0.0065D0 * SMLOG )
510 DCN2 = DC0 * RDS * TRAP
511 IF ( DEBUG ) WRITE(MDEBUG,*) ' RDS,TRAP,DCN2=',
512 * SNGL(RDS),SNGL(TRAP),SNGL(DCN2)
513C AMPLITUDE OF GAUSSIAN 1ST AND 2ND STRING
514 ATG2 = FNCH2 / (5.0132566D0 * WIDC2)
515C NEW DEFINITION OF POSITION BASED ON SEMI INCLUSIVE DATA
516 SQ2 = 2.D0 * ATG2 / DCN2
517C FINAL POSITION OF GAUSSIAN; WIDTH WIDC2 IS UNCHANGED
518 IF ( SQ2 .GT. 1.D0 ) POSC2 = WIDC2 * SQRT( 2.D0*LOG(SQ2) )
519C DENSITY OF CHARGED IN EXCESS FROM TARGET IN CENTER OF RAPIDITY
520 DCN3 = 0.5D0 * (GNU - 1.D0) * DCN2
521 IF (DEBUG) WRITE(MDEBUG,*) ' SQ2,POSC2,DCN3=',
522 * SNGL(SQ2),SNGL(POSC2),SNGL(DCN3)
523 IF ( DCN3 .GT. 0.D0 ) THEN
524C AMPLITUDE 3RD GAUSSIAN
525 ATG3 = FNCH3 / (5.0132566D0 * WIDC3)
526C AMPLITUDE IS DIVIDED BY DENSITY FOR GETTING CENTER OF 3RD GAUSSIAN
527 SQ3 = 2.D0 * ATG3 / DCN3
528C CHECK IF ADDITIVE MULTIPLICITY IS TOO LOW
529 IF ( SQ3 .GT. 1.D0 ) POSC3 = WIDC3 * SQRT( 2.D0*LOG(SQ3) )
530 IF (DEBUG) WRITE(MDEBUG,*)' SQ3,POSC3=',SNGL(SQ3),SNGL(POSC3)
531 ENDIF
532
533C NFLPI0 .EQ. 0 MEANS TREAT PI(0) RAPIDITY ACCORDING TO COLLIDER DATA
534 IF ( NFLPI0 .EQ. 0 ) THEN
535C RATIO OF RAPIDITY DENSITY TO MEAN PSEUDORAPIDITY AT CENTER WITH Z<1.5
536 IF ( ZG .LT. 1.5D0 ) THEN
537 RDG = (0.24396D0 + 0.70150424D0 * ZG)**2
538 ELSE
539 RDG = (0.55685D0 + 0.48664753D0 * ZG)**2
540 ENDIF
541C GAMMAS USE RATIO TRAG TO CALCULATE RATIO OF RAPIDITY TO
542C PSEUDO RAPIDITY DENSITY IN CENTER (TRAG = 1.1 * 0.5 ).
543C FACTOR 0.5 COMES FROM RATIO NEUTRAL/CHARGED, AS WE USE DC0, WHICH
544C IS AVERAGE PSEUDORAPIDITY DENSITY FOR CHARGED PIONS
545 TRAG = 0.55D0
546 IF ( IDIF .EQ. 0 ) THEN
547 IF ( ECMDPM .GT. 19.4D0 )
548 * TRAG = MAX( 0.4D0, 0.6658D0 - 0.01954D0 * SMLOG )
549 IF ( ECMDPM .LE. 50.D0 ) THEN
550 DCG = DC0 * RDG * TRAG
551 ELSEIF ( ECMDPM .LE. 200.D0 ) THEN
552 DCG = DC0 * RDG * TRAG * (1.D0 + 0.18D0 * LOG(ECMDPM/50.D0))
553 ELSE
554 DCG = DC0 * RDG * TRAG * 1.25D0
555 ENDIF
556 ELSE
557 DCG = DC0 * RDG * TRAG
558 ENDIF
559C DEFINE WIDTH OF STRINGS FOR NEUTRAL PIONS AND ETAS
560 WIDN2 = WIDC2 * MIN( 1.D0, 1.12275D0 - 0.0208D0 * RSLOG )
561C NEW DEFINITION OF CENTER OF GAUSSIAN BASED ON SEMI INCLUSIVE DATA
562C USING AMPLITUDE OF THE GAUSSIAN FOR NEUTRALS
563 AUXIL = 2.D0 / (5.0132566D0 * WIDN2 * DCG)
564C TOTAL MULTIPLICITY USED FOR 1ST AND 2ND STRING OF PI(0) AND ETA
565C IS GIVEN BY THEIR NUMBERS. ANALOGOUS FOR 3RD STRING
566 SP2 = DBLE ( NPIZER(2)+NETAS(2)) * AUXIL
567C FINAL CENTER OF GAUSSIANS FOR PI(0) AND ETA (WIDC2 IS UNCHANGED)
568 IF ( SP2 .GT. 1.D0 ) THEN
569 POSN2 = WIDN2 * SQRT( 2.D0 * LOG(SP2) )
570 ELSE
571 POSN2 = POSC2
572 ENDIF
573 WIDN3 = WIDN2
574 SP3 = DBLE(NPIZER(3)+NETAS(3)) * AUXIL
575 IF ( SP3 .GT. 1.D0 ) THEN
576 POSN3 = WIDN3 * SQRT( 2.D0 * LOG(SP3) )
577 ELSE
578 POSN3 = POSC3
579 ENDIF
580 ELSE
581C NFLPI0 .EQ. 1 MEANS RAPIDITY OF PI(0) AND ETA SAME AS THAT OF CHARGED
582 POSN2 = POSC2
583 WIDN2 = WIDC2
584 POSN3 = POSC3
585 WIDN3 = WIDC3
586 ENDIF
587 IF ( DEBUG ) WRITE(MDEBUG,*)
588 * ' ZG,RDG,DCG,SP2,SP3,POSN2,POSN3,WIDN2 =',
589 * SNGL(ZG),SNGL(RDG),SNGL(DCG),SNGL(SP2),SNGL(SP3),SNGL(POSN2),
590 * SNGL(POSN3),SNGL(WIDN2)
591
592C-----------------------------------------------------------------------
593 NREPR1 = 0
594C RETURN POINT. NUMBERS OF PARTICLES REMAIN UNCHANGED FOR NEXT TRY,
595C BUT INDIVIDUAL RAPIDITIES GET NEW VALUES.
596C START FROM BEGINNING IF NO MATCH AFTER 20 TRIES
597 30 CONTINUE
598 NREPR1 = NREPR1 + 1
599 IF ( NREPR1 .GT. 20 ) THEN
600 IF ( IDIF .EQ. 1 .AND. NREPRD .LE. 10 ) GOTO 1919
601 GOTO 1
602 ENDIF
603
604C FOR TOTAL NUMBER OF PARTICLES ADD 2 FOR LEADER AND ANTILEADER
605 NTOT = NTOTEM + 2
606
607C PRODUCTION OF INDIVIDUAL RAPIDITIES FOR ALL SECONDARY PARTICLES
608 CALL PARRAP
609CC IF ( DEBUG ) THEN
610CC WRITE (MDEBUG,*) ' RAPIDITIES:'
611CC WRITE (MDEBUG,134) (I,YR(I), I=3,NTOT)
612C134 FORMAT(' ',1P, (1X, I4, 5X, G13.6 ))
613CC ENDIF
614
615
616C CALCULATION OF CENTRAL RAPIDITY WITHOUT (ANTI)LEADER
617C MULTIPLICITY IN CENTER OF RAPIDITY DISTRIBUTION
618 IZN = 0.D0
619 IF ( IDIF .EQ. 0 ) THEN
620 DO 111 I = 3,NTOT
621 IF ( ABS(YR(I)) .LT. DELRAP ) IZN = IZN + 1
622 111 CONTINUE
623 IF ( IZN .LT. 1 ) THEN
624 IF ( ISEL .EQ. 0 ) GOTO 30
625C IN CASE OF FEW PARTICLES, SET PARTICLE NUMBER IN PLATEAU TO 1
626 IZN = 1
627 ENDIF
628C CENTRAL RAPIDITY DENSITY FOR CHARGED PARTICLES
629 IF ( NTOTEM .GE. 1 ) THEN
630 ZNC = MAX( 1.1D0, DBLE(NCH)*IZN/(DBLE(NTOTEM)*2.D0*DELRAP) )
631 ELSE
632 ZNC = 1.1D0
633 ENDIF
634 ELSE
635C DIFFRACTION: SHIFT RAPIDITIES + TAKE CENT.RAP.DENS. FROM PARAMETRISAT
636 DO 112 I = 3,NTOT
637 YR(I) = YR(I) + YY0
638 112 CONTINUE
639 ZNC = MAX( 1.1D0, DCN2 )
640 ENDIF
641
642C ZN ACCOUNTS FOR THE RISE OF PT WITH CENTRAL RAP.DENSITY. THE FORMULA
643C IS A FIT TO UA1 VALUES OF ARNISON ET AL, PHYS.LETT.B118(1982)167
644C REGARD, THAT OUR ZN IS DEFINED DIFFERENT FROM LITERATURE N BY 1
645C - - - - - -
646C MODIFICATION AFTER J.N. CAPDEVIELLE, (DEC.96)
647 IF ( ECMDPM .LE. 500.D0 ) THEN
648 ZN = MAX( 1.00001D0, 3.669D0 / ZNC**0.435D0 + 6.4D0 )
649 ELSE
650C TAKE INTO ACCOUNT THE RESULTS OF UA1/MIMI EXPERIMENT
651 IF ( ZNC .GE. 3.D0 ) THEN
652 ZN = 1.D0 /(ZNC*0.004111D0 + 0.145D0)+3.D0
653 ELSE
654C FOR ROCH < 3.00 (MIMI) (TO BE USED IN PTRAM)
655 ZM = 0.0033D0 * (ZNC-1.56D0)**2 + 0.406D0
656 ZN = 2.64D0/ZM + 3.D0
657 ENDIF
658 ENDIF
659C - - - - - -
660C NOW SET PARTICLE TYPE AND TRANSV. MOMENTA FOR NEW PARTICLES IN PPARAM
661C SET ALSO TRANSVERSE MASS FOR ALL PARTICLES (INCL. LEADER+ANTILEADER)
662 CALL PPARAM
663
664 IF ( IDIF .EQ. 0 ) THEN
665C NOW SET THE RAPIDITY OF THE ANTILEADER ACCORDING TO THE DISTRIBUTION
666C IN THE FEYNMAN X VARIABLE; SET THE RAPIDITY OF LEADER TO CONSUME
667C THE REMAINDER OF ENERGY
668 CALL LEDENY( LEDEFL )
669 IF ( LEDEFL .NE. 0 ) THEN
670 IF ( DEBUG ) WRITE(MDEBUG,*) ' LEDEFL=',LEDEFL
671 GOTO 30
672 ENDIF
673
674C CALCULATE FOR SINGLE COLLISION SYSTEM C.M. ENERGY + RAPIDITY SHIFT
675 IF ( GNU .LE. 1.D0 ) THEN
676 JGNU = 0.D0
677 DYGNU = 0.D0
678 ECMJAD = ECMDPM
679 ELSE
680C MULTIPLE COLLISION IN TARGET
681 JGNU = NINT(GNU-1.D0)
682C ADD ADDITIONALLY INTERACTING
683C TARGET NUCLEONS TO GET CORRECT JADACH FILTERING
684C CHOSE RANDOMLY WHETHER PROTON OR NEUTRON
685 CALL RMMAR( RD,JGNU,1 )
686 IPR = 0
687 INE = 0
688 TARMAS = PAMA(ITYP(2))
689 DO 114 I = 1,JGNU
690 NTOT = NTOT + 1
691 IF ( RD(I) .LE. .5D0 ) THEN
692 ITYP(NTOT) = 13
693 INE = INE + 1
694 ELSE
695 ITYP(NTOT) = 14
696 IPR = IPR + 1
697 ENDIF
698 TMAS(NTOT) = PAMA(ITYP(NTOT))
699 TARMAS = TARMAS + TMAS(NTOT)
700 EA(NTOT) = TMAS(NTOT)
701 PX(NTOT) = 0.D0
702 PY(NTOT) = 0.D0
703 PT2(NTOT) = 0.D0
704 114 CONTINUE
705
706C CALCULATE C.M. ENERGY + RAPIDITY SHIFT
707* YCMGNU = 0.5D0 * LOG( (ELAB+TARMAS+PLAB)/(ELAB+TARMAS-PLAB) )
708 YCMGNU = 0.5D0 * LOG( (EPLUSP**2 +TARMAS*EPLUSP)/
709 * (PAMA(ITYPE)**2+TARMAS*EPLUSP) )
710 DYGNU = YCM - YCMGNU
711
712C CALCULATE RAPIDITIES OF ADDITIONALLY INTERACTING TARGET NUCLEONS
713C IN THE CM SYSTEM OF NUCLEON-NUCLEON SYSTEM
714 DO 115 I = NTOT-JGNU+1,NTOT
715 YR(I) = - YCM
716 115 CONTINUE
717C SHIFT RAPIDITIES INTO CM SYSTEM OF GNU+1 MASSES
718 DO 113 I = 1,NTOT
719 YR(I) = YR(I) + DYGNU
720 113 CONTINUE
721
722C CENTER OF MASS ENERGY OF 1 PROJECTILE AND GNU TARGET NUCLEONS TO
723C BE USED IN THE JADACH FILTER.
724 ECMJAD = SQRT( PAMA(ITYPE)**2 + TARMAS**2 + 2.D0*TARMAS*ELAB )
725
726 ENDIF
727
728 ELSE
729C IN CASE OF DIFFRACTION SET THE RAPIDITY OF LEADER AND ANTILEADER
730C IN SUBROUTINE LEADDF
731 DYGNU = 0.D0
732 ECMJAD = ECMDPM
733 CALL LEADDF( IFLGLD )
734 IF ( IFLGLD .NE. 0 ) THEN
735 IF ( DEBUG ) WRITE(MDEBUG,*) ' IFLGLD=',IFLGLD
736 GOTO 30
737 ENDIF
738 ENDIF
739
740C CORRECT THE RAPIDITIES TO CONSERVE LONGITUDINAL MOMENTA AND ENERGY
741C USING THE ALGORITHM OF JADACH (SIMPLIFIED VERSION BY R. ATTALLAH)
742 CALL JADACH( ECMJAD,JADFLG )
743 IF ( JADFLG .NE. 0 ) THEN
744 IF ( DEBUG ) WRITE(MDEBUG,*) ' JADFLG=', JADFLG
745 IF ( JADFLG .GT. 0 ) GOTO 30
746 IF ( JADFLG .LT. 0 ) THEN
747 NREPRD = NREPRD + 1
748 IF ( NREPRD .GT. 10 ) GOTO 1
749 GOTO 1919
750 ENDIF
751 ENDIF
752
753
754C CALCULATE LAB ENERGIES OF SECONDARY PARTICLES FROM THE RAPIDITIES
755C INCLUDING THE ADDITIONAL TARGET NUCLEONS
756 ETOT = 0.D0
757 DO 135 I = 1,NTOT
758 YR(I) = YR(I) + YCM - DYGNU
759 EA(I) = TMAS(I) * COSH( YR(I) )
760 ETOT = ETOT + EA(I)
761 135 CONTINUE
762
763 IF ( DEBUG ) WRITE(MDEBUG,136)
764 * (I,ITYP(I),PX(I),PY(I),YR(I),EA(I),I=1,NTOT)
765 136 FORMAT(' NO ITYP PX PY YR EA'/
766 * (' ',I4,I3,1X,1P,4G13.6) )
767
768C-----------------------------------------------------------------------
769C LOOP OVER ALL SECONDARY PARTICLES AND THE LEADING PARTICLE
770C PUT THEM ON THE STACK
771 DO 139 LK = 5,8
772 SECPAR(LK) = CURPAR(LK)
773 139 CONTINUE
774
775C PROCESS LOOP
776 DO 140 J = 1,NTOT
777C REJECTION OF BACKWARD GOING PARTICLES
778 IF ( YR(J) .LE. 0.D0 ) THEN
779 IF ( DEBUG ) WRITE(MDEBUG,*)'HDPM : YR REJECT PARTICLE ',J
780 GOTO 140
781 ENDIF
782C CALCULATE THE PROPERTIES OF ALL SECONDARIES
783C PARTICLE TYPE
784 SECPAR(1) = ITYP(J)
785C CALCULATE GAMMA FACTOR
786 SECPAR(2) = EA(J) / PAMA(ITYP(J))
787C TOTAL MOMENTUM SQUARED
788 PTM = EA(J)**2 - PAMA(ITYP(J))**2
789 IF ( PT2(J) .GT. PTM ) THEN
790 IF ( DEBUG ) WRITE(MDEBUG,*)'HDPM : PT REJECT PARTICLE ',J
791 GOTO 140
792 ENDIF
793C EMISSION ZENITH ANGLE AGAINST TRAJECTORY OF PROJECTILE
794 IF ( PTM .EQ. 0.D0 ) THEN
795 COSTET = 1.D0
796 ELSE
797 COSTET = SQRT( 1.D0 - PT2(J) / PTM )
798 ENDIF
799C EMISSION AZIMUTH ANGLE
800 IF ( PX(J) .NE. 0.D0 .OR. PY(J) .NE. 0.D0 ) THEN
801 PHIJ = ATAN2( PY(J), PX(J) )
802 ELSE
803 PHIJ = 0.D0
804 ENDIF
805 CALL ADDANG( COSTHE,PHI, COSTET,PHIJ, SECPAR(3),SECPAR(4) )
806 IF ( SECPAR(3) .LT. C(29) ) THEN
807C OMIT UPWARD GOING PARTICLES
808 IF (DEBUG) WRITE(MDEBUG,*)'HDPM : ANGLE REJECT PARTICLE ',J
809 GOTO 140
810 ENDIF
811C PUT SECONDARY PARTICLES ON STACK, IF NOT GOING UPWARDS
812 IF ( J .GT. 2 ) THEN
813 CALL TSTACK
814 ELSE
815C PUT LEADER OR ANTI-LEADER ON STACK, IF NOT GOING UPWARDS
816 IF ( ITYP(J) .GT. 50 ) THEN
817C LEADER OR ANTI LEADER ARE RESONANCES AND DECAY
818 IRESPAR = IRESPAR + 1
819 IF ( IRESPAR .GE. 1000 ) THEN
820 WRITE(MONIOU,*) 'STACK OF RESDEC RANDOM NUMBERS FULL'
821 IRESPAR = 999
822 ENDIF
823 RESRAN(IRESPAR) = RDRES(J)
824C COUNTER FOR ENERGY-MULTIPLICITY MATRIX
825 MSMM = MSMM + 1
826 ENDIF
827 CALL TSTACK
828
829C CALCULATE ELASTICITY FROM ENERGY OF LEADER (REST OF RESONANCE DECAY)
830 IF ( J. EQ. 1 ) ELASTI = SECPAR(2)*PAMA(NINT(SECPAR(1)))/ELAB
831 ENDIF
832C COUNTERS FOR FIRST INTERACTION
833 IF ( FIRSTI ) THEN
834 IF ( SECPAR(1) .EQ. 7.D0 .OR. SECPAR(1) .EQ. 8.D0
835 * .OR. SECPAR(1) .EQ. 9.D0 ) THEN
836 IFINPI = IFINPI + 1
837 ELSEIF ( SECPAR(1) .EQ. 13.D0 .OR. SECPAR(1) .EQ. 14.D0
838 * .OR. SECPAR(1) .EQ. 15.D0 .OR. SECPAR(1) .EQ. 25.D0 ) THEN
839 IFINNU = IFINNU + 1
840 ELSEIF ( SECPAR(1) .EQ. 10.D0 .OR. SECPAR(1) .EQ. 11.D0
841 * .OR. SECPAR(1) .EQ. 12.D0 .OR. SECPAR(1) .EQ. 16.D0 ) THEN
842 IFINKA = IFINKA + 1
843 ELSEIF ( SECPAR(1) .GE. 71.D0 .AND. SECPAR(1) .LE. 74.D0) THEN
844 IFINET = IFINET + 1
845 ELSEIF ((SECPAR(1) .GE. 18.D0 .AND. SECPAR(1) .LE. 24.D0)
846 * .OR. (SECPAR(1) .GE. 26.D0 .AND. SECPAR(1) .LE. 32.D0))THEN
847 IFINHY = IFINHY + 1
848 ENDIF
849 ENDIF
850
851 140 CONTINUE
852
853C COUNTER FOR ENERGY-MULTIPLICITY MATRIX
854 MSMM = MSMM + NTOT - 2
855
856C FILL ELASTICITY IN MATRICES
857 MEL = MIN ( 1.D0+10.D0* MAX( 0.D0, ELASTI ) , 11.D0 )
858 MEN = MIN ( 4.D0+ 3.D0*LOG10(MAX( .1D0, EKINL )), 37.D0 )
859 IELDPM(MEN,MEL) = IELDPM(MEN,MEL) + 1
860 IELDPA(MEN,MEL) = IELDPA(MEN,MEL) + 1
861 IF ( ELASTI .LT. 1.D0 ) THEN
862 ELMEAN(MEN) = ELMEAN(MEN) + ELASTI
863 ELMEAA(MEN) = ELMEAA(MEN) + ELASTI
864 ENDIF
865
866 IF ( FIRSTI ) THEN
867 ELAST = ELASTI
868 FIRSTI = .FALSE.
869 ENDIF
870
871 IF ( DEBUG ) WRITE(MDEBUG,*) 'HDPM : ELAST=',SNGL(ELASTI),
872 * SNGL(ETOT),SNGL(ELAB)
873
874 RETURN
875 END
Note: See TracBrowser for help on using the repository browser.