]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/fluka/corevt.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Patch from Piergiorgio, concerning Dead andoes and the like.
[u/mrichter/AliRoot.git] / GEANT321 / fluka / corevt.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1995/10/24 10:19:54 cernlib
6* Geant
7*
8*
9#include "geant321/pilot.h"
10*CMZ : 3.21/02 29/03/94 15.41.41 by S.Giani
11*-- Author :
12*$ CREATE COREVT.FOR
13*COPY COREVT
14* *
15*=== corevt ===========================================================*
16* *
17 SUBROUTINE COREVT ( ZZTAR, BBTAR, KPROJ, PPROJ, EKE )
18
19#include "geant321/dblprc.inc"
20#include "geant321/dimpar.inc"
21#include "geant321/iounit.inc"
22*
23*----------------------------------------------------------------------*
24* *
25* Last change on 07-may-92 by Alfredo Ferrari, INFN-Milan *
26* *
27* This version has been developed by A. Ferrari trying to take into *
28* account a few papers about correlations between projectile colli- *
29* sions and secondary collisions. An improved (but slower) version *
30* is going on. *
31* This subroutine correlates intranuclear cascade with the number *
32* of projectile collisions sampled in the nucleus *
33* *
34* input parameters: *
35* bbtar - atomic number of the target *
36* kproj - type of the projectile *
37* pproj - momentum of the projectile *
38* *
39* output parameters (in common corinc) *
40* frainc - reduction factor for intran. cascade energy, *
41* inc. correlat. *
42* nsea+1 - number of high energy collisions in the nucleus *
43* *
44*----------------------------------------------------------------------*
45*
46#include "geant321/balanc.inc"
47#include "geant321/corinc.inc"
48#include "geant321/nucdat.inc"
49#include "geant321/parevt.inc"
50#include "geant321/part.inc"
51#include "geant321/qquark.inc"
52#include "geant321/resnuc.inc"
53*
54 PARAMETER ( PIO2 = 0.5D+00 * PIPIPI )
55 PARAMETER ( SQPI = 1.772453850905516 D+00 )
56 PARAMETER ( SQPIT3 = 3.D+00 * SQPI )
57 PARAMETER ( TWOSQT = 2.828427124746190 D+00 )
58 PARAMETER ( ECUTRF = 100.0 D+00 )
59 PARAMETER ( UMOREF = 13.83 D+00 )
60 PARAMETER ( FINCFR = 1.0 D+00 )
61 PARAMETER ( ANUC00 = 1.4 D+00 )
62 PARAMETER ( RAOB00 = 3.0 D+00 )
63 PARAMETER ( ACPAR0 = 0.5D+00 - 0.5D+00 * ( ANUC00 - 1.D+00 )
64 & / RAOB00 )
65*
66 COMMON / FKNUCO / HELP (2), HHLP (2), FTVTH (2), FINCX (2),
67 & EKPOLD (2), BBOLD, ZZOLD, SQROLD, ASEASQ,
68 & FSPRED, FEX0RD
69 DIMENSION KPA (39)
70 DIMENSION FRA(2,110)
71 REAL RNDM(3)
72 SAVE BBONE, FRA, KPA, XIXIXI, ANCOSQ
73 LOGICAL LNUCNW, LUFFA
74 DATA BBONE / 1.D+00 /
75 DATA KPA/1,1,3,3,3,3,3,2,2,3,3,3,3,3,3,3,2,2,3,1,1,2,3,3,3,3,
76 & 3,3,3,3,1,2,1,2,2,1,1,1,1/
77* Reduction factors for intran. cascade energy, taken from Alsmiller
78* Incoming baryons
79 DATA (FRA(1,I), I=1,107) / .048D0,.076D0,0.10D0,.12D0,
80 * .14D0,.16D0,.17D0,.19D0,.21D0,.22D0,.24D0,
81 1 .25D0,.26D0,.28D0,.29D0,.30D0,.315D0,.33D0,.34D0,
82 * .35D0,.36D0,.38D0,.39D0,.40D0,.41D0,.42D0,
83 2 .43D0,.44D0,.46D0,.47D0,.48D0,.49D0,.50D0,.51D0,
84 * .52D0,.53D0,.54D0,.55D0,.56D0,.57D0,.58D0,.59D0,
85 3 .60D0,.61D0,.62D0,.63D0,.635D0,.64D0,.65D0,.66D0,
86 * .67D0,.675D0,.68D0,.69D0,.70D0,.71D0,.715D0,
87 4 .72D0,.725D0,.73D0,.735D0,.74D0,.75D0,.755D0,
88 * .76D0,.767D0,.77D0,.77D0,.775D0,.78D0,.783D0,
89 5 .786D0,.79D0,.795D0,.80D0,.805D0,.81D0,.812D0,
90 * .815D0,.82D0,.822D0,.824D0,.825D0,.825D0,.83D0,
91 6 .832D0,.834D0,.836D0,.838D0,.84D0,.843D0,.836D0,
92 * .849D0,.852D0,.855D0,.856D0,.857D0,.858D0,
93 7 .859D0,.860D0,.862D0,.864D0,.866D0,.868D0,.870D0,
94 * .872D0,.874D0/
95* Incoming mesons
96 DATA (FRA(2,I), I=1,63) / .048D0,.076D0,0.10D0,.12D0,
97 * .14D0,.16D0,.17D0,.19D0,.21D0,.22D0,.24D0,
98 1 .25D0,.262D0,.274D0,.285D0,.295D0,.305D0,.315D0,
99 * .327D0,.340D0,.345D0,.350D0,.360D0,
100 2 .370D0,.380D0,.385D0,.390D0,.398D0,.406D0,.415D0,
101 * .420D0,.425D0,.430D0,.435D0,.440D0,.446D0,
102 3 .452D0,.458D0,.464D0,.470D0,.474D0,.478D0,.482D0,
103 * .486D0,.490D0,.495D0,.500D0,.505D0,.510D0,
104 4 .515D0,.519D0,.523D0,.526D0,.520D0,.533D0,.537D0,
105 * .540D0,.543D0,.547D0,.550D0,.555D0,.559D0,
106 5 .5625D0 /
107*
108* The calculation of number of collisions inside nucleus is moved
109* from subroutine nucevt; nsea is stored in common corinc
110*
111 IBPROJ = IBAR (KPROJ)
112* Get the "paprop" index
113 IJ = KPTOIP (KPROJ)
114 IJJ = KPA (IJ)
115* +-------------------------------------------------------------------*
116* | Incoming baryons
117 IF ( IBPROJ .NE. 0 ) THEN
118* | +----------------------------------------------------------------*
119* | |
120 IF ( IBTAR .LE. 107 ) THEN
121 FRAC = FRA ( 1, IBTAR )
122* | |
123* | +----------------------------------------------------------------*
124* | |
125 ELSE IF (IBTAR .LE. 206) THEN
126 FRAC = 0.739D+00 + 0.00126D+00 * BBTAR
127* | |
128* | +----------------------------------------------------------------*
129* | |
130 ELSE
131 FRAC = 1.D+00
132 END IF
133* | |
134* | +----------------------------------------------------------------*
135* | +----------------------------------------------------------------*
136* | | This is a very simple patch to slightly increase the cascade
137* | | yield at low energies for antibaryons. The factor 0.5 to
138* | | take into account very roughly that not all interactions will
139* | | result in complete annihilation
140 IF ( IBPROJ .LT. 0 ) THEN
141 EKEFF = EKE + 0.5D+00 * AMDISC (KPROJ)
142* | |
143* | +----------------------------------------------------------------*
144* | |
145 ELSE
146 EKEFF = EKE
147 END IF
148* | |
149* | +----------------------------------------------------------------*
150 STRRED = 0.D+00
151* | +----------------------------------------------------------------*
152* | | Apply a 50 % reduction in the intranuclear cascade probability
153* | | for each s/sbar quark of the projectile
154 DO 100 IQ = 1,3
155 STRRED = STRRED + 0.1666666666666667D+00
156 & * ABS ( IQSCHR ( MQUARK (IQ,IJ) ) )
157 100 CONTINUE
158* | |
159* | +----------------------------------------------------------------*
160* | Make the strangeness reduction effective only at low energies
161* | (where secondaries are few ==> the s composition of the projecti-
162* | le is relevant)
163 STRRED = STRRED * AM (KPROJ) / ( EKE + AM (KPROJ) )
164 FRAC = ( 1.D+00 - STRRED ) * FRAC
165* |
166* +-------------------------------------------------------------------*
167* | Incoming mesons
168 ELSE
169* | +----------------------------------------------------------------*
170* | |
171 IF ( IBTAR .LE. 63 ) THEN
172 WEIGH1 = 1.D+00 / ( 1.D+00 + ( 60.D+00 / EKE ) )
173 FRAC = FRA ( 2, IBTAR ) * ( ( 1.D+00 + 0.1333D+00 * MAX (
174 & 0.D+00, BBTAR - 35.D+00 ) / 28.D+00 ) * WEIGH1 +
175 & 1.D+00 - WEIGH1 )
176* | |
177* | +----------------------------------------------------------------*
178* | |
179 ELSE IF ( IBTAR .LE. 107 ) THEN
180 WEIGH1 = 1.D+00 / ( 1.D+00 + ( 60.D+00 / EKE ) )
181 FRAC = ( 0.75D+00 + WEIGH1 * 0.1D+00 ) * FRA ( 1, IBTAR )
182* | |
183* | +----------------------------------------------------------------*
184* | |
185 ELSE IF ( IBTAR .LE. 206 ) THEN
186 WEIGH1 = 1.D+00 / ( 1.D+00 + ( 60.D+00 / EKE ) )
187 FRAC = ( 0.6279D+00 + 0.001077D+00 * BBTAR ) * WEIGH1
188 & + ( 1.D+00 - WEIGH1 ) * ( 0.554D+00 + 0.00095D+00
189 & * BBTAR )
190* | |
191* | +----------------------------------------------------------------*
192* | |
193 ELSE
194 WEIGH1 = 1.D+00 / ( 1.D+00 + ( 60.D+00 / EKE ) )
195 FRAC = 0.75D+00 + 0.1D+00 * WEIGH1
196 END IF
197* | |
198* | +----------------------------------------------------------------*
199 EKEFF = EKE
200 STRRED = 0.D+00
201* | +----------------------------------------------------------------*
202* | | Apply a 50 % reduction in the intranuclear cascade probability
203* | | for each s/sbar quark of the projectile
204 DO 150 IQ = 1,2
205 STRRED = STRRED + 0.25D+00
206 & * ABS ( IQSCHR ( MQUARK (IQ,IJ) ) )
207 150 CONTINUE
208* | |
209* | +----------------------------------------------------------------*
210* | Make the strangeness reduction effective only at low energies
211* | (where secondaries are few ==> the s composition of the projecti-
212* | le is relevant)
213 STRRED = STRRED * AM (KPROJ) / ( EKE + AM (KPROJ) )
214 FRAC = ( 1.D+00 - STRRED ) * FRAC
215 END IF
216* |
217* +-------------------------------------------------------------------*
218* The following is a reduction factor to bring Hannes's parametri-
219* zations for the average energy carried by cascade nucleons in
220* better agreement with experimental data
221 FSPRED = FSPRD0 + ( 1.D+00 - FSPRD0 ) * MAX ( 10.D+00 - EKE,
222 & 0.D+00 ) / 10.D+00
223 FEX0RD = FSPRD0**ABS(IBPROJ)
224* +-------------------------------------------------------------------*
225* | Check whether it is the same target nucleus of the last call
226 IF ( BBTAR .NE. BBOLD .OR. ZZTAR .NE. ZZOLD ) THEN
227 LNUCNW = .TRUE.
228 SQRAMS = SQRT ( BBTAR )
229 ATO1O3 = BBTAR**0.3333333333333333D+00
230* ZTO1O3 = ZZTAR**0.3333333333333333D+00
231 BBOLD = BBTAR
232 ZZOLD = ZZTAR
233 SQROLD = SQRAMS
234 HKAP = BBTAR**2 / ( ZZTAR**2 + ( BBTAR - ZZTAR )**2 )
235 HHLP (1) = ( HKAP * ZZTAR )**0.3333333333333333D+00 / ATO1O3
236 HHLP (2) = ( HKAP * ( BBTAR - ZZTAR ) )
237 & **0.3333333333333333D+00 / ATO1O3
238 RDSNUC = R0NUCL * ATO1O3
239 RDSCOU = RCCOUL * ATO1O3
240 FLKCOU = DOST ( 1, ZZTAR )
241* | Coulomb barrier "a la Evap"
242 VEFFNU (1) = FLKCOU * COULPR * ZZTAR / RDSCOU
243 VEFFNU (2) = 0.D+00
244 AMRCAV = MAX ( 1.D+00, ( BBTAR - 2.D+00 ) ) * AMUC12
245 AMRCSQ = AMRCAV * AMRCAV
246* | +----------------------------------------------------------------*
247* | |
248 DO 3000 I = 1, 2
249 PFRMMX (I) = HHLP (I) * APFRMX
250 P2HELP = PFRMMX (I)**2
251 EFRMMX (I) = SQRT ( P2HELP + AMNUSQ (I) ) - AMNUCL (I)
252 EFRMAV (I) = 0.3D+00 * P2HELP / AMNUCL (I) * ( 1.D+00
253 & - P2HELP / ( 5.6D+00 * AMNUSQ (I) ) )
254 ERCLAV (I) = 0.3D+00 * P2HELP / AMRCAV * ( 1.D+00
255 & - P2HELP / ( 5.6D+00 * AMRCSQ ) )
256 V0WELL (I) = EFRMMX (I) + EBNDNG (I)
257 VEFFNU (I) = VEFFNU (I) + V0WELL (I)
258 ERCLMX = 0.5D+00 * P2HELP / AMRCAV * ( 1.D+00
259 & - 0.25D+00 * P2HELP / AMRCSQ )
260 EKMNNU (I) = VEFFNU (I) + EBNDNG (I)
261 & - EFRMMX (I) + ERCLMX
262 EKMXNU (I) = VEFFNU (I) + EBNDNG (I)
263 EKMNAV (I) = VEFFNU (I) + EBNDNG (I)
264 & - EFRMAV (I) + ERCLAV (I)
265 ESWELL (I) = EBNDNG (I) + V0WELL (I)
266 & - EFRMAV (I) + ERCLAV (I)
267 FTVTH (I) = ( EKMNNU (I) / EFRMMX (I) )**3
268 FTVTH (I) = 0.4D+00 * SQRT ( FTVTH (I) )
269 FINCUP (I) = AKEKA ( I, EKEFF, BBTAR ) * FSPRED
2703000 CONTINUE
271* | |
272* | +----------------------------------------------------------------*
273* |
274* +-------------------------------------------------------------------*
275* | It is the same target nucleus of the last call to Nucevv/Nucriv
276 ELSE
277 LNUCNW = .FALSE.
278 SQRAMS = SQROLD
279 FINCUP (1) = AKEKA ( 1, EKEFF, BBTAR ) * FSPRED
280 FINCUP (2) = AKEKA ( 2, EKEFF, BBTAR ) * FSPRED
281 END IF
282* |
283* +-------------------------------------------------------------------*
284 CALL NIZL (IJ, BBTAR, EKE, PPROJ, SIHA, ZLA)
285 CALL NIZL (IJ, BBONE, EKE, PPROJ, SIHN, ZLP)
286 ANUAV = BBTAR * SIHN / SIHA
287* Anuav= average number of collisions in nucleus
288 ANUSEA = ANUAV - 1.D+00
289* +-------------------------------------------------------------------*
290* | Call the function sampling the distribution for the
291* | number of high energy collisions
292 IF ( ANUSEA .GT. 0.D+00 .AND. PPROJ .GT. 5.D+00 ) THEN
293 EXPLAM = EXP ( -ANUSEA )
294 NSEA = NUDISV ( ANUAV, IBPROJ, EXPLAM, ASEASQ, APOWER,
295 & PRZERO ) - 1
296 NSEA = MIN ( NSEA, IBTAR - 1 )
297 LUFFA = .FALSE.
298* |
299* +-------------------------------------------------------------------*
300* |
301 ELSE
302 ASEASQ = 0.D+00
303 ANUSEA = 0.D+00
304 APOWER = 1.D+00
305 ANUAV = 1.D+00
306 PRZERO = 1.D+00
307 ANUCSQ = 1.D+00
308 NSEA = 0
309 RATOLD = 0.D+00
310 ACFACT = 0.D+00
311 BCFACT = 2.D+00
312 ANUC11 = ANUC00
313 ANUCOR = ANUC00
314 ACPARM = 0.D+00
315 RRPCO = 0.D+00
316 LUFFA = .TRUE.
317 END IF
318* |
319* +-------------------------------------------------------------------*
320 NSEALD = NSEA
321* +-------------------------------------------------------------------*
322* | Check if the parameterized intranuclear cascade is requested
323 IF ( LINCTV ) THEN
324 IF ( LUFFA ) GO TO 3500
325 CALL GRNDM(RNDM,1)
326 RRPCR = RNDM (1)
327* | +----------------------------------------------------------------*
328* | | Check for negative <n_sea^2>
329 IF ( ASEASQ .LT. 0.D+00 ) THEN
330 ANUCOR = ANUC00 / ( 1.D+00 + 2.D+00 / PPROJ )
331 ASEASQ = ANUCOR * ANUAV**2 - 2.D+00 * ANUAV + 1.D+00
332 ANUCSQ = ANUCOR * ANUAV**2
333 ACPARM = ACPAR0
334 ANUC11 = ANUC00
335 PRZERO = EXPLAM
336* | |
337* | +----------------------------------------------------------------*
338* | |
339 ELSE
340 ANUCSQ = ASEASQ + 2.D+00 * ANUAV - 1.D+00
341 ANUC11 = ANUCSQ / ANUAV**2
342* | | +-------------------------------------------------------------*
343* | | | Check if a power different from 2 must be used for the
344* | | | cascade particle distributions
345* | | | Power = 2:
346 IF ( .NOT. LPOWER ) THEN
347 APOWER = ANUCSQ
348* | | |
349* | | +-------------------------------------------------------------*
350* | | |
351 ELSE
352 IPOWER = NINT (-DPOWER)
353 IF ( IPOWER .EQ. 17 ) THEN
354 FEX0RD = FEX0RD * FSPRD0
355 ELSE IF ( IPOWER .EQ. 8 .OR. IPOWER .EQ. 1 .OR. IPOWER
356 & .GE. 12 ) THEN
357 FEX0RD = FEX0RD * FSPRD0
358 END IF
359 END IF
360* | | |
361* | | +-------------------------------------------------------------*
362 ANUC11 = MIN ( ANUC11, ANUC00 )
363 ANUCOR = ANUC11 / ( 1.D+00 + 5.D+00 * ( ANUC11 - 1.D+00 )
364 & / PPROJ )
365 ANUCOR = MAX ( ANUCOR, ONEONE )
366 BBCOF = ( ANUC11 * ANUAV**2 - 2.D+00 * ANUCOR * ANUAV**2
367 & + 2.D+00 * ANUCOR * ANUAV - 1.D+00 )
368 BBCOF = ABS (BBCOF)
369 AACOF = ANUCOR * ( ANUAV - 1.D+00 )**2
370 CCCOF = - ANUAV**2 * ( ANUC11 - ANUCOR )
371 IF ( AACOF - CCCOF .GT. ANGLGB ) THEN
372 RRPCO = 0.5D+00 * ( - BBCOF + SQRT ( BBCOF**2 - 4.D+00
373 & * AACOF * CCCOF ) ) / AACOF
374 ELSE
375 RRPCO = 0.D+00
376 END IF
377 RRPCO = MIN ( ONEONE, RRPCO )
378 RRPCO = MAX ( ZERZER, RRPCO )
379 RRPCO = 1.D+00 - RRPCO
380 END IF
381* | |
382* | +----------------------------------------------------------------*
383* | Supply the fraction of the total kinetic energy to be
384* | used for intranuclear cascade nucleons
385 3500 CONTINUE
386 EKUPNU (1) = AINEL (IJJ, 6, EKEFF, BBTAR, SQRAMS) * EKEFF
387 & * FSPRED
388 EKUPNU (2) = AINEL (IJJ, 7, EKEFF, BBTAR, SQRAMS) * EKEFF
389 & * FSPRED
390 EINCP = FRAC * EKUPNU (1)
391 EINCN = FRAC * EKUPNU (2)
392 FRAMAX = FRAC
393* | +----------------------------------------------------------------*
394* | | Now start Fincup (i) calculation!!!!
395 DO 4000 I = 1, 2
396 EKPOLD (I) = EKUPNU (I)
397 EKUPNU (I) = MAX ( FRAMAX * EKUPNU (I), EKMXNU (I)/TWOTHI
398 & , FOUFOU * FINCUP (I) )
399 EKUPNU (I) = MIN ( EKUPNU (I), FOUFOU * FINCUP (I) )
400 AHELP = - ( EKUPNU (I) - EKMNAV (I) ) / FINCUP (I)
401 CHELP = EKMNAV (I) / FINCUP (I)
402 DHELP = EKUPNU (I) / FINCUP (I)
403 FINC0 = FINCFR + ESWELL (I) / FINCUP (I)
404 BHELP = EXP ( AHELP / FINC0 )
405 FINCA = FINC0
406 FUNCA = - ( CHELP - DHELP * BHELP ) / ( 1.D+00 - BHELP )
407 FINCB = 2.D+00 * FINC0
408 BHELP = EXP ( AHELP / FINCB )
409 FUNCB = FINC0 - FINCB - ( CHELP - DHELP * BHELP ) /
410 & ( 1.D+00 - BHELP )
411 ICOU = 1
412 FINCLD = FINC0
413* | | +-------------------------------------------------------------*
414* | | |
4153800 CONTINUE
416 FINCX (I) = FINCA - FUNCA * ( FINCB - FINCA ) /
417 & ( FUNCB - FUNCA )
418* | | | +----------------------------------------------------------*
419* | | | |
420 IF ( ABS ( FINCX (I) - FINCLD ) .GT. 0.03D+00 .AND.
421 & ICOU .LT. 20 ) THEN
422 ICOU = ICOU + 1
423 FINCLD = FINCX (I)
424* | | | | +-------------------------------------------------------*
425* | | | | |
426 IF ( ABS (FUNCA) .LT. ABS (FUNCB) ) THEN
427 FINCB = FINCLD
428 BHELP = EXP ( AHELP / FINCB )
429 FUNCB = FINC0 - FINCB - ( CHELP - DHELP * BHELP
430 & ) / ( 1.D+00 - BHELP )
431* | | | | |
432* | | | | +-------------------------------------------------------*
433* | | | | |
434 ELSE
435 FINCA = FINCLD
436 BHELP = EXP ( AHELP / FINCA )
437 FUNCA = FINC0 - FINCA - ( CHELP - DHELP * BHELP
438 & ) / ( 1.D+00 - BHELP )
439 END IF
440* | | | | |
441* | | | | +-------------------------------------------------------*
442 GO TO 3800
443* | | | |
444* | | +-<|--<--<--<--<--<
445 END IF
446* | | | |
447* | | | +----------------------------------------------------------*
448* | | |
449* | | +-------------------------------------------------------------*
450* | | |
451* | | +-------------------------------------------------------------*
452 ESLOPE (I) = FINCUP (I) * FINCX (I)
453 AHELP = EXP ( - EKPOLD (I) / FINCUP (I) )
454 EKNOLD = FINCUP (I) - EKPOLD (I) * AHELP /
455 & ( 1.D+00 - AHELP )
456 EXMNAV (I) = EXP ( - EKMNAV (I) / ESLOPE (I) )
457 EXMNNU (I) = EXP ( - EKMNNU (I) / ESLOPE (I) )
458 EXUPNU (I) = EXP ( - EKUPNU (I) / ESLOPE (I) )
459 EKINAV (I) = ESLOPE (I) + ( EKMNAV (I) * EXMNAV (I) -
460 & EKUPNU (I) * EXUPNU (I) ) / ( EXMNAV (I) -
461 & EXUPNU (I) )
462 EKPOLD (I) = EKPOLD (I) * FRAC
463 FINCUP (I) = EKINAV (I) / EKNOLD
4644000 CONTINUE
465* | |
466* | +----------------------------------------------------------------*
467 EINCP = EINCP * FINCUP (1)
468 EINCN = EINCN * FINCUP (2)
469 AGREYP = EINCP / EKINAV (1)
470 AGREYN = EINCN / EKINAV (2)
471 AGREYT = AGREYP + AGREYN
472 CALL GRNDM(RNDM,1)
473 RNPCO = RNDM ( 1 )
474* | +----------------------------------------------------------------*
475* | | Check if we have to sample from a distribution with <n_casc>
476* | | prop. <nu^power> or to <nu>
477 IF ( RNPCO .LT. RRPCO ) THEN
478* | | +-------------------------------------------------------------*
479* | | | Power .ne. 2:
480 IF ( LPOWER ) THEN
481* | | | +----------------------------------------------------------*
482* | | | | 1**Y = 1
483 IF ( NSEA .EQ. 0 ) THEN
484 IF ( IPOWER .LT. 7 ) THEN
485 ANCOSQ = 1.D+00
486 ELSE
487 ANCOSQ = FPOWER ( IPOWER, 1, ANUAV )
488 END IF
489* | | | |
490* | | | +----------------------------------------------------------*
491* | | | | The exponent has a fixed value
492 ELSE IF ( DPOWER .GT. 0.D+00 ) THEN
493 ANCOSQ = ( NSEA + 1 )**DPOWER
494* | | | |
495* | | | +----------------------------------------------------------*
496* | | | | The function Fpower supplies the exponent
497 ELSE
498 IF ( IPOWER .LT. 7 ) THEN
499 ANCOSQ = ( NSEA + 1 )**FPOWER ( IPOWER,
500 & NSEA + 1, ANUAV )
501 ELSE IF ( IPOWER .LT. 11 ) THEN
502 ANCOSQ = ( NSEA + 1 ) * FPOWER ( IPOWER,
503 & NSEA + 1, ANUAV )
504 ELSE
505 ANCOSQ = ( NSEA + 1 )**FPOWER ( IPOWER,
506 & NSEA + 1, ANUAV )
507 END IF
508 END IF
509* | | | |
510* | | | +----------------------------------------------------------*
511* | | |
512* | | +-------------------------------------------------------------*
513* | | | Power = 2:
514 ELSE
515 ANCOSQ = ( NSEA + 1 )**2
516 END IF
517* | | |
518* | | +-------------------------------------------------------------*
519 XIXIXI = AGREYP / ( AGREYP + APOWER )
520* | |
521* | +----------------------------------------------------------------*
522* | |
523 ELSE
524 ANCOSQ = NSEA + 1
525 XIXIXI = AGREYP / ( AGREYP + ANUAV )
526 END IF
527* | |
528* | +----------------------------------------------------------------*
529 RRPC0 = ( 1.D+00 - XIXIXI )**ANCOSQ
530 RRPCR = RRPC0
531 CALL GRNDM(RNDM,1)
532 RNPCR = RNDM (1)
533 IF ( RRPCR .GE. RNPCR ) THEN
534 NGREYP = 0
535 ELSE
536 DO 4600 I = 1, ICHTAR
537 RRPC0 = RRPC0 * ( I - 1 + ANCOSQ ) * XIXIXI
538 & / I
539 RRPCR = RRPCR + RRPC0
540 IF ( RNPCR .LE. RRPCR ) GO TO 4700
541 4600 CONTINUE
542 I = I - 1
543 IF ( BBTAR .GT. 12.D+00 )
544 & WRITE (LUNERR,*)' *** RRPCR,I,NSEA,ANCOSQ*XIXIXI ***',
545 & RRPCR,I,NSEA,ANCOSQ*XIXIXI
546 4700 CONTINUE
547 NGREYP = I
548 END IF
549* | +----------------------------------------------------------------*
550* | | Check if we have to sample from a distribution with <n_casc>
551* | | prop. <nu^power> or to <nu>
552 IF ( RNPCO .LT. RRPCO ) THEN
553 XIXIXI = AGREYN / ( AGREYN + APOWER )
554* | |
555* | +----------------------------------------------------------------*
556* | |
557 ELSE
558 ANCOSQ = NSEA + 1
559 XIXIXI = AGREYN / ( AGREYN + ANUAV )
560 END IF
561* | |
562* | +----------------------------------------------------------------*
563 RRNC0 = ( 1.D+00 - XIXIXI )**ANCOSQ
564 RRNCR = RRNC0
565* | Correlate cascade neutrons to cascade protons:
566* | No correlation
567 CALL GRNDM(RNDM,1)
568 RNNCR = RNDM (1)
569 RN1GSC = RNPCR
570 RN2GSC = RNNCR
571 IF ( RRNCR .GE. RNNCR ) THEN
572 NGREYN = 0
573 ELSE
574 DO 4650 I = 1, IBTAR - ICHTAR
575 RRNC0 = RRNC0 * ( I - 1 + ANCOSQ ) * XIXIXI
576 & / I
577 RRNCR = RRNCR + RRNC0
578 IF ( RNNCR .LE. RRNCR ) GO TO 4750
579 4650 CONTINUE
580 I = I - 1
581 IF ( BBTAR .GT. 12.D+00 )
582 & WRITE (LUNERR,*)' *** RRNCR,I,NSEA,ANCOSQ*XIXIXI ***',
583 & RRNCR,I,NSEA,ANCOSQ*XIXIXI
584 4750 CONTINUE
585 NGREYN = I
586 END IF
587 NGREYT = NGREYN + NGREYP
588 NGREYN = 0
589 NGREYP = 0
590 PROBPR = AGREYP / ( AGREYP + AGREYN )
591 DO 4760 I = 1, NGREYT
592 CALL GRNDM(RNDM,1)
593 RNDPPR = RNDM (1)
594 IF ( RNDPPR .LT. PROBPR .AND. NGREYP .LT. ICHTAR )
595 & THEN
596 NGREYP = NGREYP + 1
597 ELSE IF ( NGREYN .LT. IBTAR - ICHTAR ) THEN
598 NGREYN = NGREYN + 1
599 ELSE
600 NGREYP = NGREYP + 1
601 END IF
602 4760 CONTINUE
603 FRAMAX = FRAC
604* | The number of grey particles is now correlated to the number
605* | of high energy collisions
606 EINCP = NGREYP * EKINAV (1)
607 EINCN = NGREYN * EKINAV (2)
608 TVTENT = ( NSEA + 1 ) * ( AV0WEL - AEFRMA )
609* |
610* +-------------------------------------------------------------------*
611* |
612 ELSE
613 EINCP = 0.D+00
614 EINCN = 0.D+00
615 EKUPNU (1) = 1.D+01
616 EKUPNU (2) = 1.D+01
617 ESLOPE (1) = 0.D+00
618 ESLOPE (2) = 0.D+00
619 EKINAV (1) = 1.D+10
620 EKINAV (2) = 1.D+10
621 FRAC = 0.D+00
622 PCR = 1.D+00
623 FRAINC = 0.D+00
624 TVTENT = 0.D+00
625 END IF
626* |
627* +-------------------------------------------------------------------*
628 PCROLD = PCR
629 EKTRIA = EKE - EINCP - EINCN - TVTENT
630 IF ( EKTRIA .LE. 0.5D+00 * EKE ) THEN
631 EKTRIA = 0.5D+00 * EKE
632 END IF
633 ETRIAL = EKTRIA + AM (KPROJ)
634* +-------------------------------------------------------------------*
635* |
636 IF ( ETRIAL .LT. ECUTRF ) THEN
637 PTRIAL = SQRT ( ETRIAL**2 - AM (KPROJ)**2 )
638 XCUTFF = 0.3D+00 * ETHSEA / EKTRIA
639* XCUTFF = 0.3D+00 / PTRIAL
640* |
641* +-------------------------------------------------------------------*
642* |
643 ELSE
644 PTRIAL = ETRIAL
645 UMO = SQRT ( 2.D+00*AM(1) * ETRIAL + AM (KPROJ)**2 + AM(1)**2 )
646 XCUTFF = 0.30D+00 * ETHSEA * UMO / ( UMOREF * EKTRIA )
647 END IF
648* |
649* +-------------------------------------------------------------------*
650 200 CONTINUE
651 IF ( NSEA .EQ. 0 ) GO TO 1000
652* *** Sample X-values of nsea sea-quark-antiquark pairs
653 NC3 = 0
654* +-------------------------------------------------------------------*
655* |
656 300 CONTINUE
657* *** Sea distribution X**(-1)*(1-X)**5
658 NC3 = NC3 + 1
659* | +----------------------------------------------------------------*
660* | |
661 IF ( NC3 .GT. 10 ) THEN
662 NSEA = NSEA - 1
663 GO TO 200
664 END IF
665* | |
666* | +----------------------------------------------------------------*
667 XO = 1.D+00
668 UNOSEA = 2.0D+00
669 GAMSEA = 0.05D+00
670 XSEAMX = 1.0D+00
671* | +----------------------------------------------------------------*
672* | |
673 DO 400 I=1,NSEA
674 NCOU= 0
675 22 CONTINUE
676 NCOU = NCOU + 1
677* | | +-------------------------------------------------------------*
678* | | |
679 IF ( NCOU .LT. 11 ) THEN
680 XSEA(I) = BETRST ( GAMSEA, UNOSEA, XCUTFF, XSEAMX )
681 IF ( XSEA(I) .LT. XCUTFF ) GO TO 22
682 END IF
683* | | |
684* | | +-------------------------------------------------------------*
685 NCOU = 0
686 23 CONTINUE
687 NCOU = NCOU + 1
688* | | +-------------------------------------------------------------*
689* | | |
690 IF ( NCOU .LT. 11 ) THEN
691 XASEA(I) = BETRST ( GAMSEA, UNOSEA, XCUTFF, XSEAMX )
692 IF ( XASEA(I) .LT. XCUTFF ) GO TO 23
693 END IF
694* | | |
695* | | +-------------------------------------------------------------*
696 XO = XO * ( 1.D+00 - XSEA (I) - XASEA(I) )
697 400 CONTINUE
698* | |
699* | +----------------------------------------------------------------*
700* | +----------------------------------------------------------------*
701* | |
702 IF ( XO * PTRIAL .LT. 4.D+00 ) THEN
703 NSEA = NSEA - 1
704 IF ( NSEA .LE. 0 ) GO TO 1000
705 GO TO 300
706* | |
707* +--+--<--<--<--<--< go to resample
708* | |
709 END IF
710* | |
711* | +----------------------------------------------------------------*
712 ILOW = 0
713* | +----------------------------------------------------------------*
714* | |
715 DO 500 I = 1, NSEA
716 EMSEA = EKTRIA * ( XSEA (I) + XASEA (I) )
717* | | +-------------------------------------------------------------*
718* | | |
719 IF ( EMSEA .LE. ETHSEA + AM (23) ) THEN
720 ILOW = ILOW + 1
721* | | | +----------------------------------------------------------*
722* | | | |
723 IF ( I .LT. NSEA ) THEN
724 II = I + 1
725 XSEA (II) = XSEA (II) + XSEA (I)
726 XASEA (II) = XASEA (II) + XASEA (I)
727* | | | |
728* | | | +----------------------------------------------------------*
729* | | | |
730 ELSE IF ( I - ILOW .GT. 0 ) THEN
731 II = I - ILOW
732 EKTRIA = EKTRIA * ( 1.D+00 - XSEA(I) - XASEA (I) )
733 XSEA (II) = XSEA (II) + XSEA (I)
734 XASEA (II) = XASEA (II) + XASEA (I)
735 END IF
736* | | | |
737* | | | +----------------------------------------------------------*
738* | | |
739* | | +-------------------------------------------------------------*
740* | | |
741 ELSE
742 XSEA (I-ILOW) = XSEA (I)
743 XASEA (I-ILOW) = XASEA (I)
744 EKTRIA = EKTRIA * ( 1.D+00 - XSEA(I) - XASEA (I) )
745 END IF
746* | | |
747* | | +-------------------------------------------------------------*
748 500 CONTINUE
749* | |
750* | +----------------------------------------------------------------*
751 NSEA = NSEA - ILOW
752* |
753* +-------------------------------------------------------------------*
7541000 CONTINUE
755* +-------------------------------------------------------------------*
756* | Now sample the target nucleons for the high energy collisions
757 DO 1200 I = 1, NSEA + 1
758 ZAPU = ZNOW / ANOW
759* | +----------------------------------------------------------------*
760* | | Wounded nucleon selection:
761* | | Kt is the index of the target nucleon (1=proton, 8=neutron)
762 CALL GRNDM(RNDM,1)
763 IF ( RNDM(1) .LE. ZAPU ) THEN
764 IJTARG (I) = 1
765 ZNOW = ZNOW - 1.D0
766 KTARP = KTARP + 1
767* | |
768* | +----------------------------------------------------------------*
769* | |
770 ELSE
771 IJTARG (I) = 8
772 KTARN = KTARN + 1
773 END IF
774* | |
775* | +----------------------------------------------------------------*
776 ANOW = ANOW - 1.D0
7771200 CONTINUE
778* |
779* +-------------------------------------------------------------------*
780 ZNCOLL = KTARP
781 ANCOLL = NSEA + 1
782 AHELP = EKMNNU (1) / ESLOPE (1)
783 AHELP = AHELP * FTVTH (1) * ( 1.D+00 + 0.3333333333333333D+00
784 & * AHELP ) / ( 1.D+00 - EXUPNU (1) )
785 AHELPP = AHELP
786 AHELP = EKMNNU (2) / ESLOPE (2)
787 AHELP = AHELP * FTVTH (2) * ( 1.D+00 + 0.3333333333333333D+00
788 & * AHELP ) / ( 1.D+00 - EXUPNU (2) )
789 AHELPN = AHELP
790 FEXTRA = 0.3D+00 * AGREYP
791 TVCHCP = EFRMMX (1) - EFRMAV (1) + EBNDNG (1)
792 TVGRYP = AHELPP * ( EINCP + FEXTRA * EKINAV (1) ) * AGREYP
793 & / ( AGREYP + FEXTRA ) * FEX0RD
794 NHELP = INT ( TVGRYP / TVCHCP )
795 TVGRE0 = NHELP * TVCHCP
796 PROB0 = ( TVGRYP - TVGRE0 ) / TVCHCP
797 CALL GRNDM(RNDM,1)
798 IF ( RNDM (1) .LT. PROB0 ) THEN
799 CALL GRNDM(RNDM,3)
800 P2HELP = ( PFRMMX (1) * MAX ( RNDM (1),
801 & RNDM (2), RNDM (3) ) )**2
802 TVGRYP = 0.5D+00 * P2HELP / AMNUCL (1) * ( 1.D+00
803 & - 0.25D+00 * P2HELP / AMNUSQ (1) )
804 TVGRYP = EFRMMX (1) - TVGRYP + EBNDNG (1) + TVGRE0
805 ELSE
806 TVGRYP = TVGRE0
807 END IF
808 FEXTRA = 0.3D+00 * AGREYN
809 TVCHCN = EFRMMX (2) - EFRMAV (2) + EBNDNG (2)
810 TVGRYN = AHELPN * ( EINCN + FEXTRA * EKINAV (2) ) * AGREYN
811 & / ( AGREYN + FEXTRA ) * FEX0RD
812 NHELP = INT ( TVGRYN / TVCHCN )
813 TVGRE0 = NHELP * TVCHCN
814 PROB0 = ( TVGRYN - TVGRE0 ) / TVCHCN
815 CALL GRNDM(RNDM,1)
816 IF ( RNDM (1) .LT. PROB0 ) THEN
817 CALL GRNDM(RNDM,3)
818 P2HELP = ( PFRMMX (2) * MAX ( RNDM (1),
819 & RNDM (2), RNDM (3) ) )**2
820 TVGRYN = 0.5D+00 * P2HELP / AMNUCL (2) * ( 1.D+00
821 & - 0.25D+00 * P2HELP / AMNUSQ (2) )
822 TVGRYN = EFRMMX (2) - TVGRYN + EBNDNG (2) + TVGRE0
823 ELSE
824 TVGRYN = TVGRE0
825 END IF
826 TVGRE0 = TVGRYP + TVGRYN
827 TVGREY = 0.D+00
828 NDIFFT = NGREYT - NINT (ANOW)
829 IF ( NINT (ZNOW) - NGREYP .LT. 0 ) THEN
830 NDIFFP = NGREYP - NINT ( ZNOW )
831 NGREYP = NGREYP - NDIFFP
832 EINCP = EINCP - NDIFFP * EKINAV (1)
833 NGREYN = NGREYN + NDIFFP
834 EINCN = EINCN + NDIFFP * EKINAV (2)
835 IF ( NINT (ANOW-ZNOW) - NGREYN .LT. 0 ) THEN
836 NDIFFN = NGREYN - NINT ( ANOW - ZNOW )
837 NGREYN = NGREYN - NDIFFN
838 EINCN = EINCN - NDIFFN * EKINAV (2)
839 END IF
840 ELSE IF ( NINT (ANOW-ZNOW) - NGREYN .LT. 0 ) THEN
841 NDIFFN = NGREYN - NINT ( ANOW - ZNOW )
842 NGREYN = NGREYN - NDIFFN
843 EINCN = EINCN - NDIFFN * EKINAV (2)
844 NGREYP = NGREYP + NDIFFN
845 EINCP = EINCP + NDIFFN * EKINAV (1)
846 IF ( NINT (ZNOW) - NGREYP .LT. 0 ) THEN
847 NDIFFP = NGREYP - NINT ( ZNOW )
848 NGREYP = NGREYP - NDIFFP
849 EINCP = EINCP - NDIFFP * EKINAV (1)
850 END IF
851 END IF
852 NGREYT = NGREYP + NGREYN
853 RETURN
854 END