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