]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/fluka/nudisv.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 / nudisv.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1995/10/24 10:19:58 cernlib
6* Geant
7*
8*
9#include "geant321/pilot.h"
10*CMZ : 3.21/02 25/07/94 15.04.20 by S.Giani
11*-- Author :
12*$ CREATE NUDISV.FOR
13*COPY NUDISV
14*
15*=== nudisv ===========================================================*
16*
17 INTEGER FUNCTION NUDISV ( ANU, KB, EXPLAM, ASEASQ, APOWER, PRZERO)
18
19#include "geant321/dblprc.inc"
20#include "geant321/dimpar.inc"
21#include "geant321/iounit.inc"
22*
23*----------------------------------------------------------------------*
24* *
25* Last change on 05-may-92 by Alfredo Ferrari, INFN-Milan *
26* *
27* Nudist91: *
28* multiplicity distribution of chains produced in hadron-nucleus *
29* collisions *
30* *
31* Input variables: *
32* anu = average number of collisions *
33* kb = baryon number of the projectile *
34* explam = exp ( - (anu-1) ) *
35* Output variables: *
36* nudisv = number of high energy collisions *
37*----------------------------------------------------------------------*
38*
39 PARAMETER ( INMAX = 30 )
40 PARAMETER ( IAMAX = 13 )
41 PARAMETER ( ANUMAX = 6.0D+00 )
42 PARAMETER ( ANUMED = 0.5D+00 * ( 3.75D+00 + ANUMAX ) )
43*
44C AB,(KB.NE.0) BARYON NUCLEUS
45C
46 DIMENSION AB (INMAX,IAMAX), ANUB (IAMAX), ANUM (9),
47 & ANSESQ (IAMAX), ANPOWR (IAMAX), NULAST (IAMAX)
48 DIMENSION AAB (IAMAX), DIST (0:INMAX), IEX (3,5)
49 REAL RNDM(1)
50 LOGICAL LFIRST
51 SAVE LFIRST, ANUB, ANUM, AB, ANSESQ, NULAST, ANPOWR
52 DATA ANUB/1.D0,1.1D0,1.25D0,1.45D0,1.65D0,1.74D0,2.09D0,2.66D0,
53 & 3.0D0,3.1D0,3.75D0,ANUMED,ANUMAX/
54 DATA ANUM/1.D0,1.45D0,1.52D0,1.77D0,2.19D0,2.42D0,2.47D0,
55 *2.90D0,5.D0/
56 DATA AB / 1.D+00, 29*0.D+00, 90*0.D+00,
57 & 0.5771D+00, 0.2565D+00, 0.1159D+00, 0.0397D+00,
58 & 0.0095D+00, 0.0013D+00, 0.0001D+00, 23*0.D+00,
59 & 0.5457D+00, 0.2548D+00, 0.1303D+00, 0.0515D+00,
60 & 0.0146D+00, 0.0028D+00, 0.0003D+00, 23*0.D+00,
61 & 0.4462D+00, 0.2445D+00, 0.1599D+00, 0.0919D+00, 0.0406D+00,
62 & 0.0132D+00, 0.0032D+00, 0.0005D+00, 0.0001D+00, 21*0.D+00,
63 & 0.3387D+00, 0.2105D+00, 0.1670D+00, 0.1291D+00, 0.0849D+00,
64 & 0.0441D+00, 0.0180D+00, 0.0059D+00, 0.0014D+00, 0.0003D+00,
65 & 0.0001D+00, 19*0.D+00,
66 & 0.2959D+00, 0.1885D+00, 0.1673D+00, 0.1344D+00, 0.1044D+00,
67 & 0.0653D+00, 0.0336D+00, 0.0144D+00, 0.0046D+00, 0.0012D+00,
68 & 0.0003D+00, 19*0.D+00,
69 & 0.2835D+00, 0.1818D+00, 0.1572D+00, 0.1357D+00, 0.1082D+00,
70 & 0.0715D+00, 0.0380D+00, 0.0161D+00, 0.0057D+00, 0.0018D+00,
71 & 0.0004D+00, 0.0001D+00, 18*0.D+00,
72 & 0.2247D+00, 0.1481D+00, 0.1342D+00, 0.1328D+00, 0.1211D+00,
73 & 0.0984D+00, 0.0700D+00, 0.0394D+00, 0.0198D+00, 0.0077D+00,
74 & 0.0028D+00, 0.0008D+00, 0.0002D+00, 17*0.D+00, 60*0.D+00/
75 DATA IEX / 1, 3, 6, 1, 2, 3, 2, 4, 5, 10, 12, 11, 11, 13, 12 /
76 DATA LFIRST / .TRUE. /
77*
78* +-------------------------------------------------------------------*
79* | First call: perform the initialization
80 IF ( LFIRST ) THEN
81 LFIRST = .FALSE.
82 IPOWER = NINT (-APOWER)
83* | +----------------------------------------------------------------*
84* | |
85 DO 100 IA = 1, IAMAX
86 AAB (IA) = 0.D+00
87* | | +-------------------------------------------------------------*
88* | | | This loop computes the normalization factor
89 DO 80 IN = 1, INMAX
90 AAB (IA) = AAB (IA) + AB (IN,IA)
91 80 CONTINUE
92* | | |
93* | | +-------------------------------------------------------------*
94* | | +-------------------------------------------------------------*
95* | | |
96 IF ( AAB (IA) .GT. 0.D+00 ) THEN
97 ANUAVE = 0.D+00
98* | | | +----------------------------------------------------------*
99* | | | |
100 DO 90 IN = 1, INMAX
101 AB (IN,IA) = AB (IN,IA) / AAB (IA)
102 ANUAVE = ANUAVE + IN * AB (IN,IA)
103 90 CONTINUE
104* | | | |
105* | | | +----------------------------------------------------------*
106 ANUB (IA) = ANUAVE
107 END IF
108* | | |
109* | | +-------------------------------------------------------------*
110 100 CONTINUE
111* | |
112* | +----------------------------------------------------------------*
113* | +----------------------------------------------------------------*
114* | |
115 DO 200 IE = 1,5
116 IA = IEX (2,IE)
117 IR1 = IEX (1,IE)
118 IR2 = IEX (3,IE)
119 ANUEN1 = 0.D+00
120 ANUEN2 = 0.D+00
121 AAB (IA) = 0.D+00
122 AEXTR1 = 0.5D+00 * MIN ( 1, IR1 - 1 )
123 AEXTR2 = 0.5D+00 * MIN ( 1, IR2 - 1 )
124* | | +-------------------------------------------------------------*
125* | | | This loop fills the ab(i,ia) array,
126* | | | extrapolating from the data for <nu>s,
127* | | | anub (ir1,r2)
128 DO 130 IN = 1, INMAX
129 ANUBE1 = ANUEN1
130 ANUBE2 = ANUEN2
131 ANUEN1 = ( AEXTR1 + IN ) * ANUB (IA)
132 & / ANUB (IR1)
133 ANUEN2 = ( AEXTR2 + IN ) * ANUB (IA)
134 & / ANUB (IR2)
135 NABE1 = INT (ANUBE1)
136 NABE2 = INT (ANUBE2)
137 NAEN1 = NINT (ANUEN1)
138 NAEN2 = NINT (ANUEN2)
139 IF ( AB (IN,IR1) .GT. 0.D+00 ) THEN
140 DO 110 INN = MAX (NABE1,1), NAEN1
141 IF (INN .GT. 1) THEN
142 ANUBE0 = 0.5D+00 + (INN-1)
143 IF (INN .GT. INMAX ) GO TO 110
144 ELSE
145 ANUBE0 = 0.D+00
146 END IF
147 ANUEN0 = 0.5D+00 + INN
148 ANUBEG = MAX ( ANUBE1, ANUBE0 )
149 ANUEND = MIN ( ANUEN1, ANUEN0 )
150 WEIGHT = AB (IN,IR1) * MAX ( ZERZER, ( ANUEND
151 & - ANUBEG ) / ( ANUEN1 - ANUBE1 ) )
152 AB (INN,IA) = AB (INN,IA) + WEIGHT
153 AAB (IA) = AAB (IA) + WEIGHT
154 110 CONTINUE
155 END IF
156 IF ( AB (IN,IR2) .GT. 0.D+00 ) THEN
157 DO 120 INN = MAX (NABE2,1), NAEN2
158 IF (INN .GT. 1) THEN
159 ANUBE0 = 0.5D+00 + (INN-1)
160 IF (INN .GT. INMAX ) GO TO 120
161 ELSE
162 ANUBE0 = 0.D+00
163 END IF
164 ANUEN0 = 0.5D+00 + INN
165 ANUBEG = MAX ( ANUBE2, ANUBE0 )
166 ANUEND = MIN ( ANUEN2, ANUEN0 )
167 WEIGHT = AB (IN,IR2) * MAX ( ZERZER, ( ANUEND
168 & - ANUBEG ) / ( ANUEN2 - ANUBE2 ) )
169 AB (INN,IA) = AB (INN,IA) + WEIGHT
170 AAB (IA) = AAB (IA) + WEIGHT
171 120 CONTINUE
172 END IF
173 130 CONTINUE
174* | | |
175* | | +-------------------------------------------------------------*
176* | | +-------------------------------------------------------------*
177* | | |
178 DO 140 IN =1, INMAX
179 AB (IN,IA) = AB (IN,IA) / AAB (IA)
180 140 CONTINUE
181* | | |
182* | | +-------------------------------------------------------------*
183 200 CONTINUE
184* | |
185* | +----------------------------------------------------------------*
186* | +----------------------------------------------------------------*
187* | | This loop simply creates a cumulative distribution
188 DO 12 IA = 1, IAMAX
189 AAB (IA) = 0.D+00
190 NULAST (IA) = INMAX
191* | | +-------------------------------------------------------------*
192* | | | This loop computes the normalization factor
193 DO 18 IN = INMAX, 1, -1
194 AAB (IA) = AAB (IA) + AB (IN,IA)
195 IF ( AB (IN,IA) .LE. 0.D+00 ) NULAST (IA) = IN - 1
196 18 CONTINUE
197* | | |
198* | | +-------------------------------------------------------------*
199 ANUAVE = AB (1,IA)
200 ANSESQ (IA) = 0.D+00
201 AB (1,IA) = AB (1,IA) / AAB (IA)
202* | | +-------------------------------------------------------------*
203* | | | Create the cumulative distribution
204 DO 13 IN = 2, INMAX
205 AB (IN,IA) = AB (IN,IA) / AAB (IA) + AB (IN-1,IA)
206 ANUAVE = ANUAVE + IN * ( AB (IN,IA)
207 & - AB (IN-1,IA) )
208 ANSESQ (IA) = ANSESQ (IA) + (IN-1)*(IN-1)
209 & * ( AB (IN,IA) - AB (IN-1,IA) )
210 13 CONTINUE
211* | | |
212* | | +-------------------------------------------------------------*
213 ANUB (IA) = ANUAVE
214* | | +-------------------------------------------------------------*
215* | | |
216 IF ( IPOWER .LT. 7 ) THEN
217 ANPOWR (IA) = AB (1,IA)
218* | | |
219* | | +-------------------------------------------------------------*
220* | | |
221 ELSE IF ( IPOWER .LT. 11 ) THEN
222 ANPOWR (IA) = AB (1,IA) * FPOWER ( IPOWER, 1, ANUAVE )
223* | | |
224* | | +-------------------------------------------------------------*
225* | | |
226 ELSE
227 ANPOWR (IA) = AB (1,IA) * FPOWER ( IPOWER, 1, ANUAVE )
228 END IF
229* | | |
230* | | +-------------------------------------------------------------*
231* | | +-------------------------------------------------------------*
232* | | | Compute < nu**y(nu) >:
233 DO 11 IN = 2, INMAX
234 IF ( IPOWER .LT. 7 ) THEN
235 IF ( APOWER .GT. 0.D+00 ) THEN
236 DPOWER = APOWER
237 ELSE
238 DPOWER = FPOWER ( IPOWER, IN, ANUAVE )
239 END IF
240 ANPOWR (IA) = ANPOWR (IA) + IN**DPOWER
241 & * ( AB (IN,IA) - AB (IN-1,IA) )
242 ELSE IF ( IPOWER .LT. 11 ) THEN
243 ANPOWR (IA) = ANPOWR (IA) + IN
244 & * FPOWER ( IPOWER, IN, ANUAVE )
245 & * ( AB (IN,IA) - AB (IN-1,IA) )
246 ELSE
247 ANPOWR (IA) = ANPOWR (IA) + IN
248 & **FPOWER ( IPOWER, IN, ANUAVE )
249 & * ( AB (IN,IA) - AB (IN-1,IA) )
250 END IF
251 11 CONTINUE
252* | | |
253* | | +-------------------------------------------------------------*
254 12 CONTINUE
255* | |
256* | +----------------------------------------------------------------*
257* | +----------------------------------------------------------------*
258* | |
259 DO 17 IA = 1, IAMAX
260 ANUAC = AB (1,IA)
261 ANUSQ = AB (1,IA)
262 ANUTH = AB (1,IA)
263* | | +-------------------------------------------------------------*
264* | | |
265 DO 15 I = 2, NULAST (IA)
266 ANUAC = ANUAC + I * ( AB (I,IA) - AB(I-1,IA) )
267 ANUSQ = ANUSQ + I*I * ( AB (I,IA) - AB(I-1,IA) )
268 ANUTH = ANUTH + I*I*I * ( AB (I,IA) - AB(I-1,IA) )
269 15 CONTINUE
270* | | |
271* | | +-------------------------------------------------------------*
272* | | +-------------------------------------------------------------*
273* | | |
274 IF ( IA .GT. 1 ) THEN
275 ANUSE2 = ( ANUSQ - 2.D+00 * ANUB (IA) + 1.D+00 ) /
276 & ( ANUB (IA) - 1.D+00 )**2
277 ANUSE3 = ( ANUTH - 3.D+00 * ANUSQ + 3.D+00 * ANUB (IA)
278 & - 1.D+00 ) / ( ANUB (IA) - 1.D+00 )**3
279* | | |
280* | | +-------------------------------------------------------------*
281* | | |
282 ELSE
283 ANUSE2 = 1.D+00
284 ANUSE3 = 1.D+00
285 END IF
286* | | |
287* | | +-------------------------------------------------------------*
288 17 CONTINUE
289* | |
290* | +----------------------------------------------------------------*
291 NUDISV = -1
292 RETURN
293 END IF
294* |
295* +-------------------------------------------------------------------*
296* +-------------------------------------------------------------------*
297* | Select the proper distributions for interpolations
298 DO 40 I = 1, IAMAX
299 IF (ANU .LT. ANUB(I)) GO TO 41
300 40 CONTINUE
301* |
302* +-------------------------------------------------------------------*
303 I=IAMAX + 1
304 41 CONTINUE
305 NANU = I - 1
306 IF (NANU.GE.IAMAX) GO TO 50
307 IF (NANU.LE.0) NANU = 1
308 NANUN = NANU + 1
309 WEIGH1 = ( ANU - ANUB (NANU) ) / ( ANUB (NANUN) - ANUB (NANU) )
310 WEIGH2 = ( ANUB (NANUN) - ANU ) / ( ANUB (NANUN) - ANUB (NANU) )
311 DIST (0) = 0.D+00
312 ASEASQ = WEIGH1 * ANSESQ (NANUN) + WEIGH2 * ANSESQ (NANU)
313 APOWER = WEIGH1 * ANPOWR (NANUN) + WEIGH2 * ANPOWR (NANU)
314* +-------------------------------------------------------------------*
315* | Create the proper cumulative distribution by interpolation
316* | ( note that since weigh1 + weigh2 = 1 it will be automatically
317* | normalized )
318 DO 20 IN = 1, NULAST (NANUN)
319 DIST (IN) = WEIGH1 * AB (IN,NANUN) + WEIGH2 * AB (IN,NANU)
320 WEIGHT = (IN-1) * ( DIST (IN) - DIST (IN-1) )
321 20 CONTINUE
322* |
323* +-------------------------------------------------------------------*
324 PRZERO = DIST (1)
325 CALL GRNDM(RNDM,1)
326 X=RNDM(1)
327* +-------------------------------------------------------------------*
328* | Compute the proper nu from the cumulative distribution
329 DO 30 IN = 1, NULAST (NANUN)
330 IF ( X .LE. DIST (IN) ) GO TO 31
331 30 CONTINUE
332* |
333* +-------------------------------------------------------------------*
334 IN = NULAST (NANUN)
335 31 CONTINUE
336 NUDISV = IN
337 RETURN
338* Come here for <nu> larger than 8
339 50 CONTINUE
340 APOWER = -1.D+00
341 ASEASQ = -1.D+00
342 NUD = KPOIS (EXPLAM)
343 NUDISV = NUD + 1
344 PRZERO = EXPLAM
345 WRITE (LUNERR,*)' *** Nudisv called with <nu> >= 8 !!',REAL(ANU),
346 & ' ***'
347 RETURN
348*=== end of function Nudisv ===========================================*
349 END