Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / GEANT321 / fluka / fkdres.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:20:05  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.45  by  S.Giani
11 *-- Author :
12 *=== dres =============================================================*
13 *                                                                      *
14       SUBROUTINE FKDRES ( M2, M3, T1, U, EREC, LOPPAR, JFISS )
15  
16 #include "geant321/dblprc.inc"
17 #include "geant321/dimpar.inc"
18 #include "geant321/iounit.inc"
19 *
20 *----------------------------------------------------------------------*
21 *                                                                      *
22 *  New version of DRES created  by A.Ferrari & P.Sala, INFN - Milan    *
23 *                                                                      *
24 *  Last change  on  10-apr-93  by  Alfredo Ferrari, INFN - Milan       *
25 *                                                                      *
26 *  Dres93: Dres91 plus the RAL fission model taken from LAHET thanks   *
27 *          to R.E.Prael                                                *
28 *  Dres91: new version from A. Ferrari and P. Sala, INFN - Milan       *
29 *          This routine has been adapted from the original one of the  *
30 *          Evap-5 module (KFA - Julich). Main modifications concern    *
31 *          with kinematics which is now fully relativistic and with    *
32 *          the treatment of few nucleons nuclei, which are now frag-   *
33 *          mented, even though in a very rough manner. Changes have    *
34 *          been made also to other routines of the Evap-5 package      *
35 *                                                                      *
36 *----------------------------------------------------------------------*
37 *
38 *----------------------------------------------------------------------*
39 *                                                                      *
40 *  Input variables:                                                    *
41 *     M2 = Mass number of the residual nucleus                         *
42 *     M3 = Atomic number of the residual nucleus                       *
43 *     T1 = Excitation energy of the residual nucleus before evaporation*
44 *     U  = Excitation energy of the residual nucleus after evaporation *
45 *     Erec = Recoil kinetic energy of the residual nucleus             *
46 *            The recoil direction is given by Coslbr (i)               *
47 *                                                                      *
48 *  Significant variables:                                              *
49 *     JA = Present mass number of the residual nucleus                 *
50 *     JZ = Present atomic number of the residual nucleus               *
51 *     Smom1 = Energy accumulators for the six types of evaporated      *
52 *             particles                                                *
53 *                                                                      *
54 *    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     *
55 *    !!!! Please note that the following variables concerning !!!!     *
56 *    !!!! with the present residual nucleus must be set before!!!!     *
57 *    !!!! entering DRES91: Ammres, Ptres                      !!!!     *
58 *    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     *
59 *                                                                      *
60 *----------------------------------------------------------------------*
61 *
62 C---------------------------------------------------------------------
63 C SUBNAME = DRES --- EVAPORATION
64 C           EVAPORATION DATA SHOUD BE READ ON INPUT STAGE
65 C---------------------------------------------------------------------
66 C*****A LA EVAP III(TWA,8-68)
67 #include "geant321/eva0.inc"
68 #include "geant321/eva1.inc"
69 #include "geant321/forcn.inc"
70 #include "geant321/higfis.inc"
71 #include "geant321/inpflg.inc"
72 #include "geant321/labcos.inc"
73 #include "geant321/nucdat.inc"
74 #include "geant321/resnuc.inc"
75       DIMENSION ZMASS (6), Z2MASS(6), C(3), Q(0:6), FLKCOU(6), CCOUL(6),
76      &          THRESH(6), SMALLA(6), R(6), S  (6), SOS   (6), STRUN(6),
77      &          EYE1  (6), EYE0  (6), SMOM1    (6), BNMASS(6)
78       DIMENSION CORRRR(6)
79       REAL RNDM(2)
80       LOGICAL LOPPAR, PENBAR, LFIRST
81       PARAMETER (EXPMIN=-88,EXPMAX=88)
82       SAVE ZMASS, Z2MASS, EMHN, EMNUM, UM, AMUMEV, AMEMEV, QBRBE8,
83      &     BNMASS, IEVEVP, NBE8, NRNEEP, LFIRST
84       DATA IEVEVP / 0 /
85       DATA LFIRST / .TRUE. /
86 *
87       IEVEVP = IEVEVP + 1
88 C-------------------------------------- 1.ST CALL INIT
89       IF ( LFIRST ) THEN
90          LFIRST = .FALSE.
91          FKEY   = ZERZER
92          NBE8   = 0
93          NRNEEP = 0
94          EXMASS(1) = 1.D+03 * ( AMNEUT - AMUAMU )
95          EXMASS(2) = FKENER ( ONEONE, ONEONE )
96          EXMASS(3) = FKENER ( TWOTWO, ONEONE )
97          EXMASS(4) = FKENER ( THRTHR, ONEONE )
98          EXMASS(5) = FKENER ( THRTHR, TWOTWO )
99          EXMASS(6) = FKENER ( FOUFOU, TWOTWO )
100          ZMASS(1) = 1.D+03 * AMUAMU + EXMASS (1)
101          ZMASS(2) = 1.D+03 * AMUAMU + EXMASS (2)
102          ZMASS(3) = 2.D+03 * AMUAMU + EXMASS (3)
103          ZMASS(4) = 3.D+03 * AMUAMU + EXMASS (4)
104          ZMASS(5) = 3.D+03 * AMUAMU + EXMASS (5)
105          ZMASS(6) = 4.D+03 * AMUAMU + EXMASS (6)
106          BNMASS (1) = 0.D+00
107          BNMASS (2) = 0.D+00
108          BNMASS (3) = ZMASS (1) + ZMASS (2) - ZMASS (3)
109          BNMASS (4) = TWOTWO * ZMASS (1) + ZMASS (2) - ZMASS (4)
110          BNMASS (5) = ZMASS (1) + TWOTWO * ZMASS (2) - ZMASS (5)
111          BNMASS (6) = TWOTWO * ( ZMASS (1) + ZMASS (2) ) - ZMASS (6)
112          DO 1234 KK = 1,6
113             Z2MASS (KK) = ZMASS (KK) * ZMASS (KK)
114 1234     CONTINUE
115          AMUMEV = 1.D+03 * AMUAMU
116          AMEMEV = 1.D+03 * AMELEC
117          QBRBE8 = FKENER ( EIGEIG, FOUFOU ) - TWOTWO * EXMASS (6)
118          EMN = 1.D+03 * AMNEUT
119          EMH = ZMASS (2)
120          TMP16 = 16.D+00
121          UM  = AMUMEV + FKENER ( TMP16, EIGEIG ) / 16.D+00
122          EMHN  = EMH - EMN
123          EMNUM = EMN - UM
124       END IF
125 *  |  End of initialization:
126 *  +-------------------------------------------------------------------*
127 C     --------------------------------- START OF PROCESS
128 *  +-------------------------------------------------------------------*
129 *  |  Initialize Npart and Smom if nothing has been already evaporated
130 *  |  for this event
131       IF ( JFISS .LE. 0 ) THEN
132          DO 775 I=1,6
133             NPART(I) = 0
134             SMOM1(I) = ZERZER
135   775    CONTINUE
136       END IF
137 *  |
138 *  +-------------------------------------------------------------------*
139       JA = M2
140       JZ = M3
141       U  = T1
142       RNMASS = 1.D+03 * AMMRES + U
143 * P2res and  Ptres are the squared momentum and the momentum of the
144 * residual nucleus (now in relativistic kinematics), Umo the
145 * invariant mass of the system!
146       UMO  = RNMASS
147       UMO2 = UMO * UMO
148       ELBTOT  = RNMASS + EREC
149       GAMCM   = ELBTOT / RNMASS
150       ETACM   = 1.D+03 * PTRES / RNMASS
151       HEVSUM  = ZERZER
152  1000 CONTINUE
153       LOPPAR = .FALSE.
154 *  +-------------------------------------------------------------------*
155 *  |                Check for starting data inconsistencies
156       IF (JA-JZ .LT. 0) THEN
157          WRITE(LUNOUT,6401)
158          WRITE(LUNERR,6401)
159  6401    FORMAT('1 Dres: cascade residual nucleus has mass no. less',
160      &       ' than Z!!')
161          RETURN
162 *  |
163 *  +-------------------------------------------------------------------*
164 *  |                Rough treatment for very few nucleon residual
165 *  |                nuclei. The basic ideas are:
166 *  |        a) as many as possible alpha particles are emitted
167 *  |        b) particles are emitted one per time leaving a residual
168 *  |           excitation energy proportional to number of nucleons
169 *  |           left in the residual nucleus (so we deal only with
170 *  |           two body kinematics)
171 *  |       T A K E   I N T O   A C C O U N T   T H A T   T H I S
172 *  |       T R E A T M E N T   I S   E X T R E M E L Y   R O U G H
173 *  |       T H E   T A S K   B E I N G   O N L Y   T O   S U P P L Y
174 *  |       S O M E T H I N G   T O   S H A R E   E N E R G Y   A N D
175 *  |       M O M E N T U M   A M O N G   A   F E W   F R A G M E N T S
176       ELSE IF ( JA .LE. 6 .OR. JZ .LE. 2 ) THEN
177 *  | 1000 continue moved above according to FCA suggestion
178 *1000    CONTINUE
179          JRESID = 0
180          IF ( JA .GT. 4 ) GO TO 2000
181 *  |  +----------------------------------------------------------------*
182 *  |  | First check we are not concerning with a couple of neutrons or
183 *  |  | protons
184          IF ( JA .EQ. 2 .AND. JZ .NE. 1 ) THEN
185             JEMISS = 1 + JZ / 2
186             JRESID = JEMISS
187             RNMASS = ZMASS (JRESID)
188             U = 0.D+00
189             DELTU = UMO - 2.D+00 * ZMASS (JEMISS)
190 *  |  |  +-------------------------------------------------------------*
191 *  |  |  |
192             IF ( DELTU .LE. 0.D+00 ) THEN
193                IF ( DELTU .LT. - 2.D+00 * ANGLGB * UMO ) THEN
194                   WRITE ( LUNERR, * )' *** Dres: insufficient Umo for',
195      &                               ' a nucleon couple', UMO,
196      &                                 2.D+00 * ZMASS (JEMISS)
197                END IF
198                UMO = ( UMO + DELTU ) * ( 1.D+00 + ANGLGB )
199             END IF
200 *  |  |  |
201 *  |  |  +-------------------------------------------------------------*
202             GO TO 2500
203          END IF
204 *  |  |
205 *  |  +----------------------------------------------------------------*
206 *  |  +----------------------------------------------------------------*
207 *  |  | Then check we are not concerning with one of the six
208 *  |  | standard particles
209          DO 1700 J = 6, 1, -1
210 *  |  |  +-------------------------------------------------------------*
211 *  |  |  |
212             IF ( JZ .EQ. IZ (J) .AND. JA .EQ. IA (J) ) THEN
213                HEVSUM = SMOM1(3) + SMOM1(5) + SMOM1(6) + SMOM1(4)
214                GO TO ( 1100, 1100, 1600, 1500, 1400, 1300 ), J
215 *  |  |  |  +----------------------------------------------------------*
216 *  |  |  |  |  Proton or neutron, nothing can be done
217  1100          CONTINUE
218                   RETURN
219 *  |  |  |  +----------------------------------------------------------*
220 *  |  |  |  |  Alpha:
221  1300          CONTINUE
222                   DEUDEU = MAX ( ZERZER, U + TWOTWO * BNMASS (3)
223      &                           - BNMASS (6) )
224                   PROTRI = MAX ( ZERZER, U + BNMASS (4) - BNMASS (6) )
225                   UEU3HE = MAX ( ZERZER, U + BNMASS (5) - BNMASS (6) )
226                   QNORM  = DEUDEU + PROTRI + UEU3HE
227 *  |  |  |  |  If we cannot split then return
228                   IF ( QNORM .LE. ZERZER ) RETURN
229                   CALL GRNDM(RNDM,1)
230                   V = RNDM (1)
231 *  |  |  |  |  +-------------------------------------------------------*
232 *  |  |  |  |  | Split or into two deuterons or a triton and a proton
233 *  |  |  |  |  | or a 3-He and a neutron: no account is made for
234 *  |  |  |  |  | Coulomb effects, probability is simply assumed
235 *  |  |  |  |  | proportional to reaction Qs
236                   IF ( V .LT. DEUDEU / QNORM ) THEN
237 *  |  |  |  |  | Two deuterons selected
238                      JEMISS = 3
239                      JRESID = 3
240                      RNMASS = ZMASS (3)
241                      U = ZERZER
242                      GO TO 2500
243 *  |  |  |  |  |
244 *  |  |  |  |  +-------------------------------------------------------*
245 *  |  |  |  |  | Split into a triton and a proton
246                   ELSE IF ( V .LT. ( DEUDEU + PROTRI ) / QNORM ) THEN
247                      JEMISS = 2
248                      JRESID = 4
249                      RNMASS = ZMASS (4)
250                      U = ZERZER
251                      GO TO 2500
252 *  |  |  |  |  |
253 *  |  |  |  |  +-------------------------------------------------------*
254 *  |  |  |  |  | Split into a 3-He and a neutron
255                   ELSE
256                      JEMISS = 1
257                      JRESID = 5
258                      RNMASS = ZMASS (5)
259                      U = ZERZER
260                      GO TO 2500
261                   END IF
262 *  |  |  |  |  |
263 *  |  |  |  |  +-------------------------------------------------------*
264 *  |  |  |  |
265 *  |  |  |  +----------------------------------------------------------*
266 *  |  |  |  |  3-He:
267  1400          CONTINUE
268                   DEUPRO = MAX ( ZERZER, U + BNMASS (3) - BNMASS (5) )
269                   PRPRNE = MAX ( ZERZER, U - BNMASS (5) )
270                   QNORM  = DEUPRO + PRPRNE
271 *  |  |  |  |  If we cannot split then return
272                   IF ( QNORM .LE. ZERZER ) RETURN
273                   CALL GRNDM(RNDM,1)
274                   V = RNDM (1)
275 *  |  |  |  |  +-------------------------------------------------------*
276 *  |  |  |  |  | Split or into a deuteron and a proton
277 *  |  |  |  |  | or into two protons and one neutron: no account is
278 *  |  |  |  |  | made for Coulomb effects, probability is simply assumed
279 *  |  |  |  |  | prportional to reaction Qs
280                   IF ( V .LT. DEUPRO / QNORM ) THEN
281 *  |  |  |  |  | A deuteron and a proton selected
282                      JEMISS = 2
283                      JRESID = 3
284                      RNMASS = ZMASS (3)
285                      U = ZERZER
286                      GO TO 2500
287 *  |  |  |  |  |
288 *  |  |  |  |  +-------------------------------------------------------*
289 *  |  |  |  |  | Split into 2 protons and 1 neutron: part of the exci-
290 *  |  |  |  |  | tation energy is conserved to allow the further
291 *  |  |  |  |  | splitting of the deuteron
292                   ELSE
293                      JEMISS = 2
294                      JRESID = 0
295                      FACT = ONEONE
296 *  |  |  |  |  |  +----------------------------------------------------*
297 *  |  |  |  |  |  | Loop to compute the residual excitation energy
298  1450                CONTINUE
299                         FACT = FACT * 0.6666666666666667D+00
300 *  |  |  |  |  |  | Erncm, Eepcm are the total energies of the residual
301 *  |  |  |  |  |  | nucleus and of the emitted particle in the CMS frame
302                         U      = FACT * PRPRNE + BNMASS (3)
303                         RNMASS = ZMASS (3) + U
304                         ERNCM  = HLFHLF * ( UMO2 + RNMASS**2
305      &                         - Z2MASS (JEMISS) ) / UMO
306                         EEPCM  = UMO - ERNCM
307                      IF ( EEPCM .LE. ZMASS (JEMISS) ) GO TO 1450
308 *  |  |  |  |  |  |
309 *  |  |  |  |  |  +----------------------------------------------------*
310                      GO TO 2600
311                   END IF
312 *  |  |  |  |  |
313 *  |  |  |  |  +-------------------------------------------------------*
314 *  |  |  |  |
315 *  |  |  |  +----------------------------------------------------------*
316 *  |  |  |  |  Triton:
317  1500          CONTINUE
318                   DEUNEU = MAX ( ZERZER, U + BNMASS (3) - BNMASS (4) )
319                   PRNENE = MAX ( ZERZER, U - BNMASS (4) )
320                   QNORM  = DEUNEU + PRNENE
321 *  |  |  |  |  If we cannot split then return
322                   IF ( QNORM .LE. ZERZER ) RETURN
323                   CALL GRNDM(RNDM,1)
324                   V = RNDM (1)
325 *  |  |  |  |  +-------------------------------------------------------*
326 *  |  |  |  |  | Split or into a deuteron and a neutron
327 *  |  |  |  |  | or into two protons and one neutron: no account is
328 *  |  |  |  |  | made for Coulomb effects, probability is simply assumed
329 *  |  |  |  |  | proportional to reaction Qs
330                   IF ( V .LT. DEUNEU / QNORM ) THEN
331 *  |  |  |  |  | A deuteron and a proton selected
332                      JEMISS = 1
333                      JRESID = 3
334                      RNMASS = ZMASS (3)
335                      U = ZERZER
336                      GO TO 2500
337 *  |  |  |  |  |
338 *  |  |  |  |  +-------------------------------------------------------*
339 *  |  |  |  |  | Split into 1 proton and 2 neutrons: part of the exci-
340 *  |  |  |  |  | tation energy is conserved to allow the further
341 *  |  |  |  |  | splitting of the deuteron
342                   ELSE
343                      JEMISS = 1
344                      JRESID = 0
345                      FACT = ONEONE
346 *  |  |  |  |  |  +----------------------------------------------------*
347 *  |  |  |  |  |  | Loop to compute the residual excitation energy
348  1550                CONTINUE
349                         FACT = FACT * 0.6666666666666667D+00
350 *  |  |  |  |  |  | Erncm, Eepcm are the total energies of the residual
351 *  |  |  |  |  |  | nucleus and of the emitted particle in the CMS frame
352                         U      = FACT * PRNENE + BNMASS (3)
353                         RNMASS = ZMASS (3) + U
354                         ERNCM  = HLFHLF * ( UMO2 + RNMASS**2
355      &                         - Z2MASS (JEMISS) ) / UMO
356                         EEPCM  = UMO - ERNCM
357                      IF ( EEPCM .LE. ZMASS (JEMISS) ) GO TO 1550
358 *  |  |  |  |  |  |
359 *  |  |  |  |  |  +----------------------------------------------------*
360                      GO TO 2600
361                   END IF
362 *  |  |  |  |  |
363 *  |  |  |  |  +-------------------------------------------------------*
364 *  |  |  |  |
365 *  |  |  |  +----------------------------------------------------------*
366 *  |  |  |  |  Deuteron:
367  1600          CONTINUE
368 *  |  |  |  |  +-------------------------------------------------------*
369 *  |  |  |  |  | Split into a proton and a neutron if it is possible
370                   IF ( U .GT. BNMASS (3) ) THEN
371                      JEMISS = 1
372                      JRESID = 2
373                      RNMASS = ZMASS (2)
374                      U = ZERZER
375                      GO TO 2500
376 *  |  |  |  |  |
377 *  |  |  |  |  +-------------------------------------------------------*
378 *  |  |  |  |  | Energy too low to split the deuteron, return
379                   ELSE
380                      WRITE (LUNERR,*)' **Dres: energy too low to split',
381      &                               ' a deuteron! M2,M3',M2,M3
382                      RETURN
383                   END IF
384 *  |  |  |  |  |
385 *  |  |  |  |  +-------------------------------------------------------*
386             END IF
387 *  |  |  |
388 *  |  |  +-------------------------------------------------------------*
389  1700    CONTINUE
390 *  |  |
391 *  |  +----------------------------------------------------------------*
392  2000    CONTINUE
393          A = JA
394          Z = JZ
395          Q (0)  = ZERZER
396          ENERG0 = FKENER (A,Z)
397 *  |  +----------------------------------------------------------------*
398 *  |  |   Note that Q(i) are not the reaction Qs but the remaining
399 *  |  |   energy after the reaction
400          DO 2100 K = 1, 6
401             JJA = JA - IA (K)
402             JJZ = JZ - IZ (K)
403             JJN = JJA - JJZ
404             IF ( JJN .LT. 0 .OR. JJZ .LT. 0 ) THEN
405                Q (K) = Q (K-1)
406                GO TO 2100
407             END IF
408             DDJJA = JJA
409             DDJJZ = JJZ
410             Q (K) = MAX ( U + ENERG0 - FKENER ( DDJJA, DDJJZ )
411      &            - EXMASS (K), ZERZER ) + Q (K-1)
412  2100    CONTINUE
413 *  |  |
414 *  |  +----------------------------------------------------------------*
415 *  |  +----------------------------------------------------------------*
416 *  |  |  If no emission channel is open then return
417          IF ( Q (6) .LE. ZERZER ) THEN
418             HEVSUM = SMOM1(3) + SMOM1(5) + SMOM1(6) + SMOM1(4)
419             RETURN
420          END IF
421 *  |  |
422 *  |  +----------------------------------------------------------------*
423          CALL GRNDM(RNDM,1)
424          V = RNDM (1)
425          FACT = ONEONE
426 *  |  +----------------------------------------------------------------*
427 *  |  |
428          DO 2200 J = 1, 6
429 *  |  |  +-------------------------------------------------------------*
430 *  |  |  |
431             IF ( V .LT. Q (J) / Q (6) ) THEN
432                JEMISS = J
433                JJA    = JA - IA (JEMISS)
434                JJZ    = JZ - IZ (JEMISS)
435 *  |  |  |  +----------------------------------------------------------*
436 *  |  |  |  |
437                DO 2150 JJ = 1, 6
438 *  |  |  |  |  +-------------------------------------------------------*
439 *  |  |  |  |  |
440                   IF ( JJA .EQ. IA (JJ) .AND. JJZ .EQ. IZ (JJ) ) THEN
441                      JRESID = JJ
442                      RNMASS = ZMASS (JRESID)
443                      ERNCM  = HLFHLF * ( UMO2 + Z2MASS (JRESID)
444      &                      - Z2MASS (JEMISS) ) / UMO
445                      EEPCM  = UMO - ERNCM
446                      U = ZERZER
447                      GO TO 2600
448                   END IF
449 *  |  |  |  |  |
450 *  |  |  |  |  +-------------------------------------------------------*
451  2150          CONTINUE
452 *  |  |  |  |
453 *  |  |  |  +----------------------------------------------------------*
454                AJJA   = JJA
455                ZJJZ   = JJZ
456                RNMAS0 = AJJA * AMUMEV + FKENER ( AJJA, ZJJZ )
457                GO TO 2300
458             END IF
459 *  |  |  |
460 *  |  |  +-------------------------------------------------------------*
461  2200    CONTINUE
462 *  |  |
463 *  |  +----------------------------------------------------------------*
464          WRITE ( LUNOUT,* )' **** error in Dres, few nucleon treatment',
465      &                     ' ****'
466          WRITE ( LUNERR,* )' **** error in Dres, few nucleon treatment',
467      &                     ' ****'
468          RETURN
469 *  |  +----------------------------------------------------------------*
470 *  |  | Loop to compute the residual excitation energy
471  2300    CONTINUE
472             FACT = FACT * AJJA / A
473             U = FACT * ( Q (JEMISS) - Q (JEMISS-1) )
474 *  |  | Erncm, Eepcm are the total energies of the residual
475 *  |  | nucleus and of the emitted particle in the CMS frame
476             RNMASS = RNMAS0 + U
477             ERNCM  = HLFHLF * ( UMO2 + RNMASS**2
478      &             - Z2MASS (JEMISS) ) / UMO
479             EEPCM  = UMO - ERNCM
480          IF ( EEPCM .LE. ZMASS (JEMISS) ) THEN
481             IF ( Q (JEMISS) - Q (JEMISS-1) .GE. 1.D-06 ) GO TO 2300
482 *  |  +--<--<--<--<--<--< Loop back
483 *  |  |  Actually there is no excitation energy available!
484             U = ANGLGB
485             RNMASS = ONEPLS * RNMAS0
486             ERNCM  = ONEPLS * RNMASS
487             EEPCM  = ONEPLS * ZMASS (JEMISS)
488          END IF
489 *  |  |
490 *  |  +----------------------------------------------------------------*
491          GO TO 2600
492 *  |  From here standard two bodies kinematics with Jemiss, Rnmass
493 *  |  (new excitation energy is U)
494  2500    CONTINUE
495 *  |  Erncm, Eepcm are the total energies of the residual
496 *  |  nucleus and of the emitted particle in the CMS frame
497          ERNCM = HLFHLF * ( UMO2 + RNMASS**2 - Z2MASS (JEMISS) ) / UMO
498          EEPCM = UMO - ERNCM
499  2600    CONTINUE
500 *  |  C(i) are the direction cosines of the emitted particle
501 *  |  (Jemiss) in the CMS frame, of course - C(i)
502 *  |  are the ones of the residual nucleus (Jresid if one of the
503 *  |  standard partcles, say the proton)
504          CALL RACO (C(1),C(2),C(3))
505          PCMS  = SQRT ( EEPCM**2 - Z2MASS (JEMISS) )
506 *  |  Now we perform the Lorentz transformation back to the original
507 *  |  frame (lab frame)
508 *  |  First the emitted particle:
509          ETAX = ETACM * COSLBR (1)
510          ETAY = ETACM * COSLBR (2)
511          ETAZ = ETACM * COSLBR (3)
512          PCMSX = PCMS * C (1)
513          PCMSY = PCMS * C (2)
514          PCMSZ = PCMS * C (3)
515          ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
516          EPS = GAMCM * EEPCM + ETAPCM - ZMASS (JEMISS)
517          PHELP = ETAPCM / ( GAMCM + ONEONE ) + EEPCM
518          PLBPX = PCMSX + ETAX * PHELP
519          PLBPY = PCMSY + ETAY * PHELP
520          PLBPZ = PCMSZ + ETAZ * PHELP
521          PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
522          COSLBP (1) = PLBPX / PHELP
523          COSLBP (2) = PLBPY / PHELP
524          COSLBP (3) = PLBPZ / PHELP
525 *  |  Then the residual nucleus ( for it c (i) --> - c (i) ):
526          EREC  = GAMCM * ERNCM - ETAPCM - RNMASS
527          EKRES = 1.D-03 * EREC
528          PHELP = - ETAPCM / ( GAMCM + ONEONE ) + ERNCM
529          PXRES = 1.D-03 * ( - PCMSX + ETAX * PHELP )
530          PYRES = 1.D-03 * ( - PCMSY + ETAY * PHELP )
531          PZRES = 1.D-03 * ( - PCMSZ + ETAZ * PHELP )
532          P2RES = PXRES * PXRES + PYRES * PYRES + PZRES * PZRES
533          PTRES = SQRT (P2RES)
534          COSLBR (1) = PXRES / PTRES
535          COSLBR (2) = PYRES / PTRES
536          COSLBR (3) = PZRES / PTRES
537 *  |  Score the emitted particle
538          NPART (JEMISS) = NPART (JEMISS) + 1
539          SMOM1 (JEMISS) = SMOM1 (JEMISS) + EPS
540          ITEMP=NPART(JEMISS)
541          EPART(ITEMP,JEMISS)=EPS
542          COSEVP(1,ITEMP,JEMISS)=COSLBP(1)
543          COSEVP(2,ITEMP,JEMISS)=COSLBP(2)
544          COSEVP(3,ITEMP,JEMISS)=COSLBP(3)
545 *  |  +----------------------------------------------------------------*
546 *  |  |  Check if the residual nucleus is one of the emitted particles
547          IF ( JRESID .GT. 0 ) THEN
548             J = JRESID
549             GO TO 6102
550          END IF
551 *  |  |
552 *  |  +----------------------------------------------------------------*
553          JA = JA - IA (JEMISS)
554          JZ = JZ - IZ (JEMISS)
555 * Umo is the invariant mass of the system!!
556          UMO  = RNMASS
557          UMO2 = UMO * UMO
558          ELBTOT  = RNMASS + EREC
559          GAMCM   = ELBTOT / RNMASS
560          ETACM   = 1.D+03 * PTRES / RNMASS
561          GO TO 1000
562       END IF
563 *  |
564 *  +-------------------------------------------------------------------*
565 * Come to 23 at the beginning and after the end of a "normal"
566 * evaporation cycle
567    23 CONTINUE
568       A = JA
569       Z = JZ
570       IF(JA-8)8486,8488,8486
571  8488 CONTINUE
572       IF(JZ-4)8486,1224,8486
573  8486 CONTINUE
574       DO 2 K=1,6
575          IF((A-FLA(K)).LE.(Z-FLZ(K))) THEN
576             Q(K)=99999.D0
577             GO TO 2
578          END IF
579          IF(A-TWOTWO*FLA(K))2,727,727
580   727    CONTINUE
581          IF(Z-TWOTWO*FLZ(K))2,728,728
582   728    CONTINUE
583          Q(K) = QNRG(A-FLA(K),Z-FLZ(K),A,Z) + EXMASS(K)
584 C 728    Q(K)=FKENER(A-FLA(K),Z-FLZ(K))-FKENER(A,Z)+EXMASS(K)
585     2 CONTINUE
586       FLKCOU(1)=ZERZER
587       FLKCOU(2)=DOST(1,Z-FLZ(2))
588       FLKCOU(3)=FLKCOU(2)+.06D+00
589       FLKCOU(4)=FLKCOU(2)+.12D+00
590       FLKCOU(6)=DOST(2,Z-FLZ(6))
591       FLKCOU(5)=FLKCOU(6)-.06D+00
592       CCOUL(1)=ONEONE
593       CCOU2=DOST(3,Z-FLZ(2))
594       CCOUL(2)=CCOU2+ONEONE
595       CCOUL(3)=CCOU2*1.5D0+THRTHR
596       CCOUL(4)=CCOU2+THRTHR
597       CCOUL(6)=DOST(4,Z-FLZ(6))*TWOTWO+TWOTWO
598       CCOUL(5)=TWOTWO*CCOUL(6)-ONEONE
599       SIGMA=ZERZER
600 *  Initialize the flag which checks for open particle decay with
601 *  zero excitation and pairing --> for particle unstable residual
602 *  nuclei
603       LOPPAR = .FALSE.
604       DO 33 J=1,6
605          IF(A-TWOTWO*FLA(J))5,725,725
606   725    CONTINUE
607          IF(Z-TWOTWO*FLZ(J))5,726,726
608   726    CONTINUE
609          MM=JA-IA(J)
610          ZZ=Z-FLZ(J)
611          AA=A-FLA(J)
612          IF(AA.LE.ZZ)GO TO 5
613 *  Energy threshold for the emission of the jth-particle
614          THRESH(J)=Q(J)+.88235D+00*FLKCOU(J)*FLZ(J)*
615      &   ZZ/(RMASS(MM)+RHO(J))
616          LOPPAR = LOPPAR .OR. THRESH (J) .GT. ZERZER
617          IAA = NINT(AA)
618          IZZ = NINT(ZZ)
619          NN  = IAA - IZZ
620 *  The residual nucleus excitation energy ranges from 0 up
621 *  to U - Q (J)
622          ILVMOD = IB0
623          UMXRES = U - THRESH (J)
624 *  This is the a lower case of the level density
625          SMALLA (J) = GETA ( UMXRES, IZZ, NN, ILVMOD, ISDUM, ASMMAX,
626      &                       ASMMIN )
627          CALL EEXLVL (IAA,IZZ,EEX1ST,EEX2ND,CORR)
628          EEX1ST = 1.D+03 * EEX1ST
629          EEX2ND = 1.D+03 * EEX2ND
630          CORR   = 1.D+03 * MAX ( CORR, ZERZER )
631          IF ( NN .EQ. 4 .AND. IZZ .EQ. 4 ) THEN
632             IF ( U - THRESH (J) - 6.1D+00 .GT. ZERZER ) THEN
633                CORR = SIXSIX
634             ELSE
635                TMPVAR = U-THRESH(J)-0.1D+00
636                CORR = MAX ( ZERZER, TMPVAR )
637             END IF
638          END IF
639          IF (NINT(FKEY).EQ.1) CORR=ZERZER
640          CORRRR(J)=CORR
641 * Standard calculation:
642          ARG=U-THRESH(J)-CORR
643          IF(ARG)5,4,4
644     5    CONTINUE
645          R(J)=ZERZER
646          S(J)=ZERZER
647          SOS(J)=ZERZER
648          GO TO 33
649     4    CONTINUE
650          S(J)=SQRT (SMALLA(J)*ARG)*TWOTWO
651          SOS(J)=TENTEN*S(J)
652    33 CONTINUE
653       N1=1
654       SES = MAX (S(1),S(2),S(3),S(4),S(5),S(6))
655       DO 1131 J=1,6
656          JS  = SOS(J) + ONEONE
657          FJS = JS
658          STRUN(J) = FJS - ONEONE
659          IF ( S(J) .GT. ZERZER ) THEN
660             MM  = JA-IA(J)
661             EXPSAS=MIN(EXPMAX,MAX(EXPMIN,S(J)-SES))
662             SAS = EXP (EXPSAS)
663             EXPSUS=MIN(EXPMAX,MAX(EXPMIN,-S(J)))
664             SUS = EXP (EXPSUS)
665             EYE1(J) = ( TWOTWO * S(J)**2 -SIXSIX * S(J)
666      &              + SIXSIX + SUS * ( S(J)**2 - SIXSIX ) )
667      &              / ( EIGEIG * SMALLA(J)**2 )
668             IF ( J .EQ. 1 ) THEN
669                EYE0(J) = ( S(J) - ONEONE + SUS ) / ( TWOTWO*SMALLA(J) )
670 * Standard calculation
671                R   (J) = RMASS(MM)**2 * ALPH(MM) * ( EYE1(J) + BET(MM)
672      &                 * EYE0(J) ) * SAS
673             ELSE
674                R   (J) = CCOUL(J) * RMASS(MM)**2 * EYE1(J) * SAS
675             END IF
676             R (J) = MAX ( ZERZER, R (J) )
677             SIGMA = SIGMA + R (J)
678          END IF
679  1131 CONTINUE
680       NCOUNT = 0
681  6202 CONTINUE
682       IF(SIGMA)9,9,10
683  9    CONTINUE
684       DO 6100 J = 1,6
685          IF(JA-IA(J))6100,6101,6100
686  6101    CONTINUE
687          IF(JZ-IZ(J))6100,6102,6100
688  6100 CONTINUE
689       GO TO 6099
690  6102 CONTINUE
691       IF ( U .GT. ANGLGB ) GO TO 1000
692       JEMISS = J
693 C*****STORE,RESIDUAL NUC IS OF EMITTED PARTICLE TYPE
694 * If we are here this means that the residual nucleus is equal to
695 * one of the six emitted particle (the j-th one). So give to it
696 * all the energy, score it and return with 0 recoil and excitation
697 * energy for the residual nucleus
698       EPS = EREC
699       NPART(JEMISS) = NPART(JEMISS)+1
700       ITEMP=NPART(JEMISS)
701       NRNEEP = NRNEEP + 1
702       SMOM1(JEMISS) = SMOM1(JEMISS) + EPS
703       ITEMP=NPART(JEMISS)
704       EPART(ITEMP,JEMISS)=EPS
705       COSEVP(1,ITEMP,JEMISS) = COSLBR(1)
706       COSEVP(2,ITEMP,JEMISS) = COSLBR(2)
707       COSEVP(3,ITEMP,JEMISS) = COSLBR(3)
708       GO TO 8002
709 *
710 *-->-->-->-->--> go to return
711  6099 CONTINUE
712       IF(JA-8)72,51,72
713    51 CONTINUE
714       IF(JZ-4)72,1224,72
715 * Come here for a "normal" step
716    10 CONTINUE
717       LOPPAR = .FALSE.
718       CALL GRNDM(RNDM,1)
719       URAN=RNDM(1)*SIGMA
720       SUM = ZERZER
721       DO 16 J=1,6
722          K   = J
723          SUM = R(J)+SUM
724          IF ( SUM - URAN .GT. 0.D+00 ) GO TO 17
725    16 CONTINUE
726    17 CONTINUE
727       JEMISS=K
728       NPART(JEMISS) = NPART (JEMISS) + 1
729       JS = SOS (JEMISS) + ONEONE
730       IF(JS-1000)4344,4345,4345
731  4345 CONTINUE
732       RATIO2=(S(JEMISS)**3-6.D0*S(JEMISS)**2+15.D0*
733      &S(JEMISS)-15.D0)/((2.D0*S(JEMISS)**2-6.D0*S(JEMISS)+6.D0)*SMALLA
734      &(JEMISS))
735       GO TO 4347
736  4344 CONTINUE
737       RATIO2=(P2(JS)+(P2(JS+1)-P2(JS))*
738      &(SOS(JEMISS)-STRUN(JEMISS)))/SMALLA(JEMISS)
739  4347 CONTINUE
740       EPSAV=RATIO2*TWOTWO
741 *  +-------------------------------------------------------------------*
742 *  |  Neutron channel selected:
743       IF (JEMISS .EQ. 1) THEN
744          MM=JA-IA(J)
745          EPSAV=(EPSAV+BET(MM))/(ONEONE+BET(MM)*EYE0(JEMISS)
746      &        /EYE1(JEMISS))
747 *  |  +----------------------------------------------------------------*
748 *  |  |  Compute the fission width relative to the neutron one:
749 *  |  |  this part is taken from subroutine EMIT of LAHET
750          IF ( IFISS .GT. 0 .AND. JZ .GE. 71 .AND. .NOT. FISINH ) THEN
751 *  |  |  +-------------------------------------------------------------*
752 *  |  |  |  Compute the correction factor for the fission width:
753             IF ( JZ .GT. 88 ) THEN
754                AGOES = ONEONE
755 *  |  |  |
756 *  |  |  +-------------------------------------------------------------*
757 *  |  |  |
758             ELSE
759                AGOES = MAX ( ONEONE, ( U-SEVSEV ) / ( EPSAV+SEVSEV ) )
760             END IF
761 *  |  |  |
762 *  |  |  +-------------------------------------------------------------*
763 *  |  |  Finally this is the relative fission width:
764 *  |  |  This is : Probfs = 1 / ( 1 + G_n / G_f )
765             PROBFS = FPROB ( Z, A, U ) / AGOES
766 *  |  |  +-------------------------------------------------------------*
767 *  |  |  |  Check if it will be fission:
768             CALL GRNDM(RNDM,1)
769             IF ( RNDM (1) .LT. PROBFS ) THEN
770                FISINH = .TRUE.
771                KFISS  = 1
772 *  |  |  |  Undo the counting of the neutron evaporation
773                NPART (JEMISS) = NPART (JEMISS) -1
774                GO TO 9260
775             END IF
776 *  |  |  |
777 *  |  |  +-------------------------------------------------------------*
778          END IF
779 *  |  |
780 *  |  +----------------------------------------------------------------*
781       END IF
782 *  |
783 *  +-------------------------------------------------------------------*
784   20  CONTINUE
785       CALL GRNDM(RNDM,2)
786       E1=-HLFHLF*LOG(RNDM(1))
787       E2=-HLFHLF*LOG(RNDM(2))
788 * Eps should be the total kinetic energy in the CMS frame
789 * Standard calculation:
790       EPS=(E1+E2)*EPSAV+THRESH(JEMISS)-Q(JEMISS)
791       AR = A - IA(JEMISS)
792       ZR = Z - IZ(JEMISS)
793 * The CMS energy is updated
794       IMASS = NINT (AR)
795       IF ( IMASS .EQ. 8 .AND. NINT (ZR) .EQ. 4 ) THEN
796          UNEW = U - EPS - Q(JEMISS)
797          UMAX = U - THRESH(JEMISS)
798          IF ( UNEW .GT. 6.D+00 ) THEN
799             UMIN = 6.D+00
800          ELSE IF ( UNEW .GT. 4.47D+00 .AND. UMAX .GT. 6.D+00 ) THEN
801             UMIN = 4.47D+00
802             UNEW = 6.D+00
803          ELSE IF ( UNEW .GT. 1.47D+00 .AND. UMAX .GT. 2.94D+00 ) THEN
804             UMIN = 1.47D+00
805             UNEW = 2.94D+00
806          ELSE
807             UMIN = -0.1D+00
808             UNEW = ANGLGB * 0.1D+00
809          END IF
810       ELSE IF ( IMASS .LE. 4 ) THEN
811          IPRO = NINT ( ZR )
812          INEU = IMASS - IPRO
813          IF ( IMASS .EQ. 1 ) THEN
814 *  Be sure that residual neutrons or protons are not left excited
815             UMIN = 0.D+00
816             UNEW = 0.D+00
817             EPS  = U - Q(JEMISS)
818          ELSE IF ( IPRO .EQ. 0 .OR. INEU .EQ. 0 ) THEN
819 *  Ipro protons or ineu neutrons arrived here!
820             UMIN = CORRRR(JEMISS)
821             UNEW = U - EPS - Q(JEMISS)
822          ELSE IF ( IMASS .LE. 2 ) THEN
823 *  Be sure that residual deuterons are not left excited!
824             UMIN = 0.D+00
825             UNEW = 0.D+00
826             EPS  = U - Q(JEMISS)
827          ELSE IF ( ABS ( INEU - IPRO ) .LE. 1 ) THEN
828 *  For the moment also residual 3-H, 3-He and 4-He are not left
829 *  excited !
830             UMIN = 0.D+00
831             UNEW = 0.D+00
832             EPS  = U - Q(JEMISS)
833          ELSE
834             UMIN = CORRRR(JEMISS)
835             UNEW = U - EPS - Q(JEMISS)
836          END IF
837       ELSE
838          UMIN = CORRRR(JEMISS)
839          UNEW = U - EPS - Q(JEMISS)
840       END IF
841 * Standard calculation
842       IF(UNEW-UMIN)6200,6220,6220
843  6220 CONTINUE
844 *or   RNMASS = AR * AMUMEV + FKENER (AR,ZR)
845       RNMASS = AR * AMUMEV + FKENER (AR,ZR) + UNEW
846       UMIN2  = ( RNMASS + ZMASS (JEMISS) )**2
847       IF ( UMIN2 .GE. UMO2 ) THEN
848          GO TO 6200
849       END IF
850       U = UNEW
851 * C(i) are the direction cosines of the evaporated particle in the CMS
852 * frame, of course - C(i) are the ones of the residual nucleus
853       CALL RACO(C(1),C(2),C(3))
854 * Erncm, Eepcm are the total energies of the residual nucleus and
855 * of the evaporated particle in the CMS frame
856       ERNCM = 0.5D+00 * ( UMO2 + RNMASS**2 - Z2MASS (JEMISS) ) / UMO
857       EEPCM = UMO - ERNCM
858       PCMS  = SQRT ( EEPCM**2 - Z2MASS (JEMISS) )
859 * Now we perform the Lorentz transformation back to the original
860 * frame (lab frame)
861 * First the evaporated particle:
862       ETAX = ETACM * COSLBR (1)
863       ETAY = ETACM * COSLBR (2)
864       ETAZ = ETACM * COSLBR (3)
865       PCMSX = PCMS * C (1)
866       PCMSY = PCMS * C (2)
867       PCMSZ = PCMS * C (3)
868       ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
869       EPS = GAMCM * EEPCM + ETAPCM - ZMASS (JEMISS)
870       PHELP = ETAPCM / (GAMCM + 1.D0) + EEPCM
871       PLBPX = PCMSX + ETAX * PHELP
872       PLBPY = PCMSY + ETAY * PHELP
873       PLBPZ = PCMSZ + ETAZ * PHELP
874       PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
875       COSLBP (1) = PLBPX / PHELP
876       COSLBP (2) = PLBPY / PHELP
877       COSLBP (3) = PLBPZ / PHELP
878 * Then the residual nucleus ( for it c (i) --> - c (i) ):
879       EREC  = GAMCM * ERNCM - ETAPCM - RNMASS
880       EKRES = 1.D-03 * EREC
881       PHELP = - ETAPCM / (GAMCM + 1.D0) + ERNCM
882       PXRES = 1.D-03 * ( - PCMSX + ETAX * PHELP )
883       PYRES = 1.D-03 * ( - PCMSY + ETAY * PHELP )
884       PZRES = 1.D-03 * ( - PCMSZ + ETAZ * PHELP )
885       P2RES = PXRES * PXRES + PYRES * PYRES + PZRES * PZRES
886       PTRES = SQRT (P2RES)
887       COSLBR (1) = PXRES / PTRES
888       COSLBR (2) = PYRES / PTRES
889       COSLBR (3) = PZRES / PTRES
890 * Check energy and momentum conservation !!
891       IF (EREC .LE. 0.D+00) THEN
892          PTRES = 0.D+00
893          EREC  = 0.D+00
894       END IF
895 * Umo is the invariant mass of the system!!
896       UMO  = RNMASS
897       UMO2 = UMO * UMO
898       ELBTOT  = RNMASS + EREC
899       GAMCM   = ELBTOT / RNMASS
900       ETACM   = 1.D+03 * PTRES / RNMASS
901       GO TO 76
902  6200 CONTINUE
903       NCOUNT = NCOUNT + 1
904       IF ( NCOUNT .LE. 10 ) GO TO 20
905       SIGMA = SIGMA - R(JEMISS)
906 * if we are here we have sampled for > 10 times a negative energy Unew
907       NPART(JEMISS)=NPART(JEMISS)-1
908       R(JEMISS) = 0.D0
909       NCOUNT = 0
910       GO TO 6202
911    76 CONTINUE
912       JAT=JA-IA(JEMISS)
913       JZT=JZ-IZ(JEMISS)
914       IF(JAT-JZT)9,9,77
915    77 CONTINUE
916       JA=JAT
917       JZ=JZT
918 C*****STORE,END OF NORMAL CYCLE
919       SMOM1(JEMISS)=SMOM1(JEMISS)+EPS
920       ITEMP=NPART(JEMISS)
921       EPART(ITEMP,JEMISS)=EPS
922       COSEVP(1,ITEMP,JEMISS)=COSLBP(1)
923       COSEVP(2,ITEMP,JEMISS)=COSLBP(2)
924       COSEVP(3,ITEMP,JEMISS)=COSLBP(3)
925 * The following card switch to the rough splitting treatment
926       IF (JA .LE. 2) GO TO 1000
927       IF(JA-8)23,1223,23
928  1223 CONTINUE
929       IF(JZ-4)23,1224,23
930 * If we are here the residual nucleus is a 8-Be one, break it into
931 * two alphas with all the available energy (U plus the Q of the breakup)
932 * , score them and return with 0 recoil and excitation energy
933  1224 CONTINUE
934       LOPPAR = .FALSE.
935       IF(U)1228,1229,1229
936  1228 CONTINUE
937       EPS=0.D0
938       GO TO 1230
939  1229 CONTINUE
940  1230 CONTINUE
941       NBE8=NBE8+1
942 * C(i) are the direction cosines of the first alpha in the CMS
943 * frame, of course - C(i) are the ones of the other
944       CALL RACO(C(1),C(2),C(3))
945 * Ecms is the total energy of the alphas in the CMS frame
946       ECMS  = 0.5D+00 * UMO
947       PCMS  = SQRT ( ECMS**2 - Z2MASS (6) )
948 * Now we perform the Lorentz transformation back to the original
949 * frame (lab frame)
950 * First alpha:
951       ETAX = ETACM * COSLBR (1)
952       ETAY = ETACM * COSLBR (2)
953       ETAZ = ETACM * COSLBR (3)
954       PCMSX = PCMS * C (1)
955       PCMSY = PCMS * C (2)
956       PCMSZ = PCMS * C (3)
957       ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
958       EPS = GAMCM * ECMS + ETAPCM - ZMASS (6)
959       PHELP = ETAPCM / (GAMCM + 1.D0) + ECMS
960       PLBPX = PCMSX + ETAX * PHELP
961       PLBPY = PCMSY + ETAY * PHELP
962       PLBPZ = PCMSZ + ETAZ * PHELP
963       PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
964 * Store the first alpha!!
965       SMOM1(6) = SMOM1(6) + EPS
966       NPART(6) = NPART(6) + 1
967       ITEMP = NPART(6)
968       EPART (ITEMP,6) = EPS
969       COSEVP (1,ITEMP,6) = PLBPX / PHELP
970       COSEVP (2,ITEMP,6) = PLBPY / PHELP
971       COSEVP (3,ITEMP,6) = PLBPZ / PHELP
972 * Then the second alpha ( for it c (i) --> - c (i) ):
973       EPS = GAMCM * ECMS - ETAPCM - ZMASS (6)
974       PHELP = - ETAPCM / (GAMCM + 1.D0) + ECMS
975       PLBPX = - PCMSX + ETAX * PHELP
976       PLBPY = - PCMSY + ETAY * PHELP
977       PLBPZ = - PCMSZ + ETAZ * PHELP
978       PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
979 * Store the second alpha !!
980       SMOM1(6) = SMOM1(6) + EPS
981       NPART(6) = NPART(6) + 1
982       ITEMP = NPART(6)
983       EPART (ITEMP,6) = EPS
984       COSEVP (1,ITEMP,6) = PLBPX / PHELP
985       COSEVP (2,ITEMP,6) = PLBPY / PHELP
986       COSEVP (3,ITEMP,6) = PLBPZ / PHELP
987  8002 CONTINUE
988       LOPPAR = .FALSE.
989       EREC   = ZERZER
990       U      = ZERZER
991       EKRES  = ZERZER
992       PTRES  = ZERZER
993    72 CONTINUE
994       HEVSUM=SMOM1(3)+SMOM1(5)+SMOM1(6)+SMOM1(4)
995       RETURN
996 *.......................................................................
997 *///// RAL FISSION ROUTINE /////
998  9260 CONTINUE
999 *  +-------------------------------------------------------------------*
1000 *  |  Record the direction cosines of the fissioning nucleus
1001       DO 9270 I=1,3
1002          COSLF0 (I) = COSLBR (I)
1003  9270 CONTINUE
1004 *  |
1005 *  +-------------------------------------------------------------------*
1006       CALL FISFRA ( JA, JZ, U, EREC, UMO, GAMCM, ETACM )
1007 *  +-------------------------------------------------------------------*
1008 *  |  Check for fission failures!!
1009       IF ( .NOT. FISINH ) THEN
1010          PENBAR = .FALSE.
1011          GO TO 23
1012       END IF
1013 *  |
1014 *  +-------------------------------------------------------------------*
1015 *  Do not pick up the fission fragments, rather go back to Evevap
1016       HEVSUM=SMOM1(3)+SMOM1(5)+SMOM1(6)+SMOM1(4)
1017       RETURN
1018       END