Use double precision constants
[u/mrichter/AliRoot.git] / GEANT321 / peanut / bimsel.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.5  2002/06/25 08:17:33  hristov
6 * Additional protection
7 *
8 * Revision 1.4  2002/06/22 10:50:18  hristov
9 * Better protection
10 *
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 *
14 * Revision 1.2  2002/05/13 12:40:58  hristov
15 * Dummy subroutines to avoid files with no code in
16 *
17 * Revision 1.1.1.1  1999/05/18 15:55:22  fca
18 * AliRoot sources
19 *
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
65 *      SAVE
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
301                DUMMY = 2.D+00 * X1 / ( RADPRO + X1 )
302                IF (ABS(DUMMY).GT.1.D0) THEN
303                   PRINT *,"Warning in GEANT321/peanut/bimsel.F "
304                   PRINT *,"Illegal ACOS argument ",DUMMY
305                   DUMMY = SIGN(1.D0, DUMMY )
306                ENDIF
307                ANGRED = ACOS ( DUMMY ) / PI
308                X1 = X1 / ( R0PROT * RPRONU * RADCO2 )
309                DSKRED = ( 0.5D+00 * X1 * X1 + X1 + 1.D+00 ) * EXP (-X1)
310      &                * ANGRED
311                IF (DSKRED.EQ.0.D0) DSKRED=1.D+00
312             ELSE
313                X1 = RADPRO + BIMPTR - RADTOT
314                DUMMY = 2.D+00 * X1 / ( RADPRO + X1 )
315                IF (ABS(DUMMY).GT.1.D0) THEN
316                   PRINT *,"Warning in GEANT321/peanut/bimsel.F "
317                   PRINT *,"Illegal ACOS argument ",DUMMY
318                   DUMMY = SIGN(1.D0, DUMMY )
319                ENDIF
320                ANGRED = ACOS ( DUMMY ) / PI
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