]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/fluka/corrin.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / GEANT321 / fluka / corrin.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.42 by S.Giani
11*-- Author :
12*$ CREATE CORRIN.FOR
13*COPY CORRIN
14* *
15*=== corrin ===========================================================*
16* *
17 SUBROUTINE CORRIN ( 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* Created on 10 may 1990 by Alfredo Ferrari & Paola Sala *
26* Infn - Milan *
27* *
28* Last change on 23-mar-93 by Alfredo Ferrari *
29* *
30* This version has been developed by A. Ferrari starting from the *
31* the one by J. M. Zazula: it includes a few differences mainly in *
32* the correlation since now we try to get the excitation energy *
33* from a correct energy balance *
34*----------------------------------------------------------------------*
35*
36#include "geant321/balanc.inc"
37#include "geant321/corinc.inc"
38#include "geant321/nucdat.inc"
39#include "geant321/parevt.inc"
40#include "geant321/part.inc"
41#include "geant321/resnuc.inc"
42 PARAMETER ( ALFMAX = 5.00D+00 )
43 PARAMETER ( FINCFR = 1.00D+00 )
44 PARAMETER ( TMNOLD = 0.02D+00 )
45 COMMON / FKNUCO / HELP (2), HHLP (2), FTVTH (2), FINCX (2),
46 & EKPOLD (2), BBOLD, ZZOLD, SQROLD, ASEASQ,
47 & FSPRED, FEX0RD
48 DIMENSION KPA (26)
49 DIMENSION FRA(2,110), ESLOLD (2), AVMULT (2), EXDUMM (2),
50 & V0WSAV (2), VEWSAV (2)
51 REAL RNDM(3)
52 LOGICAL LSTOP, LNUCNW, LLLIN, LLPOW
53 SAVE FRA, FRAC, FRAC0, FRAMAX, KPA, IJJ, LSTOP, V0WSAV, VEWSAV
54 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/
55*
56* Reduction factors for intran. cascade energy, taken from Alsmiller
57* Incoming baryons
58 DATA (FRA(1,I), I=1,107) / .048D0,.076D0,0.10D0,.12D0,
59 * .14D0,.16D0,.17D0,.19D0,.21D0,.22D0,.24D0,
60 1 .25D0,.26D0,.28D0,.29D0,.30D0,.315D0,.33D0,.34D0,
61 * .35D0,.36D0,.38D0,.39D0,.40D0,.41D0,.42D0,
62 2 .43D0,.44D0,.46D0,.47D0,.48D0,.49D0,.50D0,.51D0,
63 * .52D0,.53D0,.54D0,.55D0,.56D0,.57D0,.58D0,.59D0,
64 3 .60D0,.61D0,.62D0,.63D0,.635D0,.64D0,.65D0,.66D0,
65 * .67D0,.675D0,.68D0,.69D0,.70D0,.71D0,.715D0,
66 4 .72D0,.725D0,.73D0,.735D0,.74D0,.75D0,.755D0,
67 * .76D0,.767D0,.77D0,.77D0,.775D0,.78D0,.783D0,
68 5 .786D0,.79D0,.795D0,.80D0,.805D0,.81D0,.812D0,
69 * .815D0,.82D0,.822D0,.824D0,.825D0,.825D0,.83D0,
70 6 .832D0,.834D0,.836D0,.838D0,.84D0,.843D0,.836D0,
71 * .849D0,.852D0,.855D0,.856D0,.857D0,.858D0,
72 7 .859D0,.860D0,.862D0,.864D0,.866D0,.868D0,.870D0,
73 * .872D0,.874D0/
74* Incoming mesons
75 DATA (FRA(2,I), I=1,63) / .048D0,.076D0,0.10D0,.12D0,
76 * .14D0,.16D0,.17D0,.19D0,.21D0,.22D0,.24D0,
77 1 .25D0,.262D0,.274D0,.285D0,.295D0,.305D0,.315D0,
78 * .327D0,.340D0,.345D0,.350D0,.360D0,
79 2 .370D0,.380D0,.385D0,.390D0,.398D0,.406D0,.415D0,
80 * .420D0,.425D0,.430D0,.435D0,.440D0,.446D0,
81 3 .452D0,.458D0,.464D0,.470D0,.474D0,.478D0,.482D0,
82 * .486D0,.490D0,.495D0,.500D0,.505D0,.510D0,
83 4 .515D0,.519D0,.523D0,.526D0,.520D0,.533D0,.537D0,
84 * .540D0,.543D0,.547D0,.550D0,.555D0,.559D0,
85 5 .5625D0 /
86*
87 IBPROJ = IBAR (KPROJ)
88 IJJ = KPA (KPROJ)
89 ANUAV = 1.D+00
90 ANCOLL = 1.D+00
91 ANUSEA = 1.D+00
92 NSEA = 0
93 LSTOP = .FALSE.
94* +-------------------------------------------------------------------*
95* | Incoming baryons
96 IF ( IBPROJ .NE. 0 ) THEN
97* | +----------------------------------------------------------------*
98* | |
99 IF ( IBTAR .LE. 107 ) THEN
100 FRAC = FRA ( 1, IBTAR )
101* | |
102* | +----------------------------------------------------------------*
103* | |
104 ELSE IF (IBTAR .LE. 206) THEN
105 FRAC = 0.739D+00 + 0.00126D+00 * BBTAR
106* | |
107* | +----------------------------------------------------------------*
108* | |
109 ELSE
110 FRAC = 1.D+00
111 END IF
112* | |
113* | +----------------------------------------------------------------*
114* |
115* +-------------------------------------------------------------------*
116* | Incoming mesons
117 ELSE
118* | +----------------------------------------------------------------*
119* | |
120 IF ( IBTAR .LE. 63 ) THEN
121 FRAC = FRA ( 2, IBTAR ) * ( 1.D+00 + 0.1333D+00 * MAX (
122 & 0.D+00, BBTAR - 35.D+00 ) / 28.D+00 )
123* | |
124* | +----------------------------------------------------------------*
125* | |
126 ELSE IF ( IBTAR .LE. 107 ) THEN
127 FRAC = 0.85D+00 * FRA ( 1, IBTAR )
128* | |
129* | +----------------------------------------------------------------*
130* | |
131 ELSE IF ( IBTAR .LE. 206 ) THEN
132 FRAC = 0.6279D+00 + 0.001077D+00 * BBTAR
133* | |
134* | +----------------------------------------------------------------*
135* | |
136 ELSE
137 FRAC = 0.85D+00
138 END IF
139* | |
140* | +----------------------------------------------------------------*
141 END IF
142* |
143* +-------------------------------------------------------------------*
144* +-------------------------------------------------------------------*
145* |
146 IF ( BBTAR .NE. BBOLD .OR. ZZTAR .NE. ZZOLD ) THEN
147 LNUCNW = .TRUE.
148* | Supply the fraction of the total kinetic energy to be
149* | used for intranuclear cascade nucleons
150 SQRAMS = SQRT ( BBTAR )
151 ATO1O3 = BBTAR**0.3333333333333333D+00
152* ZTO1O3 = ZZTAR**0.3333333333333333D+00
153 HKAP = BBTAR**2 / ( ZZTAR**2 + ( BBTAR - ZZTAR )**2 )
154 HHLP (1) = ( HKAP * ZZTAR )**0.3333333333333333D+00 / ATO1O3
155 HHLP (2) = ( HKAP * ( BBTAR - ZZTAR ) )
156 & **0.3333333333333333D+00 / ATO1O3
157 RDSNUC = R0NUCL * ATO1O3
158 RDSCOU = RCCOUL * ATO1O3
159 FLKCOU = DOST ( 1, ZZTAR )
160 VEFFNU (1) = COULPR * FLKCOU * ZZTAR / RDSCOU
161 VEFFNU (2) = 0.D+00
162 AMRCAV = MAX ( 1.D+00, ( BBTAR - 2.D+00 ) ) * AMUC12
163 AMRCSQ = AMRCAV * AMRCAV
164* | +----------------------------------------------------------------*
165* | |
166 DO 3000 I = 1, 2
167 PFRMMX (I) = HHLP (I) * APFRMX
168 P2HELP = PFRMMX (I)**2
169 EFRMMX (I) = SQRT ( P2HELP + AMNUSQ (I) ) - AMNUCL (I)
170 EFRMAV (I) = 0.3D+00 * P2HELP / AMNUCL (I) * ( 1.D+00
171 & - P2HELP / ( 5.6D+00 * AMNUSQ (I) ) )
172 ERCLAV (I) = 0.3D+00 * P2HELP / AMRCAV * ( 1.D+00
173 & - P2HELP / ( 5.6D+00 * AMRCSQ ) )
174 V0WELL (I) = EFRMMX (I) + EBNDNG (I)
175 VEFFNU (I) = VEFFNU (I) + V0WELL (I)
176 V0WSAV (I) = V0WELL (I)
177 VEWSAV (I) = VEFFNU (I)
178 ERCLMX = 0.5D+00 * P2HELP / AMRCAV * ( 1.D+00
179 & - 0.25D+00 * P2HELP / AMRCSQ )
180 EKMNNU (I) = VEFFNU (I) + EBNDNG (I)
181 & - EFRMMX (I) + ERCLMX
182 EKMXNU (I) = VEFFNU (I) + EBNDNG (I)
183 EKMNAV (I) = VEFFNU (I) + EBNDNG (I)
184 & - EFRMAV (I) + ERCLAV (I)
185 ESWELL (I) = EBNDNG (I) + V0WELL (I)
186 & - EFRMAV (I) + ERCLAV (I)
187 FTVTH (I) = ( EKMNNU (I) / EFRMMX (I) )**3
188 FTVTH (I) = 0.4D+00 * SQRT ( FTVTH (I) )
1893000 CONTINUE
190* | |
191* | +----------------------------------------------------------------*
192 BBOLD = BBTAR
193 ZZOLD = ZZTAR
194 SQROLD = SQRAMS
195* |
196* +-------------------------------------------------------------------*
197* |
198 ELSE
199 V0WSAV (1) = V0WELL (1)
200 VEWSAV (1) = VEFFNU (1)
201 V0WSAV (2) = V0WELL (2)
202 VEWSAV (2) = VEFFNU (2)
203 LNUCNW = .FALSE.
204 END IF
205* |
206* +-------------------------------------------------------------------*
207 FRAC0 = FRAC
208 RETURN
209 ENTRY CORSTP ( EKEFF )
210 LSTOP = .TRUE.
211 ENTRY CORRNC ( EKEFF )
212 WEIGH1 = MIN ( 1.D+00, EKEFF / 4.D+00 )
213 FRAC = ( WEIGH1 * 2.5D+00 + ( 1.D+00 - WEIGH1 ) ) * FRAC0
214 FRAMAX = 2.5D+00
215* The following is a reduction factor to bring Hannes's parametri-
216* zations for the average energy carried by cascade nucleons in
217* better agreement with experimental data
218 FSPRED = FSPRD0 + ( 1.D+00 - FSPRD0 ) * MAX ( 8.D+00 - EKEFF,
219 & 0.D+00 ) / 8.D+00
220 FSPRED = FSPRED * MIN ( EKEFF / 0.8D+00, 1.D+00 )
221 TMPFSP = 1.25D+00 * FSPRD0
222 FSPRED = MIN ( FSPRED, TMPFSP )
223 AVMULT (1) = FRAC * BNKEKA ( 1, EKEFF, BBOLD, SQROLD )
224 AVMULT (2) = FRAC * BNKEKA ( 2, EKEFF, BBOLD, SQROLD )
225 EKUPNU (1) = BEKEKA ( 2, EKEFF, BBOLD, SQROLD )
226 EKUPNU (2) = BEKEKA ( 3, EKEFF, BBOLD, SQROLD )
227 EINCP = FRAC * EKUPNU (1)
228 EINCN = FRAC * EKUPNU (2)
229 EKMAX = EKEFF - EBNDAV
230 EKUPN0 = MAX ( EKUPNU (1), EKUPNU (2) )
231* +-------------------------------------------------------------------*
232* |
233 DO 4000 I = 1, 2
234 FINCUP (I) = BKEKA ( I, EKEFF, BBOLD )
235 RATRAT = MIN ( 1.D+00, 0.5D+00 * EKMAX / FINCUP (I) )
236 TMPFIN = 0.5D+00 * EKMAX
237 FINCUP (I) = MIN ( FINCUP (I), TMPFIN )
238 EKPOLD (I) = EKUPNU (I) * FRAC * RATRAT
239 EKUPNU (I) = FINCUP (I)
240 TMPEKU = 1.5D+00 * EKMXNU (I)
241 EKUPNU (I) = MAX ( FRAMAX * EKUPNU (I), TMPEKU )
242 TMPEKU = 0.7D+00 * EKEFF
243 EKUPNU (I) = MIN ( EKUPNU (I), TMPEKU )
244 ESLOLD (I) = MAX ( FINCUP (I), TMNOLD )
245 FINCX (I) = ( ESLOLD (I) + ESWELL (I) ) / ESLOLD (I) * FSPRED
246 ESLOPE (I) = MIN ( ESLOLD (I) * FINCX (I), EKMAX )
247 AHELP = EXP ( - EKPOLD (I) / ESLOLD (I) )
248 EKNOLD = ESLOLD (I) - EKPOLD (I) * AHELP /
249 & ( 1.D+00 - AHELP )
250 EXMNAV (I) = EXP ( - EKMNAV (I) / ESLOPE (I) )
251 EXMNNU (I) = EXP ( - EKMNNU (I) / ESLOPE (I) )
252 EXUPNU (I) = EXP ( - EKUPNU (I) / ESLOPE (I) )
253 EKINAV (I) = ESLOPE (I) + ( EKMNAV (I) * EXMNAV (I) -
254 & EKUPNU (I) * EXUPNU (I) ) / ( EXMNAV (I) -
255 & EXUPNU (I) )
256 FINCUP (I) = AVMULT (I) * EKINAV (I) / EKPOLD (I)
257* FINCUP (I) = EKINAV (I) / EKNOLD
2584000 CONTINUE
259* |
260* +-------------------------------------------------------------------*
261 EINCP = EKPOLD (1)
262 EINCN = EKPOLD (2)
263 IF ( LSTOP ) RETURN
264* Power sampling
265 LLPOW = .TRUE.
266 LLLIN = .FALSE.
267* Take into account that a fraction 1/A^2/3 is lost because of
268* peripheral collisions
269 PERCOR = ATO1O3 * ATO1O3
270 PERCOR = PERCOR / ( PERCOR - 1.D+00 )
271 EINCP0 = EINCP * FINCUP (1) * PERCOR
272 EINCN0 = EINCN * FINCUP (2) * PERCOR
273 EINCT = EINCP0 + EINCN0
274* For the other distr,:
275 EINCMX = EKEFF - EBNDAV - ( AV0WEL - AEFRMA ) / ( 1.D+00 +
276 & 4.D+00 * ( AV0WEL - AEFRMA ) / EKEFF )
277 EIUSE = 0.5D+00 * ( EINCMX + EKMAX )
278 AHNORM = 0.D+00
279 EINCM0 = - AINFNT
280* +-------------------------------------------------------------------*
281* | Compute alfa:
282 IF ( EINCMX .GT. 2.1D+00 * EINCT ) THEN
283 EINFIN = MIN ( 0.95D+00 * EINCMX, ( ALFMAX + 2.D+00 ) * EINCT )
284 ALFA = EINFIN / EINCT - 2.D+00
285* |
286* +-------------------------------------------------------------------*
287* | Energy interval too small, sample E uniformly
288 ELSE
289 LLPOW = .FALSE.
290 END IF
291* |
292* +-------------------------------------------------------------------*
293 EINCMN = 0.D+00
294* +-------------------------------------------------------------------*
295* | Energy is so small that we cannot produce cascade nucleons,
296* | sample an excitation energy and return
297 IF ( EINCMX .LE. EINCMN + 0.25D+00 * EBNDAV ) THEN
298 ISAMPL = 50
299* |
300* +-------------------------------------------------------------------*
301* |
302 ELSE
303 ISAMPL = 0
304 END IF
305* |
306* +-------------------------------------------------------------------*
307* +-------------------------------------------------------------------*
308* | Now make the linear sampling for Eincp, Eincn, on [0,Eicnmx]
309* | according to P(E)dE=A(Efin-E)^a dE
3105000 CONTINUE
311 ISAMPL = ISAMPL + 1
312 CALL GRNDM(RNDM,1)
313 IF ( LLPOW ) THEN
314 EIHELP = EINFIN * ( 1.D+00 - RNDM (1)**(1.D+00
315 & /(ALFA+1.D+00)) )
316 ELSE IF ( LLLIN ) THEN
317 EINCMX = EINCM0
318 EIHELP = EINFIN - SQRT ( EINFIN**2 - 2.D+00 * RNDM (1)
319 & / AHNORM )
320 ELSE
321 EIHELP = EINCMX * RNDM (1)
322 END IF
323* | *** end linear sampling ***
324 EINCP = EINCP0 * EIHELP / EINCT
325 EINCN = EINCN0 * EIHELP / EINCT
326 AHELP = EKMNNU (1) / ESLOPE (1)
327 AHELP = AHELP * FTVTH (1) * EINCP * ( 1.D+00 + AHELP *
328 & 0.3333333333333333D+00 ) / ( 1.D+00 - EXUPNU (1) )
329 TVGRE0 = AHELP
330 AHELP = EKMNNU (2) / ESLOPE (2)
331 AHELP = AHELP * FTVTH (2) * EINCN * ( 1.D+00 + AHELP *
332 & 0.3333333333333333D+00 ) / ( 1.D+00 - EXUPNU (2) )
333 TVGRE0 = TVGRE0 + AHELP
334* | +----------------------------------------------------------------*
335* | | Energy is so small that we cannot produce cascade nucleons,
336* | | sample an excitation energy and return
337 IF ( ISAMPL .GE. 50 ) THEN
338 CALL GRNDM(RNDM,3)
339 PRNDM = MAX ( RNDM (1), RNDM (2), RNDM (3) )
340 TVGRE0 = MAX ( AV0WEL - PRNDM**2 * AEFRMX - EBNDAV, ZERZER )
341 EINCP = 0.D+00
342 EINCN = 0.D+00
343 TVGREY = 0.D+00
344 RETURN
345 END IF
346* | |
347* | +----------------------------------------------------------------*
348* | Changed:
349 IF ( EIHELP + TVGRE0 .GE. EIUSE ) GO TO 5000
350* |
351* +--<--<--<--<--< Resampling!
352* |
353* +-------------------------------------------------------------------*
354 PCR = EIHELP / EINCT
355 FRAINC = FRAC * PCR
356 TVGREY = 0.D+00
357 AGREYP = EINCP / EKINAV (1)
358 AGREYN = EINCN / EKINAV (2)
359* ==== Discretize the distribution !!! ==== *
360 CALL GRNDM(RNDM,2)
361 NGREYP = INT ( AGREYP )
362 IF ( RNDM(1) .LT. AGREYP - NGREYP ) NGREYP = NGREYP + 1
363 NGREYN = INT ( AGREYN )
364 IF ( RNDM(2) .LT. AGREYN - NGREYN ) NGREYN = NGREYN + 1
365* ====
366 AGREYT = AGREYP + AGREYN
367 NGREYT = NGREYP + NGREYN
368 PPROCS = AGREYP / AGREYT
369 NGREYP = 0
370 NGREYN = 0
371 EINCP = 0.D+00
372 EINCN = 0.D+00
373 IF ( NGREYT .LE. 0 ) GO TO 9000
374 ARDUMM = MAX ( 1.D+00, BBOLD - 2.D+00 )
375 CALL RBKINI ( 1, .TRUE., EXDUMM, TKDUMM, TSDUMM,
376 & PSDUMM, ARDUMM, TRDUMM )
377 V0SAV1 = V0WELL (1)
378 V0SAV2 = V0WELL (2)
379 VESAV1 = VEFFNU (1)
380 VESAV2 = VEFFNU (2)
381 V0WELL (1) = V0WSAV (1)
382 V0WELL (2) = V0WSAV (2)
383 VEFFNU (1) = VEWSAV (1)
384 VEFFNU (2) = VEWSAV (2)
385 IRETRY = 0
386 DO 8000 I = 1, NGREYT
387 CALL GRNDM(RNDM,1)
388 RNDMPR = RNDM (1)
389 ARDUMM = MAX ( 1.D+00, ARDUMM - 1.D+00 )
390 7500 CONTINUE
391 IF ( RNDMPR .LT. PPROCS ) THEN
392 NGREYP = NGREYP + 1
393 CALL RBKINI ( 1, .FALSE., EXDUMM, TKDUMM, TSDUMM,
394 & PSDUMM, ARDUMM, TRDUMM )
395 EINCP = EINCP + TKDUMM
396 IF ( EINCP + EINCN .GT. EINCMX .AND. IRETRY .LT. 5 ) THEN
397 EINCP = EINCP - TKDUMM
398 CALL RBKMIN (1)
399 IRETRY = IRETRY + 1
400 GO TO 7500
401 END IF
402 ELSE
403 NGREYN = NGREYN + 1
404 CALL RBKINI ( 2, .FALSE., EXDUMM, TKDUMM, TSDUMM,
405 & PSDUMM, ARDUMM, TRDUMM )
406 EINCN = EINCN + TKDUMM
407 IF ( EINCP + EINCN .GT. EINCMX .AND. IRETRY .LT. 5 ) THEN
408 EINCN = EINCN - TKDUMM
409 CALL RBKMIN (2)
410 IRETRY = IRETRY + 1
411 GO TO 7500
412 END IF
413 END IF
414 8000 CONTINUE
415 V0WELL (1) = V0SAV1
416 V0WELL (2) = V0SAV2
417 VEFFNU (1) = VESAV1
418 VEFFNU (2) = VESAV2
419 9000 CONTINUE
420 EINCT = EINCP + EINCN
421 IF ( EINCT + TVGRE0 .GT. EIUSE ) THEN
422 TVGRE0 = 0.D+00
423 EIUSE = MIN ( ONEONE, EINCMX / EINCT )
424 EINCP = EINCP * EIUSE
425 EINCN = EINCN * EIUSE
426 EINCT = EINCP + EINCN
427 END IF
428 RETURN
429*=== End of subroutine corrin =========================================*
430 END