Use double precision constants
[u/mrichter/AliRoot.git] / GEANT321 / peanut / bimsel.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
9e58b465 5* Revision 1.5 2002/06/25 08:17:33 hristov
6* Additional protection
7*
0ba2e3bb 8* Revision 1.4 2002/06/22 10:50:18 hristov
9* Better protection
10*
9c4c8b17 11* Revision 1.3 2002/06/20 15:58:37 hristov
12* Protection in case of invalid ACOS argument. The reason for such argument has to be investigated.
13*
900dc090 14* Revision 1.2 2002/05/13 12:40:58 hristov
15* Dummy subroutines to avoid files with no code in
16*
62be6b28 17* Revision 1.1.1.1 1999/05/18 15:55:22 fca
18* AliRoot sources
19*
fe4da5cc 20* Revision 1.2 1997/10/17 10:00:03 mclareni
21* Remove SAVE statement for NT
22*
23* Revision 1.1.1.1 1995/10/24 10:22:00 cernlib
24* Geant
25*
26*
27#include "geant321/pilot.h"
28*CMZ : 3.21/02 29/03/94 15.41.45 by S.Giani
29*-- Author :
30*$ CREATE BIMSEL.FOR
31*COPY BIMSEL
32*
33*=== bimsel ===========================================================*
34*
35 SUBROUTINE BIMSEL ( JPROJ, TXX, TYY, TZZ, LBCHCK )
36
37#include "geant321/dblprc.inc"
38#include "geant321/dimpar.inc"
39#include "geant321/iounit.inc"
40*
41*----------------------------------------------------------------------*
42*----------------------------------------------------------------------*
43*
44#include "geant321/balanc.inc"
45#include "geant321/eva0.inc"
46#include "geant321/nucdat.inc"
47#include "geant321/nucgeo.inc"
48#include "geant321/part.inc"
49#include "geant321/parevt.inc"
50#include "geant321/resnuc.inc"
51*
52 PARAMETER ( FEFFEC = 1.518066780142162 D+00 )
53 PARAMETER ( BETMAX = 0.4 D+00 )
54*
55 REAL RNDM(2)
56*
57 LOGICAL LBCHCK, LFERMI, LLMDBR
58*
59 SAVE ABTAR , ZZTAR , SIGMP0, SIGMN0,
60 & AMNHLP, RHOBIM, RPRONU, RADPRP, RADPRN, DSKRED, RHRUSF,
61 & AUSFL , ZUSFL , BNDSAV, RADHLP, BFCHLP, BIMCLM, PRCOLP,
62 & PRCOLN, IBTOLD, ICTOLD, KPROJ , NTRIAL, ITFRMI
63 SAVE SBHAL0, SBSKI0, SBCEN0, SBCEN1, SBSKI1, SBHAL1
64#ifndef CERNLIB_WINNT
62be6b28 65* SAVE
fe4da5cc 66#endif
67 DATA IBTOLD, ICTOLD / 2*0 /
68*
69 KPROJ = JPROJ
70 AUSFL = IBTAR
71 ZUSFL = ICHTAR
72 RHRUSF = 1.D+00
73 BEPROJ = PNUCCO / ( EKECON + AM (KPROJ) )
74 CXIMPC = TXX
75 CYIMPC = TYY
76 CZIMPC = TZZ
77 NTRIAL = 0
78 RHOBIM = - AINFNT
79 IF ( KPROJ .EQ. 1 .OR. KPROJ .EQ. 8 ) THEN
80 IPWELL = 1 + KPROJ / 8
81 WLLRED = 1.D+00
82 BNDNUC = BNENRG (IPWELL)
83 ELSE
84 IPWELL = 0
85 IF ( IBAR (KPROJ) .EQ. 0 ) THEN
86 IF ( KPROJ .LE. 11 ) THEN
87 WLLRED = 0.D+00
88 BNDNUC = 0.D+00
89 ELSE
90 WLLRED = POTMES
91 BNDNUC = BNENRG (3)
92 END IF
93 ELSE
94 WLLRED = POTBAR
95 BNDNUC = BNENRG (3)
96 END IF
97 END IF
98 BNDSAV = BNDNUC
99 IF ( IBAR (KPROJ) .NE. 0 ) THEN
100 RPRONU = 1.D+00
101 ELSE IF ( KPROJ .NE. 7 ) THEN
102 RPRONU = 0.8164965809277260D+00
103 ELSE
104 RPRONU = 0.D+00
105 END IF
106 IF ( LBCHCK ) THEN
107 LFERMI = .FALSE.
108 EKESIG = EKECON
109 PPRSIG = PNUCCO
110 CALL SIGFER ( KPTOIP (KPROJ), EKESIG, PPRSIG, LFERMI )
111 PRCOLP = ZUSFL / AUSFL * SIGMAP
112 PRCOLN = ( AUSFL - ZUSFL ) / AUSFL * SIGMAN
113 SIGMAA = PRCOLP + PRCOLN
114 PRCOLP = PRCOLP / SIGMAA
115 PRCOLN = 1.D+00 - PRCOLP
116 END IF
117 IF ( RPRONU .GT. ANGLGB ) THEN
118 IF ( LPARWV ) THEN
119 LLMDBR = .TRUE.
120 TMP102 = 1.D-02
121 PDEBRO = MAX ( PNUCCO, TMP102 )
122 ALMBAR = PLABRC / PDEBRO
123 DEBRLM = 0.5D+00 * ALMBAR
124 RADCOR = SQRT ( (RPRONU * RMSPRO)**2 + ALMBAR**2 )
125 & / ( RMSPRO * RPRONU )
126 ELSE
127 PDEBRO = ( EKECON + BNDNUC ) * ( EKECON + BNDNUC + 2.D+00
128 & * AM (KPROJ) )
129 LLMDBR = .FALSE.
130 DEBRLM = 0.D+00
131 ALMBAR = 0.D+00
132 LLLMAX = -1
133 RADCOR = SQRT ( (RPRONU * RMSPRO)**2 + PLABRC**2 / PDEBRO
134 & ) / ( RMSPRO * RPRONU )
135 END IF
136 ELSE
137 RADCOR = 0.D+00
138 LLMDBR = .FALSE.
139 DEBRLM = 0.D+00
140 END IF
141 RADCO2 = RADCOR
142 RADPRO = MIN ( TWOTWO * RMSPRO * RPRONU * RADCOR, SKGT16 )
143 RADPRP = RADPRO
144 RADPRN = RADPRO
145 IF ( IBTAR .NE. IBTOLD .OR. ICHTAR .NE. ICTOLD ) THEN
146 IBTOLD = IBTAR
147 ICTOLD = ICHTAR
148 ABTAR = IBTAR
149 ZZTAR = ICHTAR
150 AR1O3 = RMASS (IBTAR)
151 AMNHLP = 0.5D+00 * ( AMNUCL (1) + AMNUCL (2) )
152 HKAP = ABTAR**2 / ( ZZTAR**2 + ( ABTAR - ZZTAR )**2 )
153 HHLP (1) = ( HKAP * ZZTAR )**0.3333333333333333D+00 / AR1O3
154 HHLP (2) = ( HKAP * ( ABTAR - ZZTAR ) )
155 & **0.3333333333333333D+00 / AR1O3
156 RHOCEN = RHOTAB (IBTAR)
157 ALPHAL = ALPTAB (IBTAR)
158 RADIU0 = RADTAB (IBTAR)
159 SKINDP = SKITAB (IBTAR)
160 HALODP = HALTAB (IBTAR)
161 RADIU1 = RADIU0 + SKINDP
162 RADTOT = RADIU1 + HALODP
163 RHOCOR = ONEMNS * RHOCEN
164 RHOSKN = ALPHAL * RHOCEN
165 PFRCEN (1) = HHLP (1) * PFRTAB (IBTAR)
166 PFRCEN (2) = HHLP (2) * PFRTAB (IBTAR)
167 RHOAVE = RHATAB (IBTAR)
168 PFRAVE = PFATAB (IBTAR)
169 EKFAVE = EKATAB (IBTAR)
170 OMALHL = 1.D+00 - ALPHAL
171 RAD1O2 = RADIU0 + 0.5D+00 * SKINDP / OMALHL
172 SKNEFF = SKINDP / OMALHL
173 RADSKN = RADIU0 + SKNEFF
174 EKFCEN (1) = SQRT ( AMNUSQ (1) + PFRCEN (1)**2 ) - AMNUCL (1)
175 EKFCEN (2) = SQRT ( AMNUSQ (2) + PFRCEN (2)**2 ) - AMNUCL (2)
176 IF ( PFRCEN (1) .GT. PFRCEN (2) ) THEN
177 ITNCMX = 1
178 ELSE
179 ITNCMX = 2
180 END IF
181 CALL NCLVST ( IBTAR, ICHTAR )
182 END IF
183 IF ( IPWELL .LE. 0 ) IPWELL = ITNCMX
184 CALL NCLVIN
185 IF ( EKECON .LT. 2.D+00 * GAMMIN ) THEN
186 EKECON = 0.D+00
187 PNUCCO = 0.D+00
188 LABRST = .TRUE.
189 RADPRO = 0.D+00
190 RADSIG = ( RADTOT + DEBRLM ) * BFCLMB
191 RADMAX = RADTOT
192 LLLMAX = -1
193 OPACTY = 2.D+00
194 CALL RSTSEL (KPROJ)
195 RETURN
196 END IF
197 RADMAX = RADTOT + RADPRO
198 BIMCLM = RDCLMB * BFCLMB
199 IF ( LLMDBR ) THEN
200 RADHLP = RADMAX
201 IF ( RADHLP .LE. RDCLMB ) THEN
202 BFCMAX = BFCLMB
203 BFCHLP = 0.5D+00 * CLMBBR * RDCLMB / EKECON
204 ELSE
205 BFCHLP = 0.5D+00 * CLMBBR * RDCLMB / EKECON
206 BFCMAX = SQRT ( 1.D+00 - CLMBBR * RDCLMB / EKECON / RADHLP )
207 END IF
208 BIMMAX = RADHLP * BFCMAX
209 LLLMAX = INT ( BIMMAX / ALMBAR )
210 RADSIG = ALMBAR * ( LLLMAX + 1.D+00 )
211 SIGGEO = PI * RADSIG * RADSIG
212 ELSE
213 RADHLP = RADTOT + RADPRO + DEBRLM
214 IF ( RADHLP .LE. RDCLMB ) THEN
215 BFCMAX = BFCLMB
216 ELSE
217 BFCHLP = 0.5D+00 * CLMBBR * RDCLMB / EKECON
218 BFCMAX = SQRT ( 1.D+00 - CLMBBR * RDCLMB / EKECON / RADHLP )
219 END IF
220 RADSIG = RADHLP * BFCMAX
221 END IF
222 R0TRAJ = - RADTOT
223 R1TRAJ = - R0TRAJ
224 4200 CONTINUE
225 CALL GRNDM(RNDM,2)
226 RPHI1 = 2.D+00 * RNDM (1) - 1.D+00
227 RPHI2 = 2.D+00 * RNDM (2) - 1.D+00
228 RPHI12 = RPHI1 * RPHI1
229 RPHI22 = RPHI2 * RPHI2
230 RSQ = RPHI12 + RPHI22
231 IF ( RSQ .GT. 1.D+00 ) GO TO 4200
232 SINPHI = 2.D+00 * RPHI1 * RPHI2 / RSQ
233 COSPHI = ( RPHI12 - RPHI22 ) / RSQ
234 SINT02 = 1.D+00 - CZIMPC * CZIMPC
235 IF ( SINT02 .LT. ANGLSQ ) THEN
236 UBIMPC = COSPHI
237 VBIMPC = SINPHI
238 WBIMPC = 0.D+00
239 ELSE
240 SINTH0 = SQRT ( SINT02 )
241 SINPH0 = CYIMPC / SINTH0
242 COSPH0 = CXIMPC / SINTH0
243 UBIMPC = COSPHI * COSPH0 * CZIMPC - SINPHI * SINPH0
244 VBIMPC = COSPHI * SINPH0 * CZIMPC + SINPHI * COSPH0
245 WBIMPC = - COSPHI * SINTH0
246 END IF
247 GO TO 4500
248 ENTRY BIMNXT ( LBCHCK )
249 IF ( EKECON .LT. 2.D+00 * GAMMIN ) THEN
250 LABRST = .TRUE.
251 CALL RSTNXT
252 RETURN
253 END IF
254 4300 CONTINUE
255 BNDNUC = BNDSAV
256 SIGMAP = SIGMP0
257 SIGMAN = SIGMN0
258 4400 CONTINUE
259 CALL GRNDM(RNDM,1)
260 ANMFP = - LOG ( 1.D+00 - RNDM (1) ) / DSKRED
261 IF ( SBRES * SIGMAA .GT. ANMFP ) THEN
262 GO TO 6000
263 END IF
264 4500 CONTINUE
265 5000 CONTINUE
266 SBUSED = 0.D+00
267 NTRIAL = NTRIAL + 1
268 IF ( LLMDBR ) THEN
269 ALLMAX = LLLMAX + 1.D+00
270 CALL GRNDM(RNDM,2)
271 RNDLLL = ALLMAX * MAX ( RNDM (1), RNDM (2) )
272 LLLACT = INT (RNDLLL)
273 BIMPTR = RNDLLL * ALMBAR
274 BIMPTR = ABS (BIMPTR)
275 IF ( BIMPTR .LE. BIMCLM ) THEN
276 BFCEFF = BFCLMB
277 ELSE
278 HLPHLP = BFCHLP / BIMPTR
279 BFCEFF = 1.D+00 / ( HLPHLP + SQRT ( HLPHLP * HLPHLP
280 & + 1.D+00 ) )
281 END IF
282 BIMPTR = BIMPTR / BFCEFF
283 IF ( BIMPTR .GT. RADHLP ) GO TO 5000
284 ELSE
285 CALL GRNDM(RNDM,2)
286 BIMPTR = RADSIG * MAX ( RNDM (1), RNDM (2) )
287 IF ( BIMPTR .LE. BIMCLM ) THEN
288 BFCEFF = BFCLMB
289 ELSE
290 HLPHLP = BFCHLP / BIMPTR
291 BFCEFF = 1.D+00 / ( HLPHLP + SQRT ( HLPHLP * HLPHLP
292 & + 1.D+00 ) )
293 END IF
294 BIMPTR = BIMPTR / BFCEFF
295 END IF
296 BIMMEM = BIMPTR
297 IF ( BIMPTR .GT. RADTOT - RADPRO ) THEN
298 BIMPCT = 0.5D+00 * ( RADTOT + BIMPTR - RADPRO )
299 IF ( BIMPTR .GE. RADTOT ) THEN
300 X1 = BIMPTR - RADTOT
900dc090 301 DUMMY = 2.D+00 * X1 / ( RADPRO + X1 )
9e58b465 302 IF (ABS(DUMMY).GT.1.D0) THEN
900dc090 303 PRINT *,"Warning in GEANT321/peanut/bimsel.F "
304 PRINT *,"Illegal ACOS argument ",DUMMY
9e58b465 305 DUMMY = SIGN(1.D0, DUMMY )
900dc090 306 ENDIF
9c4c8b17 307 ANGRED = ACOS ( DUMMY ) / PI
fe4da5cc 308 X1 = X1 / ( R0PROT * RPRONU * RADCO2 )
309 DSKRED = ( 0.5D+00 * X1 * X1 + X1 + 1.D+00 ) * EXP (-X1)
310 & * ANGRED
9e58b465 311 IF (DSKRED.EQ.0.D0) DSKRED=1.D+00
fe4da5cc 312 ELSE
313 X1 = RADPRO + BIMPTR - RADTOT
900dc090 314 DUMMY = 2.D+00 * X1 / ( RADPRO + X1 )
9e58b465 315 IF (ABS(DUMMY).GT.1.D0) THEN
900dc090 316 PRINT *,"Warning in GEANT321/peanut/bimsel.F "
317 PRINT *,"Illegal ACOS argument ",DUMMY
9e58b465 318 DUMMY = SIGN(1.D0, DUMMY )
900dc090 319 ENDIF
9c4c8b17 320 ANGRED = ACOS ( DUMMY ) / PI
fe4da5cc 321 X1 = X1 / ( R0PROT * RPRONU * RADCO2 )
322 DSKRED = 1.D+00 - ( 0.5D+00 * X1 * X1 + X1 + 1.D+00 )
323 & * EXP (-X1) * ANGRED
324 END IF
325 ELSE
326 DSKRED = 1.D+00
327 BIMPCT = BIMPTR
328 END IF
329 IF ( .NOT. LBCHCK ) THEN
330 RHOSAV = RHOBIM
331 RHOBIM = FRHONC ( BIMPCT )
332 IF ( RHOBIM .EQ. RHOSAV ) GO TO 5500
333 PFRBIM = FPFRNC ( RHOBIM, ITNCMX )
334 EKFBIM = FEKFNC ( PFRBIM, ITNCMX )
335 RHOHLP = FRHONC ( BIMPTR )
336 PFRHLP = FPFRNC ( RHOHLP, IPWELL )
337 PFRHLP = 0.5D+00 * PFRHLP * PFRHLP / AMNUSQ (IPWELL)
338 IF ( BIMPTR .GT. RADTOT ) BNDNUC = BNDNUC * ( 1.D+00
339 & - ( BIMPTR - RADTOT ) / ( RADHLP - RADTOT ) )
340 VPRBIM = WLLRED * ( AMNUCL (IPWELL) * PFRHLP
341 & * ( 1.D+00 - 0.5D+00 * PFRHLP ) + BNDNUC )
342 LFERMI = .TRUE.
343 EKESIG = EKECON
344 PPRSIG = PNUCCO
345 CALL SIGFER ( KPTOIP (KPROJ), EKESIG, PPRSIG, LFERMI )
346 PRCOLP = ZUSFL / AUSFL * SIGMAP
347 PRCOLN = ( AUSFL - ZUSFL ) / AUSFL * SIGMAN
348 SIGMAA = PRCOLP + PRCOLN
349 PRCOLP = PRCOLP / SIGMAA
350 PRCOLN = 1.D+00 - PRCOLP
351 5500 CONTINUE
352 ELSE
353 RHOBIM = FRHONC ( BIMPCT )
354 END IF
355 XBIMPC = UBIMPC * BIMPCT
356 YBIMPC = VBIMPC * BIMPCT
357 ZBIMPC = WBIMPC * BIMPCT
358 CALL GRNDM(RNDM,1)
359 ANMFP = - LOG ( 1.D+00 - RNDM (1) ) / DSKRED
360 IF ( BIMPCT .GT. RAD1O2 ) THEN
361 SBTTSQ = 4.D+00 * ( RADTOT**2 - BIMPCT**2 ) * RHOBIM**2
362 IF ( SBTTSQ .LE. ( ANMFP / SIGMAA )**2 ) GO TO 5000
363 END IF
364 CALL SBCOMP ( SBHAL0, SBSKI0, SBCEN0, SBCEN1, SBSKI1, SBHAL1 )
365 SBTOT = SBHAL0 + SBSKI0 + SBCEN0 + SBCEN1 + SBSKI1 + SBHAL1
366 SBTOT = RHRUSF * SBTOT
367 IF ( SBTOT * SIGMAA .LE. ANMFP ) GO TO 5000
368 6000 CONTINUE
369 SBUSED = SBUSED * RHRUSF + ANMFP / SIGMAA
370 SBRES = SBTOT - SBUSED
371 SBUSED = SBUSED / RHRUSF
372 LELSTC = .TRUE.
373 NTARGT = 1
374 CALL GRNDM(RNDM,1)
375 IF ( RNDM (1) .LT. PRCOLP ) THEN
376 KNUCIM = 1
377 ITFRMI = 1
378 ELSE
379 KNUCIM = 8
380 ITFRMI = 2
381 END IF
382 IPRTYP = KPROJ * 10 + KNUCIM
383 CALL RSCOMP ( SBHAL0, SBSKI0, SBCEN0, SBCEN1, SBSKI1, SBHAL1 )
384 BNDNUC = BNDSAV
385 RHOIMP = FRHONC ( ABS (RIMPCT) )
386 PFRIMP = FPFRNC ( RHOIMP, ITFRMI )
387 EKFIMP = FEKFNC ( PFRIMP, ITFRMI )
388 RIMHLP = ABS (RIMPTR)
389 RHOIMT = FRHONC ( RIMHLP )
390 PFRPRO = FPFRNC ( RHOIMT, IPWELL )
391 EKFPRO = FEKFNC ( PFRPRO, IPWELL )
392 IF ( RIMHLP .GT. RADTOT ) BNDNUC = BNDNUC * ( 1.D+00 - ( RIMHLP
393 & - RADTOT ) / ( RADHLP - RADTOT ))
394 VPRWLL = WLLRED * ( EKFPRO + BNDNUC )
395 EKEWLL = EKECON + VPRWLL
396 EPSWLL = EKEWLL + AM (KPROJ)
397 IF ( .NOT. LBCHCK ) THEN
398 PPRWLL = SQRT ( EKEWLL * ( EKEWLL + 2.D+00 * AM (KPROJ) ) )
399 CALL PHDWLL ( UBIMPC, VBIMPC, WBIMPC )
400 PNFRMI = PFNCLV ( ITFRMI, .TRUE. )
401 IF ( PNFRMI .LT. -100.D+00 ) GO TO 4400
402 CALL RACO ( PXFERM, PYFERM, PZFERM )
403 PXFERM = PXFERM * PNFRMI
404 PYFERM = PYFERM * PNFRMI
405 PZFERM = PZFERM * PNFRMI
406 ERES = EKEWLL + AM (KPROJ) + AM (KNUCIM) + EKFERM
407 PXRES = PXPROJ + PXFERM
408 PYRES = PYPROJ + PYFERM
409 PZRES = PZPROJ + PZFERM
410 PTRES2 = PXRES**2 + PYRES**2 + PZRES**2
411 UMO2 = ERES**2 - PTRES2
412 EKESIG = 0.5D+00 * ( UMO2 - AM (KPROJ)**2 - AM (KNUCIM)**2 )
413 & / AM (KNUCIM) - AM (KPROJ)
414 EKFIMP = MAX ( EKFERM, EKFIMP )
415 ELSE
416 EKESIG = EPSWLL - AM (KPROJ)
417 END IF
418 PPRSIG = SQRT ( EKESIG * ( EKESIG + 2.D+00 * AM (KPROJ) ) )
419 SIGMN0 = SIGMAN
420 SIGMP0 = SIGMAP
421 LFERMI = .FALSE.
422 CALL SIGFER ( KPTOIP (KPROJ), EKESIG, PPRSIG, LFERMI )
423 IF ( KNUCIM .EQ. 1 ) THEN
424 SIGMAR = SIGMAP / SIGMP0
425 ELSE
426 SIGMAR = SIGMAN / SIGMN0
427 END IF
428 SIGMAR = MIN ( SIGMAR, ONEONE )
429 CALL GRNDM(RNDM,1)
430 RNDREJ = RNDM(1)
431 IF ( RNDREJ .GE. SIGMAR ) GO TO 4300
432 IF ( LBCHCK ) THEN
433 ZITA = 0.5D+00 * ( EKFIMP + EKFPRO ) / EKEWLL
434 IF ( ZITA .LE. 0.5D+00 ) THEN
435 PZITA = 1.D+00 - 1.4D+00 * ZITA
436 ELSE
437 PZITA = 1.D+00 - 1.4D+00 * ZITA + 0.4D+00 * ZITA
438 & * ( 2.D+00 - 1.D+00 / ZITA )**2.5D+00
439 END IF
440 RNDREJ = RNDREJ / SIGMAR
441 IF ( RNDREJ .GE. PZITA ) GO TO 4300
442 ELSE
443 PZITA = 1.D+00
444 END IF
445 OPACTY = 1.D+00 / NTRIAL
446 RETURN
447*=== End of subroutine Bimsel =========================================*
448 END