This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / peanut / nucnuc.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:22:01  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.46  by  S.Giani
11 *-- Author :
12 *$ CREATE NUCNUC.FOR
13 *COPY NUCNUC
14 *
15 *=== nucnuc ===========================================================*
16 *
17       SUBROUTINE NUCNUC ( IKPMX , KRFLIN, WEE   , ERECMN, LBIMPC,
18      &                    LBCHCK, ICYCL , NHOLE , NPROT , NNEUT ,
19      &                    LEXIT , LNWINT )
20  
21 #include "geant321/dblprc.inc"
22 #include "geant321/dimpar.inc"
23 #include "geant321/iounit.inc"
24 *
25 *----------------------------------------------------------------------*
26 *----------------------------------------------------------------------*
27 *
28 #include "geant321/balanc.inc"
29 #include "geant321/finuc.inc"
30 #include "geant321/nucdat.inc"
31 #include "geant321/nucgeo.inc"
32 #include "geant321/parevt.inc"
33 #include "geant321/parnuc.inc"
34 #include "geant321/part.inc"
35 #include "geant321/resnuc.inc"
36 *
37       REAL RNDM(2)
38       LOGICAL LBCHCK, LBIMPC, LTROUB, LEXIT, LNWINT
39 *
40       NPNCLD = NPNUC
41  1000 CONTINUE
42       IF ( LELSTC ) THEN
43          LELSTC = .FALSE.
44          LNWINT = .FALSE.
45       ELSE
46          LEXIT  = .FALSE.
47          LNWINT = .TRUE.
48          RETURN
49       END IF
50       NHOLE = NHOLE + 1
51       ICYCL = ICYCL + 1
52       IF ( NUSCIN .LE. 0 ) THEN
53          PXHLP  = PXFERM + PXPROJ - CXIMPC * PPRWLL
54          PYHLP  = PYFERM + PYPROJ - CYIMPC * PPRWLL
55          PZHLP  = PZFERM + PZPROJ - CZIMPC * PPRWLL
56          ERECMN = 0.5D+00 * ( PXHLP**2 + PYHLP**2 + PZHLP**2 ) / AMNTAR
57          ERECMN = ERECMN * ( 1.D+00 - 0.25D+00 * ERECMN )
58          ERECMN = MAX ( ERECMN, EKFERM * AM (KNUCIM) / AMNTAR )
59       ELSE
60          ERECMN = EKFERM * AM (KNUCIM) / AMNTAR
61       END IF
62       ERECMN = 0.D+00
63       POTINC = EKEWLL - EKECON + EKFERM
64       IF ( NUSCIN .EQ. 0 .AND. KPRIN .NE. KNUCIM ) THEN
65          ZAFT   = ZNOW + ICH (KPRIN) - ICH (KNUCIM)
66          AMMAFT = ANOW * AMUAMU + 0.001D+00 * FKENER ( ANOW, ZAFT )
67          AMNAFT = AMMAFT - ZAFT * AMELEC + ELBNDE (NINT(ZAFT))
68          IF ( KPRIN .EQ. 1 ) THEN
69             DEFRPR = DEFPRO
70          ELSE
71             DEFRPR = DEFNEU
72          END IF
73          IF ( KNUCIM .EQ. 8 ) THEN
74             DEFRNU = MAX ( DEFNEU, HLFHLF * EEXANY + ERECMN
75      &             + EKFERM - EKFIMP )
76          ELSE
77             DEFRNU = MAX ( DEFPRO, HLFHLF * EEXANY + ERECMN
78      &             + EKFERM - EKFIMP )
79          END IF
80       ELSE IF ( NUSCIN .EQ. 0 ) THEN
81          IF ( KPRIN .EQ. 1 ) THEN
82             DEFRPR = MAX ( DEFPRO, HLFHLF * EEXANY + ERECMN
83      &             + EKFERM - EKFIMP )
84             DEFRNU = DEFRPR
85          ELSE
86             DEFRPR = MAX ( DEFNEU, HLFHLF * EEXANY + ERECMN
87      &             + EKFERM - EKFIMP )
88             DEFRNU = DEFRPR
89          END IF
90       ELSE
91          ADDPRO = 0.D+00
92          ADDTAR = 0.D+00
93          IF ( KPRIN .EQ. 8 ) THEN
94             DEFRPR = DEFNEU + ADDPRO
95          ELSE
96             DEFRPR = DEFPRO + ADDPRO
97          END IF
98          IF ( KNUCIM .EQ. 8 ) THEN
99             DEFRNU = DEFNEU + ADDTAR
100          ELSE
101             DEFRNU = DEFPRO + ADDTAR
102          END IF
103       END IF
104       NPLSIN = 2
105       UMO2   = ERES**2 - PTRES2
106       UMO    = SQRT (UMO2)
107       GAMCM = ERES  / UMO
108       ETAX  = PXRES / UMO
109       ETAY  = PYRES / UMO
110       ETAZ  = PZRES / UMO
111       ETAPCM = ETAX * PXPROJ + ETAY * PYPROJ + ETAZ * PZPROJ
112       ECMSPR = GAMCM * ( EKEWLL + AM (KPRIN) ) - ETAPCM
113       PHELP  = EKEWLL + AM (KPRIN) - ETAPCM / ( GAMCM + 1.D+00 )
114       PCMSX  = PXPROJ - ETAX * PHELP
115       PCMSY  = PYPROJ - ETAY * PHELP
116       PCMSZ  = PZPROJ - ETAZ * PHELP
117       ETAPCM = ETAX * PXFERM + ETAY * PYFERM + ETAZ * PZFERM
118       ECMSNU = GAMCM * ( EKFERM + AM (KNUCIM) ) - ETAPCM
119       PCMS2  = PCMSX**2 + PCMSY**2 + PCMSZ**2
120       PCMS   = SQRT ( PCMS2 )
121       IF ( KPRIN .EQ. KNUCIM ) THEN
122          EKESAM = 0.5D+00 * ( UMO2 - AM (KPRIN)**2 - AM (KNUCIM)**2 )
123      &          / AM (KNUCIM) - AM (KPRIN)
124          CALL SAMCST ( 1, EKESAM, COSTHE )
125       ELSE IF ( KPRIN .EQ. 8 ) THEN
126          EKESAM = 0.5D+00 * ( UMO2 - AM (KPRIN)**2 - AM (KNUCIM)**2 )
127      &          / AM (KNUCIM) - AM (KPRIN)
128          CALL SAMCST ( 8, EKESAM, COSTHE )
129       ELSE
130          EKESAM = 0.5D+00 * ( UMO2 - AM (KNUCIM)**2 - AM (KPRIN)**2 )
131      &          / AM (KPRIN) - AM (KNUCIM)
132          CALL SAMCST ( 8, EKESAM, COSTHE )
133          COSTHE = - COSTHE
134       END IF
135       SINTHE = SQRT ( ( 1.D+00 - COSTHE ) * ( 1.D+00 + COSTHE ) )
136  2000 CONTINUE
137          CALL GRNDM(RNDM,2)
138          RPHI1 = 2.D+00 * RNDM (1) - 1.D+00
139          RPHI2 = 2.D+00 * RNDM (2) - 1.D+00
140          RPHI12 = RPHI1 * RPHI1
141          RPHI22 = RPHI2 * RPHI2
142          RSQ = RPHI12 + RPHI22
143       IF ( RSQ .GT. 1.D+00 ) GO TO 2000
144       SINPHI = 2.D+00 * RPHI1 * RPHI2 / RSQ
145       COSPHI = ( RPHI12 - RPHI22 ) / RSQ
146       CZAXCM = PCMSZ / PCMS
147       SINT02 = ( 1.D+00 - CZAXCM ) * ( 1.D+00 + CZAXCM )
148       IF ( SINT02 .LT. ANGLSQ ) THEN
149          CXCMS = COSPHI * SINTHE
150          CYCMS = SINPHI * SINTHE
151          CZCMS = CZAXCM * COSTHE
152       ELSE
153          SINTH0 = SQRT ( SINT02 )
154          UPRIME = SINTHE * COSPHI
155          VPRIME = SINTHE * SINPHI
156          CXAXCM = PCMSX / PCMS
157          CYAXCM = PCMSY / PCMS
158          COSPH0 = CXAXCM / SINTH0
159          SINPH0 = CYAXCM / SINTH0
160          CXCMS = UPRIME * COSPH0 * CZAXCM - VPRIME * SINPH0
161      &         + COSTHE * CXAXCM
162          CYCMS = UPRIME * SINPH0 * CZAXCM + VPRIME * COSPH0
163      &         + COSTHE * CYAXCM
164          CZCMS = - UPRIME * SINTH0 + COSTHE * CZAXCM
165       END IF
166       PCMSX  = PCMS * CXCMS
167       PCMSY  = PCMS * CYCMS
168       PCMSZ  = PCMS * CZCMS
169       NPNUC = NPNUC + 1
170       KPNUCL (NPNUC) = KPRIN
171       KRFNUC (NPNUC) = KRFLIN + 1
172       ETAPCM = ETAX * PCMSX + ETAY * PCMSY + ETAZ * PCMSZ
173       PHELP  = ECMSPR + ETAPCM / ( GAMCM + 1.D+00 )
174       ENNUC  (NPNUC) = GAMCM * ECMSPR + ETAPCM
175       IF ( ENNUC (NPNUC) - AM (KPRIN) .LE. EKFPRO + DEFRPR ) THEN
176          NPNUC  = NPNUC - 1
177          LBCHCK = .FALSE.
178          IF ( LBIMPC ) THEN
179             CALL BIMNXT ( LBCHCK )
180             RHOMEM = 0.5D+00 * ( RHOIMP + RHOIMT )
181             EKFMEM = 0.5D+00 * ( EKFIMP + EKFPRO )
182          ELSE
183             CALL NWINXT ( LBCHCK )
184             IF ( BIMPCT .GT. RADTOT ) THEN
185                NHOLE = NHOLE - 1
186                ICYCL = ICYCL - 1
187                CALL PHDSET ( IKPMX )
188                IBRES = IBRES - IBAR (KPRIN)
189                ICRES = ICRES - ICH  (KPRIN)
190                BBRES = IBRES
191                ZZRES = ICRES
192                AMMRES = BBRES * AMUAMU + 0.001D+00 * FKENER
193      &                ( BBRES, ZZRES)
194                AMNRES = AMMRES - ZZRES * AMELEC + ELBNDE ( ICRES )
195                LTROUB = .FALSE.
196                CALL UMOFIN ( IKPMX, BBRES, ZZRES, LTROUB )
197                IF ( LTROUB ) THEN
198                   KPNUCL (IKPMX) = 0
199                   UMO2  = ERES**2 - PTRES2
200                   UMO = SQRT (UMO2)
201                   WRITE ( LUNOUT,* )' 0_P:UMO,AMNRES',UMO,AMNRES
202                   IF ( KPRIN .EQ. 1 ) THEN
203                      NPROT = NPROT + 1
204                   ELSE
205                      NNEUT = NNEUT + 1
206                   END IF
207                   LEXIT = .TRUE.
208                   RETURN
209                END IF
210                EKNNUC = ENNUC (IKPMX) - AM (KPRIN)
211                NP = NP + 1
212                TKI   (NP) = ENNUC  (IKPMX) - AM (KPRIN)
213                KPART (NP) = KPRIN
214                PLR   (NP) = PNUCL  (IKPMX)
215                CXR   (NP) = PXNUCL (IKPMX) / PLR (NP)
216                CYR   (NP) = PYNUCL (IKPMX) / PLR (NP)
217                CZR   (NP) = PZNUCL (IKPMX) / PLR (NP)
218                WEI   (NP) = WEE
219                KPNUCL (IKPMX) = 0
220                IGREYP = IGREYP + ICH (KPRIN)
221                IGREYN = IGREYN + 1 - ICH (KPRIN)
222                PXINTR = PXINTR + PXNUCL (IKPMX)
223                PYINTR = PYINTR + PYNUCL (IKPMX)
224                PZINTR = PZINTR + PZNUCL (IKPMX)
225                EINTR  = EINTR  + ENNUC  (IKPMX)
226                IBINTR = IBINTR + IBAR   (KPART(NP))
227                ICINTR = ICINTR + ICH    (KPART(NP))
228                LEXIT  = .TRUE.
229                RETURN
230             END IF
231             XSTNUC (IKPMX) = XIMPTR
232             YSTNUC (IKPMX) = YIMPTR
233             ZSTNUC (IKPMX) = ZIMPTR
234             RSTNUC (IKPMX) = ABS (RIMPTR)
235          END IF
236          NHOLE = NHOLE - 1
237          ICYCL = ICYCL - 1
238          GO TO 1000
239       END IF
240       EKFNUC (NPNUC) = EKFPRO
241       PXNUCL (NPNUC) = PCMSX + ETAX * PHELP
242       PYNUCL (NPNUC) = PCMSY + ETAY * PHELP
243       PZNUCL (NPNUC) = PCMSZ + ETAZ * PHELP
244       PNUCL  (NPNUC) = SQRT ( PXNUCL (NPNUC)**2 + PYNUCL (NPNUC)**2
245      &                      + PZNUCL (NPNUC)**2 )
246       XSTNUC (NPNUC) = XIMPTR
247       YSTNUC (NPNUC) = YIMPTR
248       ZSTNUC (NPNUC) = ZIMPTR
249       RSTNUC (NPNUC) = ABS (RIMPTR)
250       RHNUCL (NPNUC) = RHOIMT
251       NPNUC = NPNUC + 1
252       KPNUCL (NPNUC) = KNUCIM
253       KRFNUC (NPNUC) = KRFLIN + 1
254       ETAPCM = - ETAPCM
255       PHELP  = ECMSNU + ETAPCM / ( GAMCM + 1.D+00 )
256       ENNUC  (NPNUC) = GAMCM * ECMSNU + ETAPCM
257       IF ( ENNUC (NPNUC) - AM (KNUCIM) .LE. EKFIMP + DEFRNU ) THEN
258          NPNUC  = NPNUC - NPLSIN
259          LBCHCK = .FALSE.
260          IF ( LBIMPC ) THEN
261             CALL BIMNXT ( LBCHCK )
262             RHOMEM = 0.5D+00 * ( RHOIMP + RHOIMT )
263             EKFMEM = 0.5D+00 * ( EKFIMP + EKFPRO )
264          ELSE
265             CALL NWINXT ( LBCHCK )
266             IF ( BIMPCT .GT. RADTOT ) THEN
267                NHOLE = NHOLE - 1
268                ICYCL = ICYCL - 1
269                CALL PHDSET ( IKPMX )
270                IBRES = IBRES - IBAR (KPRIN)
271                ICRES = ICRES - ICH  (KPRIN)
272                BBRES = IBRES
273                ZZRES = ICRES
274                AMMRES = BBRES * AMUAMU + 0.001D+00 * FKENER
275      &                ( BBRES, ZZRES)
276                AMNRES = AMMRES - ZZRES * AMELEC + ELBNDE ( ICRES )
277                LTROUB = .FALSE.
278                CALL UMOFIN ( IKPMX, BBRES, ZZRES, LTROUB )
279                IF ( LTROUB ) THEN
280                   KPNUCL (IKPMX) = 0
281                   UMO2  = ERES**2 - PTRES2
282                   UMO = SQRT (UMO2)
283                   WRITE ( LUNOUT,* )' 0_T:UMO,AMNRES',UMO,AMNRES
284                   IF ( KPRIN .EQ. 1 ) THEN
285                      NPROT = NPROT + 1
286                   ELSE
287                      NNEUT = NNEUT + 1
288                   END IF
289                   LEXIT = .TRUE.
290                   RETURN
291                END IF
292                EKNNUC = ENNUC (IKPMX) - AM (KPRIN)
293                NP = NP + 1
294                TKI   (NP) = ENNUC  (IKPMX) - AM (KPRIN)
295                KPART (NP) = KPRIN
296                PLR   (NP) = PNUCL  (IKPMX)
297                CXR   (NP) = PXNUCL (IKPMX) / PLR (NP)
298                CYR   (NP) = PYNUCL (IKPMX) / PLR (NP)
299                CZR   (NP) = PZNUCL (IKPMX) / PLR (NP)
300                WEI   (NP) = WEE
301                KPNUCL (IKPMX) = 0
302                IGREYP = IGREYP + ICH (KPRIN)
303                IGREYN = IGREYN + 1 - ICH (KPRIN)
304                PXINTR = PXINTR + PXNUCL (IKPMX)
305                PYINTR = PYINTR + PYNUCL (IKPMX)
306                PZINTR = PZINTR + PZNUCL (IKPMX)
307                EINTR  = EINTR  + ENNUC  (IKPMX)
308                IBINTR = IBINTR + IBAR   (KPART(NP))
309                ICINTR = ICINTR + ICH    (KPART(NP))
310                LEXIT = .TRUE.
311                RETURN
312             END IF
313             XSTNUC (IKPMX) = XIMPTR
314             YSTNUC (IKPMX) = YIMPTR
315             ZSTNUC (IKPMX) = ZIMPTR
316             RSTNUC (IKPMX) = ABS (RIMPTR)
317          END IF
318          NHOLE = NHOLE - 1
319          ICYCL = ICYCL - 1
320          GO TO 1000
321       END IF
322       EKFNUC (NPNUC) = EKFIMP
323       PXNUCL (NPNUC) = -PCMSX + ETAX * PHELP
324       PYNUCL (NPNUC) = -PCMSY + ETAY * PHELP
325       PZNUCL (NPNUC) = -PCMSZ + ETAZ * PHELP
326       PNUCL  (NPNUC) = SQRT ( PXNUCL (NPNUC)**2 + PYNUCL (NPNUC)**2
327      &                      + PZNUCL (NPNUC)**2 )
328       XSTNUC (NPNUC) = XIMPCT
329       YSTNUC (NPNUC) = YIMPCT
330       ZSTNUC (NPNUC) = ZIMPCT
331       RSTNUC (NPNUC) = ABS (RIMPCT)
332       RHNUCL (NPNUC) = RHOIMP
333       POTOUT = ENNUC (NPNUC) - AM (KPRIN) + ENNUC (NPNUC-1) - AM(KNUCIM)
334      &       - EKECON
335       LBIMPC = .FALSE.
336       LEXIT  = .FALSE.
337       NUSCIN = NUSCIN + 1
338       ISCTYP (NUSCIN) = KPRIN * 10 + KNUCIM
339       NHLEXP = NHLEXP + 1
340       HOLEXP (NHLEXP) = EKFIMP - EKFERM
341       RHOEXP = RHOEXP + 0.5D+00 * ( RHOIMP + RHOIMT )
342       EKFEXP = EKFEXP + 0.5D+00 * ( EKFIMP + EKFPRO )
343       CALL NCLVFX
344       DO 3000 KP = NPNCLD+1, NPNUC
345          KPNUC = KPNUCL (KP)
346          IF ( AM (KPNUC) .LE. 0.D+00 ) THEN
347             TAUTAU = RZNUCL / PNUCL (KP)
348          ELSE
349             TAUEFF = 0.5D+00 * TAUFOR * AM (13) / AM (KPNUC)
350             CALL GRNDM(RNDM,1)
351             TAUTAU = - TAUEFF / AM (KPNUC) * LOG ( 1.D+00 - RNDM
352      &             (1) )
353             TAUTAU = MAX ( TAUTAU, RZNUCL / PNUCL (KP) )
354          END IF
355          XSTNUC (KP) = XSTNUC (KP) + PXNUCL (KP) * TAUTAU
356          YSTNUC (KP) = YSTNUC (KP) + PYNUCL (KP) * TAUTAU
357          ZSTNUC (KP) = ZSTNUC (KP) + PZNUCL (KP) * TAUTAU
358          RSTNUC (KP) = SQRT ( XSTNUC (KP)**2 + YSTNUC (KP)**2
359      &               + ZSTNUC (KP)**2 )
360  3000 CONTINUE
361       RETURN
362 *=== End of subroutine Nucnuc =========================================*
363       END