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