GEANT321/fluka/nucriv.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:19:57  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.43  by  S.Giani
  11 *-- Author :
  12 *$ CREATE NUCRIV.FOR
  13 *COPY NUCRIV
  14 *
  15 *=== nucriv ===========================================================*
  16 *
  17       SUBROUTINE NUCRIV (ITTTT,ELAB,CX,CY,CZ,BBTAR,ZZTAR,RHOO)
  18
  19 #include "geant321/dblprc.inc"
  20 #include "geant321/dimpar.inc"
  21 #include "geant321/iounit.inc"
  22 *
  23 *----------------------------------------------------------------------*
  24 *                                                                      *
  25 *  Last change  on  16-sep-92    by    Alfredo Ferrari, Infn - Milan   *
  26 *                                                                      *
  27 *  Nucrin90 by A. Ferrari: tentative version to overcome many troubles *
  28 *                          in energy, momentum and charges conservation*
  29 *                                                                      *
  30 *----------------------------------------------------------------------*
  31 *
  32 C
  33 C--------------------------------------------------
  34 C*** FINUC: NUCRIN FINAL STATE PARTICLE LIST WITH KINEM. VARIABLES
  35 C*** FINLSP HADRIN FINAL STATE PARTICLE LIST WITH KINEM. VARIABLES
  36 C*** (NUMBER OF PARTICLES,PARTICLE TYPE INDEX,DIRECTION COSINES,ENERGY
  37 C*** ABSOLUTE MOMENTUM AND IN FINUC IN ADDIT. EXCITATION ENERGY
  38 C--------------------------------------------------
  39 #include "geant321/balanc.inc"
  40 #include "geant321/corinc.inc"
  41 #include "geant321/finuc2.inc"
  42 #include "geant321/finlsp2.inc"
  43 #include "geant321/nucdat.inc"
  44 #include "geant321/resnuc.inc"
  45 *     PARAMETER ( TVEPSI = 1.D-05 )
  46       COMMON / FKNUCF / DELEFT, EKRECL, V0EXTR, ITTA, ITJ, LVMASS
  47       DIMENSION EXSOP (2)
  48       REAL RNDM(1),RANGS1,RANGS2
  49       COMMON /FKFERM/ELABKE,A,TPNVK,ELATE,IFERT
  50 C
  51 C--------------------------------------------------
  52 C***
  53 C***
  54 C***  SAMPLING OF A HADRON-NUCLEUS COLLISION EVENT
  55 C***
  56 C   ITTTT - TYPE OF INCOMING HADRON
  57 C
  58 C--------------------------------------------------
  59 C***  POSSIBLE: N, AN,PI,K,AK,Y,AY BY THE INDEX ITTTT=IT=1,2,8-25
  60 C   PLAB - MOMENTUM OF INCOMING HADRON
  61 C
  62 C***  RANGE UP TO 5. GEV/C (FROM ABOUT 0. ... 0.1 GEV/C)
  63 C   ELAB - TOTAL ENERGY OF INCOMING HADRON IN GEV
  64 C   CX,CY,CZ - DIRECTION COSINES
  65 C***  ANUC,ZNUC = NUMBERS OF NUCLEONS AND PROTONS
  66 C*** RHOO = MATERIAL DENSITY (G/CM**3)
  67 C
  68 C*** FINAL STATE PARTICLE CHARACTERISTICS TABLE IN /FINUC/:
  69 C*** IRN   - NUMBER OF FINAL STATE PARTICLES
  70 C*** ITRN(I) - FINAL STATE PARTICLE TYPE INDEX
  71 C*** CXRN,CYRN,CZRN (I) - DIRECTION COSINES OF F.S.P. (LAB SYST.)
  72 C*** ELR,PLR (I) - LAB.ENERGY AND MOMENTUM OF F.S.P. (GEV, GEV/C)
  73 C*** TV    - EXCITATION ENERGY (GEV)
  74 C
  75 C--------------------------------------------------
  76       DIMENSION THRESR(30)
  77       SAVE THRESR,ITPRF,IAMC,JAMC,INS,IXPI,JNUC
  78       DATA THRESR/1.9D0, 0.D0, 5*9.D0, 1.9D0, 0.D0, 3*9.D0, 1.08D0,
  79      * 1.08D0, 1.44D0, 1.08D0, 6*9.D0, 1.08D0,
  80      * 1.44D0, 1.08D0, 5*9.D0/
  81       DIMENSION IXPI (30)
  82 *
  83 *----------------------------------------------------------------------*
  84 *                                                                      *
  85 *     Ixpi :                                                           *
  86 *             = -1 for antinucleons ( pbar, nbar, and for Lambdabar,   *
  87 *                  see below )                                         *
  88 *             =  0 for pi-, K- ( stopping particles )                  *
  89 *             =  1 for other mesons and particles (p, n, pi+, pi0, K+, *
  90 *                  K0, K0bar )                                         *
  91 *             =  2 for leptons (e-, e+, nu, nubar, mu+, mu-, photons,  *
  92 *                  Klong, Kshort )                                     *
  93 *             =  4 for hyperons ( Lambda, Lambdabar, Sigma-, Sigma+,   *
  94 *                  Sigma0 ). Actually for Lambdabar the inxpi flag in  *
  95 *                  the code, which corresponds to the ixpi one of the  *
  96 *                  projectile, is then set to - 1 as for antinucleons  *
  97 *                                                                      *
  98 *----------------------------------------------------------------------*
  99 *
 100       DATA IXPI/1,-1,5*2,1,-1,3*2,1,0,1,0,4,4,2,3*4,1,1,1,5*2/
 101       COMMON /FKPERC/ IPERCO
 102       COMMON /FKENCO/ ETEST,TNKTE
 103 C
 104 C--------------------------------------------------
 105 C***  PARTICLE CHARACTERISTICS: MASSES, DECAY WIDTH, LIFE TIME,ELECT.AND
 106 C***  BARYONIC CHARGE, DECAY CHANNEL INDICEES
 107 C--------------------------------------------------
 108 *
 109 *  /Abltis/ common is initialized in Hadden. Masses are the same as in
 110 *  /Part/ common at least up to particle n. 94, the same for the baryonic
 111 *  charge and the electric charge : what about particle 26????????????
 112 *
 113       COMMON / FKABLT / AM (110), GA (110), TAU (110), ICH (110),
 114      &                  IBAR (110), K1(110), K2(110)
 115       COMMON / FKNUCT / ETHR, PTHR
 116       DIMENSION INS(30)
 117       DIMENSION AMHH(15),IAMC(15),JAMC(30)
 118       DATA IAMC/1,2,8,9,12,15,16,17,18,19,20,21,22,24,25/
 119       DATA JAMC/2*1,5*0,2*1,2*0,1,2*0,16*0/
 120       DATA INS/13,13,5*32,14,14,3*32,10,10
 121      *,12,11,7*32,15,15,5*32/
 122       DATA JNUC/0/
 123       DIMENSION ITPRF(110)
 124       DATA ITPRF/-1,-1,5*1,-1,-1,1,1,1,-1,-1,-1,-1,6*1,-1,-1,-1,85*1/
 125       LOGICAL LEMINU, LPRONU, LPCASC, LNCASC, LCORIN, LDELTX, LVMASS,
 126      &        LHYPER
 127 *
 128       LCORIN = .FALSE.
 129       LVMASS = .FALSE.
 130       ANUC = BBTAR
 131       ZNUC = ZZTAR
 132       TPNVK  = 0.D+00
 133       IRNTOT = 0
 134       INNUC=1
 135       ITNUCR=ITTTT
 136       SICO   = 1.D0
 137       ELPH   = 0.D0
 138       PLABCO = 1.D0
 139 C
 140 C--------------------------------------------------
 141 C*** SICO,PLABCO - EFFECTIVE CROSS SECTION- AND EFFECT. LABMOMENTUM
 142 C***      CORRECTION FACTORS FOR Y-A- AND AY-A-COLLISIONS
 143 C*** ITNUCR - HYPERON NUCLEUS COLL. PARTICLE TYPE FOR USE IN NIZL,
 144 C***      IT - USED IN HADRIN
 145 C--------------------------------------------------
 146       IBAT = 0
 147 C--------------------------------------------------
 148 C*** INXPI:
 149 C*** ORDERING INDICEES FOR STOPPING PARTICLES (0), OTHER MESONS (1),
 150 C*** ANTINUCLEONS (-1), HYPERONS (4),LEPTONS,KO,AKO(2),OTHER PART. (1)
 151 C
 152 C--------------------------------------------------
 153 C*** INCLUSION OF HYPERON-NUCLEUS-COLLISIONS
 154 C*** PRELIMINARY FOR SIGMA+ NO STRANGENES-CONSERVATION
 155 C
 156 C--------------------------------------------------
 157 *
 158 *  From here it is the number of the incoming nucleon
 159 *
 160       IT = ITTTT
 161 *
 162 *  Nxpi = 0 means hadrin call is possible
 163 *
 164       NXPI  = 0
 165       INXPI = IXPI(IT)
 166       IF ( IT .EQ. 23 .AND. ELAB - AM(IT) .LE. 0.05D+00 ) INXPI = 0
 167 *  +-------------------------------------------------------------------*
 168 *  |
 169       IF ( INXPI .EQ. 4 ) THEN
 170          CALL HYPERO ( IT, ITNUCR, SICO, PLABCO )
 171          LHYPER = .TRUE.
 172 *  |
 173 *  +-------------------------------------------------------------------*
 174 *  |
 175       ELSE
 176          LHYPER = .FALSE.
 177       END IF
 178 *  |
 179 *  +-------------------------------------------------------------------*
 180 C
 181 C--------------------------------------------------
 182 C*** CUT OFF ENERGY CONSTANTS:
 183 C--------------------------------------------------
 184       ETHR  = 0.001D+00
 185 C
 186 C--------------------------------------------------
 187 C*** CDDT,CEET PARAMETERS FOR GAUSSIAN WIDTH'
 188 C--------------------------------------------------
 189       N = ITTTT
 190 *
 191 *  Particle 18 is Alambda
 192 *
 193       IF ( ITTTT .EQ. 18 ) INXPI = -1
 194 C
 195 C--------------------------------------------------
 196 C*** ETEST = ENERGY CONSERVATION TEST VARIABLE (SHOULD BE 0 AT RETURN)
 197 C*** TLAB = KINETIC ENERGY
 198 C--------------------------------------------------
 199       JJ    = JNUC
 200       CDDT  = 0.5D0
 201       IHACA = 0
 202       IN    = INS(N)
 203 C
 204 C--------------------------------------------------
 205 C*** JJ=1:
 206 C*** OPTION 1: GO TO NIZL, CROSS SECTION CALCULATION
 207 C*** JJ=2:
 208 C*** OPTION 2: EVENT SAMPLING
 209 C--------------------------------------------------
 210       TLAB  = ELAB - AM (IT)
 211 *  +-------------------------------------------------------------------*
 212 *  |  Decide which is the target nucleon ( flag Itta ) : a proton .....
 213       CALL GRNDM(RNDM,1)
 214       IF ( RNDM (1) .LT. ZNUC / ANUC ) THEN
 215          ITTA = 1
 216          ITJ  = 1
 217 *  |
 218 *  +-------------------------------------------------------------------*
 219 *  |                                            .... or a neutron
 220       ELSE
 221          ITTA = 8
 222          ITJ  = 2
 223       END IF
 224 *  |
 225 *  +-------------------------------------------------------------------*
 226       V0OLD = V0WELL (ITJ)
 227       PLAB = SQRT ( TLAB * ( ELAB + AM (IT) ) )
 228 *  +-------------------------------------------------------------------*
 229 *  |
 230       IF ( ABS ( PLAB - 5.D0 ) .GE. 4.99999D0 ) THEN
 231          WRITE(LUNOUT,99996) PLAB
 232          WRITE(LUNERR,99996) PLAB
 233 *        STOP         Commented out A. Fasso' 1989
 234          IRN = NP0
 235          RETURN
 236 99996    FORMAT (3(5H ****/),' Nucrin: projectile momentum',
 237      &          ' outside of the allowed region, plab, kproj:',
 238      *   1P,E15.5,3X,0P,I3,/,3(5H ****/))
 239       END IF
 240 *  |
 241 *  +-------------------------------------------------------------------*
 242 *     IF ( JJ .EQ. 1 ) GO TO 1000
 243 *-->-->-->-->--> go to cross section calculation if this option was se-
 244 *                lected
 245 C
 246 C--------------------------------------------------
 247 C*** IF LOW KINETIC ENERGY ( TLAB < ETHR ): ONLY EXCITATION
 248 C*** STOPPING PARTICLES
 249 C*** AND ANNIHILATION , INXPI=IXPI(IT) LE 0
 250 C--------------------------------------------------
 251 *  +-------------------------------------------------------------------*
 252 *  |  Inxpi > 0: p, n, pi+, K+, Lambda, Sigma+, Sigma-, Sigma0
 253 *  |  ( no lepton or pi0 projectile of course, now also pi0
 254 *  |    are possible projectiles )
 255       IF (INXPI .GT. 0) THEN
 256 C
 257 C--------------------------------------------------
 258 C*** ETHRR=THRESHOLD VALUE, FOR TLAB<ETHRR IS NO H-N COLLISION POSSIBLE,
 259 C*** ONLY EXCITATION + CASCADE AND/OR ANNIHILATION CORRESP TO THE HADRN
 260 *  |  +----------------------------------------------------------------*
 261 *  |  |  Raised to 10 MeV: A. Ferrari, 11-apr-1991
 262          IF (TLAB .LE. 10.D+00 * ETHR) THEN
 263 *  |  |   Here particles not stopping and not to be annihilated below
 264 *  |  |   threshold
 265             ANOW  = ANOW  + IBAR (IT)
 266             KTARN = KTARN - IBAR (IT) + ICH (IT)
 267             ZNOW  = ZNOW  + ICH  (IT)
 268             KTARP = KTARP - ICH  (IT)
 269             AMNRES = AMUAMU * ANOW + 1.D-03 * FKENER ( ANOW, ZNOW ) -
 270      &               ZNOW * AMELEC + ELBNDE ( NINT (ZNOW) )
 271             TV  = ETTOT - AMNRES
 272             IRN = NP0
 273             RETURN
 274          END IF
 275 *  |  |
 276 *  |  +----------------------------------------------------------------*
 277 *  |
 278 *  +-------------------------------------------------------------------*
 279 *  |    Here for ap, an, pi-, k-, alambda
 280       ELSE
 281 *  |  +----------------------------------------------------------------*
 282 *  |  |  Threshold is raised for particles with ixpi < 0 (ap, an,
 283 *  |  |  alambda) to 80 MeV
 284          IF (INXPI .LT. 0) THEN
 285             ETHRR = 0.08D+00
 286 *  |  |
 287 *  |  +----------------------------------------------------------------*
 288 *  |  |  for pi-, K- the threshold is Ethr = 1 MeV
 289          ELSE
 290             ETHRR = ETHR
 291          END IF
 292 *  |  |
 293 *  |  +----------------------------------------------------------------*
 294 C
 295 C--------------------------------------------------
 296 C*** UP TO ETHRR=.08 GEV NO AN, AP IN FINAL STATE
 297 C--------------------------------------------------
 298          AIO  = 10.D0
 299          XPI  = 1.D0 - AIO * ( TLAB - ETHRR )
 300 *  |  Patch for pseudo pi0
 301          IF ( IT .EQ. 23 ) XPI = 1.D+00
 302 *  |  For pi- and K- :
 303 *  |                   Xpi < 0 for Tlab > 1/Aio + Ethrr = 101 MeV
 304 *  |                   Xpi > 1 for Tlab < Ethrr = 1 MeV
 305 *  |  For antibaryons:
 306 *  |                   Xpi < 0 for Tlab > 1/Aio + Ethrr = 180 MeV
 307 *  |                   Xpi > 1 for Tlab < Ethrr = 80 MeV
 308 *  |  +----------------------------------------------------------------*
 309 *  |  |   Nxpi=0 means hadrin-call is possible, else impossible
 310 *  |  |   So Hadrin call is always possible for pi- and K- if
 311 *  |  |   Tlab > 101 MeV, and for antibaryons if Tlab > 180 MeV.
 312 *  |  |   Hadrin call is always impossible for pi- and K- if
 313 *  |  |   Tlab < 1 MeV (never possible) and for antibaryons if
 314 *  |  |   Tlab < 80 MeV. Actually Nxpi > 1 for antibaryons means forced
 315 *  |  |   annihilation, but Hadrin can be called the same, even though
 316 *  |  |   only final states with annihilation will be accepted
 317          CALL GRNDM(RNDM,1)
 318          IF ( RNDM(1) .GT. XPI ) THEN
 319             NXPI = 0
 320 *  |  |                       particle non stopping (annihilating)
 321 *  |  +----------------------------------------------------------------*
 322 *  |  |                       particle stopping (annihilating)
 323          ELSE
 324             NXPI = 1
 325          END IF
 326 *  |  |
 327 *  |  +----------------------------------------------------------------*
 328       END IF
 329 *  |
 330 *  +-------------------------------------------------------------------*
 331 *  +-------------------------------------------------------------------*
 332 *  |    Incident pi- or K- stopping in the nucleus, no possibility
 333 *  |    to call hadrin: it must be updated using the new stopping
 334 *  |    particle module as soon as it will be available!!
 335       IF ( NXPI .GT. 0 .AND. IBAR (IT) .EQ. 0 ) THEN
 336 *  |  Give the total lab energy of stopping pi-, K- into cascade and
 337 *  |  exit: label 405 is where this is performed.
 338 *  |  Before change a proton into a neutron to account for charge
 339 *  |  conservation and adjust the total energy with the mass
 340 *  |  difference (Note that for K- no strangeness conservation occurs
 341          IF ( IT .NE. 23 ) THEN
 342             ZNOW  = ZNOW - 1.D+00
 343             KTARP = KTARP + 1
 344             KTARN = KTARN - 1
 345          END IF
 346          TPK  = 0.D0
 347          TNK  = 0.D0
 348          TV   = 0.D0
 349          TOEFF = ELAB
 350          LCORIN = .TRUE.
 351          CALL CORSTP ( TOEFF )
 352          CALL GRANOR(RANGS1, RANGS2)
 353          TNKTE = MAX ( 0.003D+00, EKMNNU (2)
 354      &         * ( 0.5D+00 - 0.25D+00 * RANGS1 ) )
 355          TPKTE = MAX ( 0.003D+00, EKMNNU (1)
 356      &         * ( 0.5D+00 - 0.25D+00 * RANGS2 ) )
 357          GO TO 405
 358 *  |
 359 *  +-->-->-->-->--> go to the cascade simulation
 360       END IF
 361 *  |
 362 *  +----------------------------------------------------------------*
 363 C
 364 C--------------------------------------------------
 365 C*** ELABKE = INVARIANT KINETIC ENERGY OF THE H-A-SYSTEM+PROJ.H.MASS
 366 C*** A = TOTAL ENERGY OF THE PRIMARY PARTICLE IN LABSYSTEM
 367 C*** IFORBI,IFERT LOOP COUNTING INDICEES FOR CHOICE OF FERMION MOMENTUM
 368 C*** AMNTAR - NUCLEUS MASS
 369 C*** UMOJAN - H-A-C.M.S.-ENERGY
 370 C*** ELABKE - KINETIC H-A-ENERGY IN C.M.S
 371 C--------------------------------------------------
 372 *
 373       UMOJAN = SQRT ( AM(IT)**2 + AMNTAR**2 + 2.D+00 * ELAB * AMNTAR )
 374 C
 375 C--------------------------------------------------
 376 C***  START WITH THE EVENT SIMULATION
 377 C--------------------------------------------------
 378       ELABKE = UMOJAN - AMNTAR
 379       A = ELAB
 380       IFORBI = 0
 381       IFERT  = 0
 382 *  Call Ferset for setting up anything for the Fermi motion
 383       CALL FERSET
 384       AMMIT = AM (IT)
 385       AMMTA = AM (ITTA)
 386       AIT2  = AMMIT**2
 387       ATA2  = AMMTA**2
 388 *  Ecmo is the square of the invariant mass of the projectile-nucleon
 389 *  system
 390 *     ECMO = AIT2 + ATA2 + 2.D+00 * AMMTA * ELAB
 391 *  Put here the correct value for Ecmo taking into account
 392 *  Fermi motion etc.
 393 *  Compute Ekrecl in the most favourable situation, with final
 394 *  excitation energy zero (tvepsi) and all the energy given to the
 395 *  projectile
 396 *  Dex = - Tveuz, Ekres = sqrt(amnres**2 + pfrmi**2) - amnres,
 397 *  Ekrecl = Ekres + Dex
 398 *  **** Now the checks are done using atomic masses directly **** *
 399       EKREC0 = EKRECL
 400       TVEUZ0 = TVEUZ
 401       PFRMSQ = PXFRM**2 + PYFRM**2 + PZFRM**2
 402       EKRES  = SQRT ( ( AMMRES + 2.D+00 * TVEPSI )**2 + PFRMSQ ) -
 403      &         AMMRES - 2.D+00 * TVEPSI
 404 *  -TVEUZ to reduce the excitation energy not the recoil one
 405 * (it would have been more exact to set
 406 *  Efrm = Efrm + Tveuz, but it is better to conserve Efrm
 407 *  coherent with the original value)
 408       EKRECL = EKRES - TVEUZ + 2.D+00 * TVEPSI
 409       EZERO  = EFRM - V0WELL (ITJ) - EKRECL - EBNDNG (ITJ)
 410       ECHCK  = ELAB + EZERO
 411       PLABSQ = PLAB * PLAB
 412       PFRSCA = PXFRM * CX + PYFRM * CY + PZFRM * CZ
 413       ECMO   = ECHCK**2 - PLABSQ - PFRMSQ - 2.D+00 * PLAB * PFRSCA
 414       UMIN2  = ( AMMIT + AMMTA )**2
 415 *  +----------------------------------------------------------------*
 416 *  |   Check if the Hadrin call is energetically possible, else
 417 *  |   give the total energy to cascade and excitation:
 418       IF ( ECMO .LT. UMIN2 .AND. NXPI .LE. 0 ) THEN
 419          V0WELL (ITJ) = V0WELL (ITJ) - V0EXTR
 420          ANOW  = ANOW  + IBAR (IT)
 421          KTARN = KTARN - IBAR (IT) + ICH (IT)
 422          ZNOW  = ZNOW  + ICH  (IT)
 423          KTARP = KTARP - ICH  (IT)
 424          TPK  = 0.D+00
 425          TNK  = 0.D+00
 426          TV   = 0.D+00
 427          TPNVK  = 0.D+00
 428          TOEFF  = ELAB - IBAR (IT) * AMMIT
 429          LCORIN = .TRUE.
 430          CALL CORSTP ( TOEFF )
 431          CALL GRANOR(RANGS1, RANGS2)
 432          TNKTE = MAX ( 0.003D+00, EKMNNU (2)
 433      &         * ( 0.5D+00 - 0.25D+00 * RANGS1 ) )
 434          TPKTE = MAX ( 0.003D+00, EKMNNU (1)
 435      &         * ( 0.5D+00 - 0.25D+00 * RANGS2 ) )
 436          GO TO 405
 437 *  |
 438 *  |-->-->-->-->--> go to the cascade simulation
 439       END IF
 440 *  |
 441 *  +-------------------------------------------------------------------*
 442       EKRECL = EKREC0
 443 *  ******  The check for peripheral collisions is now moved here: *****
 444 *   The old comment was:
 445 *   " I.E.,H-COLLISION WITH N IN THE MOST OUTSIDE NUCLEUS SHELL
 446 *     WITH THE THICKNESS OF ONE NUCLEON RADIUS
 447 *     RATIO OF THE SHELL VOLUME TO NUCLEUS VOLUME IS ELANU3
 448 *     THE H-N-COLL. WILL BE SAMPLED FIRSTLY WITH THE TOTAL PRIMARY
 449 *     ENERGY ELAB, THE REMAINING ENERGY IS USED IN THE ABOVE SAMPLED
 450 *     RATIO FOR CASCADE AND EXCITATION SWITCH VARIABLE IS IPERCO=1 "
 451 *   The old coding was: ....
 452 *     ELRAN  = RNDM (V)
 453 *     ELANU3 = ANUC**(-0.6666666666666667D+00)
 454 *     IF ( ELRAN .LE. ELANU3 .AND. ABS (ELPP-T3) .GT. 0.0001D0)
 455 *    ..... but
 456 *    in my opinion we must perform peripheral collisions checking
 457 *    the ratio between the total area (r**2) and the circular
 458 *    corona corresponding to 1 nucleon radius, so checking:
 459 *    (now to switch off peripheral collisions we need Elanu3 = 1)
 460 *    Furthermore no cascade and excitation energy is computed since
 461 *    it is extremely likely that the remaining energy (not more than
 462 *    Tveuz) will go thouroly in excitation one: anyway it is divided
 463 *    among the 3 contributions in the usual way
 464 *     ELRAN  = RNDM (V)
 465 *     ELANU3 = ( ( ANUC - 1.D+00 ) / ANUC )**0.6666666666666667D+00
 466 *    We changed again to the original idea, which can be geometrically
 467 *    explained as a corona of dr = r0 , corresponding roughly to
 468 *    1 nucleon. Anyway we need a much more detailed insight of this
 469 *    problem with possibly checks with experimental data
 470       CALL GRNDM(RNDM,1)
 471       ELRAN  = RNDM (1)
 472       ELANU3 = 1.D+00 - 1.D+00 / ANUC**0.6666666666666667D+00
 473 *  +-------------------------------------------------------------------*
 474 *  |  Check for peripheral collisions: note that in this way peripheral
 475 *  |  collisions are switched off for stopping antibaryons!!
 476       IF ( ELRAN .GE. ELANU3 .AND. NXPI .EQ. 0 ) THEN
 477          IPERCO = 1
 478          TOEFF  = ELAB - AMMIT
 479          LCORIN = .TRUE.
 480          CALL CORSTP ( TOEFF )
 481          CALL GRANOR(RANGS1, RANGS2)
 482          TNKTE = MAX ( 0.003D+00, EKMNNU (2)
 483      &         * ( 0.5D+00 - 0.25D+00 * RANGS1 ) )
 484          TPKTE = MAX ( 0.003D+00, EKMNNU (1)
 485      &         * ( 0.5D+00 - 0.25D+00 * RANGS2 ) )
 486          EKRECL = EKREC0
 487          ELPP   = ELAB
 488          PLPP   = PLAB
 489          DEX    = 0.D+00
 490          PSPESQ = PLABSQ + PFRMSQ + 2.D+00 * PLAB * PFRSCA
 491          IUMO   = 0
 492 *  |  +----------------------------------------------------------------*
 493 *  |  |
 494 1500     CONTINUE
 495             IUMO   = IUMO   + 1
 496             EEX    = TVEUZ  + DEX
 497 *  |  |  Now for iperco we are working with atomic masses
 498 *           AMSTAR = AMNRES + EEX
 499             AMSTAR = AMMRES + EEX
 500             EKRES  = SQRT ( AMSTAR**2 + PFRMSQ ) - AMSTAR
 501 *  |  |   Note +DEX to reduce actually the excitation energy, not
 502 *  |  |   the recoil one (it would have been more exact to set
 503 *  |  |   Efrm = Efrm - Dex, but it is better to conserve Efrm
 504 *  |  |   coherent with the original value)
 505             EKRECL = EKRES + DEX
 506             EZERO  = EFRM  - V0WELL (ITJ) - EKRECL - EBNDNG (ITJ)
 507             ESPENT = ELPP  + EZERO
 508 *           ELEFT  = ETTOT - ESPENT
 509             ELEFT  = ETTOT - ESPENT + DELEFT
 510             UMO2   = ESPENT**2 - PSPESQ
 511             DELTU2 = UMIN2 - UMO2
 512 *  |  |  +-------------------------------------------------------------*
 513 *  |  |  |  We need more invariant mass for the collision!!!
 514             IF ( DELTU2 .GT. 0.D+00 ) THEN
 515 *  |  |  |  +----------------------------------------------------------*
 516 *  |  |  |  |
 517                IF ( IUMO .GT. 3 ) THEN
 518                   V0WELL (ITJ) = V0WELL (ITJ) - V0EXTR
 519                   ANOW  = ANOW  + IBAR (IT)
 520                   KTARN = KTARN - IBAR (IT) + ICH (IT)
 521                   ZNOW  = ZNOW  + ICH  (IT)
 522                   KTARP = KTARP - ICH  (IT)
 523                   GO TO 405
 524                END IF
 525 *  |  |  |  |
 526 *  |  |  |  +----------------------------------------------------------*
 527                DELTEX = 0.5D+00 * DELTU2 * ELEFT / ( ESPENT * AMSTAR )
 528                DELTEX = MIN ( DELTEX, EEX )
 529                DEX    = DEX - DELTEX
 530             GO TO 1500
 531 *  |  |  |
 532 *  |  +-<|--<--<--<--<--< go to compute the new invariant mass!!
 533 *  |  |  |
 534 *  |  |  +-------------------------------------------------------------*
 535 *  |  |  |
 536             ELSE
 537                TVEUZ = EEX
 538             END IF
 539 *  |  |  |
 540 *  |  |  +-------------------------------------------------------------*
 541 *  |  |
 542 *  |  +----------------------------------------------------------------*
 543 *  |  we are ready for the interaction !!!!!
 544          IHACA = 0
 545 *  |  +----------------------------------------------------------------*
 546 *  |  |
 547 1600     CONTINUE
 548             ITFRHD = IT
 549             IHACA  = IHACA + 1
 550             IF ( PLPP .LT. 1.D-04 ) THEN
 551                WRITE (LUNERR,*)
 552      &              ' Nucriv: iperco=1,plpp,elpp,it',
 553      &                PLPP,ELPP,IT
 554                WRITE (LUNERR,*)'      IHACA,LVMASS,ELAB,PLAB',
 555      &                                IHACA,LVMASS,ELAB,PLAB
 556             END IF
 557             CALL FERHAV ( IT, ELPP, PLPP, CX, CY, CZ )
 558 *  |  |  +-------------------------------------------------------------*
 559 *  |  |  | Check if we have reached the maximum numbers of calls
 560             IF ( IHACA .GT. 15 ) THEN
 561                GO TO 1700
 562 *  |  |  |-->-->-->-->--> give up !!!
 563             END IF
 564 *  |  |  |
 565 *  |  |  +-------------------------------------------------------------*
 566 *  |  | Ir is the number of secondaries produced in the
 567 *  |  | Ferhad/Hadrin call
 568 *  |  |  +-------------------------------------------------------------*
 569 *  |  |  |
 570             IF ( IR .EQ. 0 ) THEN
 571                V0WELL (ITJ) = V0WELL (ITJ) + V0EXTR
 572                GO TO 1600
 573 *  |  |  |
 574 *  |  +-<|-<--<--<--<--< go to resampling the interaction
 575             END IF
 576 *  |  |  |
 577 *  |  |  +-------------------------------------------------------------*
 578 *  |  | Of course ir=1 means that the only outgoing particle is
 579 *  |  | still the projectile !! If we have an incoming antibaryon
 580 *  |  | we need at least 2 particles in the final state (why???)
 581 *  |  |  +-------------------------------------------------------------*
 582 *  |  |  |
 583             IF ( IR .LT. 2 .AND. INXPI .LT. 0 ) THEN
 584                V0WELL (ITJ) = V0WELL (ITJ) + V0EXTR
 585                GO TO 1600
 586 *  |  |  |
 587 *  |  +-<|-<--<--<--<--< go to resampling the interaction
 588             END IF
 589 *  |  |  |
 590 *  |  |  +-------------------------------------------------------------*
 591 1700     CONTINUE
 592          TVEUZ = TVEUZ + EKRECL
 593          TV  = TVEUZ
 594          TPK = 0.D+00
 595          TNK = 0.D+00
 596          GO TO 2800
 597 *  |  |                           End loop for Ferhad/hadrin call
 598 *  |  +----------------------------------------------------------------*
 599 *  |    Yes peripheral collison
 600 *  +-------------------------------------------------------------------*
 601 *  |    No peripheral collision
 602       ELSE
 603          IPERCO = -1
 604       END IF
 605 *  |
 606 *  +-------------------------------------------------------------------*
 607       I1651 = 0
 608 *  +-------------------------------------------------------------------*
 609 *  |  It is possible to come here for resampling !!!!
 610 1651  CONTINUE
 611          TVEUZ = TVEUZ0
 612          IRN = NP0
 613          DDT = CDDT * TLAB
 614          EMD = 0.D+00
 615 C
 616 C--------------------------------------------------
 617 C*** EFFECTIVE COLLISION ENERGY CALCULATION FOR USE IN CASCADE
 618 C*** AND EXCITATION ENERGY CALCULATION FOR ANNIHILATION
 619 C--------------------------------------------------
 620 C
 621 C--------------------------------------------------
 622 C*** I653 LOOP COUNTER INTEGER FOR BAD ENERGY CORRECTION IN ANNIHILATION
 623 C--------------------------------------------------
 624          I653 = 0
 625 *  |  +----------------------------------------------------------------*
 626 *  |  |  It is possible to come here for resampling !!!!
 627  653     CONTINUE
 628 C
 629 C--------------------------------------------------
 630 C*** TOEFF = EFFECTIVE KINETIC COLLISION ENERGY
 631 C--------------------------------------------------
 632             I653 = I653 + 1
 633             TOEFF= TLAB
 634 C
 635 C--------------------------------------------------
 636 C*** TOEFF ONLY FOR ANNIHILATION NE TLAB POSSIBLE, ELSE JUST SELECTION
 637 C*** OF ENERGY FRACTIONS
 638 C--------------------------------------------------
 639 *  |  |  +-------------------------------------------------------------*
 640 *  |  |  |   If we are here with a pi- or a K- this means Nxpi = 0,
 641 *  |  |  |   ==> no stopping (annihilating) particle but Hadrin call
 642 *  |  |  |   for pi-, K-, while for ap, an, alambdas may be that Nxpi=1
 643 *  |  |  |   ==> no Hadrin call (more exactly, Hadrin can be called but
 644 *  |  |  |   the particle is stopping ( annihilating ) )
 645             IF ( INXPI .LE. 0 ) THEN
 646 *  |  |  |   Note, Ibar (itta) = 1 of course, Ibar (it) maybe 0 or -1
 647 *  |  |  |   so Ibl = 0 for pi-, K- and Ibl = 1 for ap, an, alambdas
 648                IBL = IBAR (ITTA) * IBAR (IT)
 649                IBL = ABS  (IBL)
 650 *  |  |  |   Icl = -1 for pi-, K-, ap, Icl = 0 for an and alambdas
 651                ICL = ICH (IT)
 652 *  |  |  |   ( Ibar (itta) - Ibar (it) ) * Ibl = 2 for ap, an, alambdas
 653 *  |  |  |   ( Ibar (itta) - Ibar (it) ) * Ibl = 0 for pi-, K-
 654 *  |  |  |   ( 1 - Ibl ) * |Icl| = 0 for ap, an, alambdas
 655 *  |  |  |   ( 1 - Ibl ) * |Icl| = 1 for K-, pi-
 656                EMD = AM (IT) * ( ( ( IBAR (ITTA) - IBAR(IT) ) * IBL
 657      &             + ( 1 - IBL ) * ABS (ICL) ) )
 658 *  |  |  |   Finally Emd = 2 Am(it) for antibaryons and Am (it) for
 659 *  |  |  |   pi- and K-. Other particles cannot enter the if, however
 660 *  |  |  |   for baryons it is 0 regardless of the charge and for
 661 *  |  |  |   positive charged mesons is Am (it): the conclusion is
 662 *  |  |  |   that who wrote the code is a ......... . Even not
 663 *  |  |  |   taking into account the "if" all the stuff could be:
 664 *              EMD = AM (IT) * ( 1 - IBAR (IT) )
 665 C
 666 C--------------------------------------------------
 667 C*** EMD=2*AM(IT) FOR ANNIHILATION, =AM(IT) FOR MESONS, =0ELSE
 668 C--------------------------------------------------
 669                THRES = THRESR (IT)
 670 *  |  |  | Ecms is the total available energy in the Cms frame
 671 *  |  |  | and of course is the invariant mass of the projectile-
 672 *  |  |  | nucleon system
 673                ECMS  = SQRT  (ECMO)
 674                IF (ECMS .LT. THRES) EMD = 0.D0
 675 *  |  |  | This card means Emd = 0 always for pi-, K- !! and also
 676 *  |  |  | Emd = 0 for ap, an and alambdas if Hadrin will be called
 677 *  |  |  | (they are not stopping!!)
 678                IF (NXPI .LT. 1) EMD = 0.D0
 679 C
 680 C--------------------------------------------------
 681 C*** EMD=0 FOR MESONS AND BARYONS,IF NO STOPPING OF PARTICLES IF H-N-CMS
 682 C*** -ENERGY LT TABULATED THRESHOLD
 683 C--------------------------------------------------
 684 *  |  |  | This card means Emd = 0 always for pi-, K- (and baryons, but
 685 *  |  |  | for them we don't enter the if)
 686                IF ( IBAR(IT) .GE. 0 ) EMD = 0.D0
 687 C--------------------------------------------------
 688 C*** ECI = TOTAL EFFECTIVE ENERGY IN H-N-CMS-SYSTEM
 689 C--------------------------------------------------
 690 *  |  |  | Eci is the total energy of the projectile in the frame
 691 *  |  |  | were the target is at rest, if the "effective" invariant
 692 *  |  |  | mass id Ecms + Emd (see kinematics inside Ferevv and
 693 *  |  |  | Difevv )
 694                ECI   = 0.5D+00 * ( ( ECMS + EMD )**2 - AIT2 - ATA2 ) /
 695      &                 AMMTA
 696                TOEFF = ECI - AMMIT
 697 *  |  |  |  +----------------------------------------------------------*
 698 *  |  |  |  |
 699                IF ( TOEFF .LT. TLAB ) THEN
 700                   TOEFF = TLAB
 701                END IF
 702 *  |  |  |  |
 703 *  |  |  |  +----------------------------------------------------------*
 704             END IF
 705 *  |  |  |
 706 *  |  |  +-------------------------------------------------------------*
 707             LCORIN = .TRUE.
 708 *  |  | Now the entry Corrnc in Corrin set the correct values for the
 709 *  |  | excitation energy Tv and the cascade energies Tpk and Tnk
 710             CALL CORRNC ( TOEFF )
 711             CALL GRANOR(RANGS1, RANGS2)
 712             TNKTE = MAX ( 0.003D+00, EKMNNU (2)
 713      &            * ( 0.5D+00 - 0.25D+00 * RANGS1 ) )
 714             TPKTE = MAX ( 0.003D+00, EKMNNU (1)
 715      &            * ( 0.5D+00 - 0.25D+00 * RANGS2 ) )
 716             TV  = TVGRE0
 717             TPK = EINCP
 718             TNK = EINCN
 719 C--------------------------------------------------
 720 C*** ELPP=TOTAL PROJECTILE ENERGY,PLPP=ABSOLUT MOMENTUM OF THE PROJECTIL
 721 C*** IN THE LABSYSTEM
 722 C
 723 C--------------------------------------------------
 724             T3   = AM (IT)
 725             TPNVK = TV + TNK + TPK
 726             ELPP  = ELAB - TPNVK
 727 *  |  |  +-------------------------------------------------------------*
 728 *  |  |  | Check for stopping pi- and K-: in my mind this point can
 729 *  |  |  | never be reached by stopping pi-, K- since they should be
 730 *  |  |  | already gone to 405. Anyway...
 731             IF ( INXPI .EQ. 0 .AND. NXPI .NE. 0 ) THEN
 732 *  |  |  |  +----------------------------------------------------------*
 733 *  |  |  |  |    Again in my mind Inxpi = 0 <==> Ibar (itttt) = 0,
 734 *  |  |  |  |    anyway ....
 735                IF ( IBAR (ITTTT) .GE. 0 ) THEN
 736                   TVEUZ = 0.D+00
 737                   V0WELL (ITJ) = V0WELL (ITJ) - V0EXTR
 738                   ANOW  = ANOW  + IBAR (IT)
 739                   KTARN = KTARN - IBAR (IT) + ICH (IT)
 740                   ZNOW  = ZNOW  + ICH  (IT)
 741                   KTARP = KTARP - ICH  (IT)
 742                   GO TO 405
 743 *  |  |  |  |-->-->-->-->--> go to the cascade simulation !!
 744                END IF
 745 *  |  |  |  |
 746 *  |  |  |  +----------------------------------------------------------*
 747             END IF
 748 *  |  |  |
 749 *  |  |  +-------------------------------------------------------------*
 750 C
 751 C--------------------------------------------------
 752 C*** IF NO PRIMARY PI-,K- OR IF HADRIN FOR PI-,K- WILL BE CALLED
 753 C*** (NO STOPPING PARTICLE CASE) GO TO 403
 754 C--------------------------------------------------
 755 C
 756 C--------------------------------------------------
 757 C*** CASE OF ENERGY CORRECTION IF LOOP IN ANNIHILATION
 758 C--------------------------------------------------
 759 *  |  |  +-------------------------------------------------------------*
 760 *  |  |  |
 761             IF ( IBAR (IT) .GE. 0 .AND. NXPI .EQ. 0 ) THEN
 762                ELIM = T3 + TVEPSI
 763 *  |  |  |
 764 *  |  |  +-------------------------------------------------------------*
 765 *  |  |  |
 766             ELSE
 767                ELIM = AM (13)
 768             END IF
 769 *  |  |  |
 770 *  |  |  +-------------------------------------------------------------*
 771  6530       CONTINUE
 772 *  |  |  +-------------------------------------------------------------*
 773 *  |  |  |   If Elpp < 0 for less then 10 iterations (I653), start
 774 *  |  |  |   again with effective collision energy calculation, else
 775 *  |  |  |   correct TPK, TNK and TV. Now the maximum number of itera-
 776 *  |  |  |   tion has been decreased to 4, and we check Elpp < T3.
 777 *  |  |  |   We have no enough energy or we have looped too many
 778 *  |  |  |   times, so give up and force the "effective" energy to
 779 *  |  |  |   be 0, reducing cascade and excitation energy
 780             IF ( I653 .GE. 4 .AND. ELPP .LT. ELIM ) THEN
 781 *  |  |  |   The following cards simply take away the negative
 782 *  |  |  |   "Elpp" to cascade and excitation in the same proportion
 783 *  |  |  |   they were
 784 *  |  |  |   Modified to take into account the Fermi energy:
 785                A65 = 1.D+00 + ( ELPP - ELIM ) / TPNVK
 786                TPK = TPK * A65
 787                TNK = TNK * A65
 788                TV  = TV  * A65
 789                TPNVK  = TPNVK * A65
 790                TVGRE0 = TV
 791                ELPP   = ELIM
 792 *  |  |  |
 793 *  |  |  +-------------------------------------------------------------*
 794 *  |  |  |
 795             ELSE IF ( ELPP .LT. 0.D+00 ) THEN
 796                GO TO 653
 797 *  |  |  |
 798 *  |  +-<|--<--<--<--< go to resampling Tv, Tnk, Tpk
 799             END IF
 800 *  |  |  |
 801 *  |  |  +-------------------------------------------------------------*
 802
 803 *  |  |  +-------------------------------------------------------------*
 804 *  |  |  |              If the available energy is larger than the
 805 *  |  |  |              mass energy of the projectile check the inva-
 806 *  |  |  |              riant mass of the system
 807             IF ( ELPP .GT. T3 .AND. NXPI .EQ. 0 ) THEN
 808                IUMO = 0
 809                DEX  = 0.D+00
 810 *  |  |  |  +----------------------------------------------------------*
 811 *  |  |  |  |
 812 2500           CONTINUE
 813 *  |  |  |  |   Now also here work with atomic masses
 814                   IUMO = IUMO + 1
 815                   PLPPSQ = ( ELPP + AMMIT ) * ( ELPP - AMMIT )
 816                   PLPP   = SQRT ( PLPPSQ )
 817                   PPLUS  = PLAB - PLPP
 818                   EPLUS  = ELAB - ELPP
 819                   EEX    = TVEUZ  + DEX
 820                   AMSTAR = AMMRES + EEX
 821                   EKRES  = SQRT ( AMSTAR**2 + PPLUS**2 + PFRMSQ - 2.D+00
 822      &                   * PFRSCA * PPLUS ) - AMSTAR
 823 *  |  |  |  |   Note +DEX to reduce actually the excitation energy, not
 824 *  |  |  |  |   the recoil one (it would have been more exact to set
 825 *  |  |  |  |   Efrm = Efrm - Dex, but it is better to conserve Efrm
 826 *  |  |  |  |   coherent with the original value)
 827                   EKRECL = EKRES - EPLUS + DEX
 828 *  |  |  |  |  +-------------------------------------------------------*
 829 *  |  |  |  |  |  Check the recoil energy:
 830                   IF ( EKRECL .LT. 0.D+00 .AND. IUMO .LE. 2 ) THEN
 831                      EKRECL = EKREC0
 832                      LDELTX = .FALSE.
 833 *  |  |  |  |  |
 834 *  |  |  |  |  +-------------------------------------------------------*
 835 *  |  |  |  |  |
 836                   ELSE
 837                      LDELTX = .TRUE.
 838                   END IF
 839 *  |  |  |  |  |
 840 *  |  |  |  |  +-------------------------------------------------------*
 841                   EZERO  = EFRM  - V0WELL (ITJ) - EKRECL - EBNDNG (ITJ)
 842                   ESPENT = ELPP  + EZERO
 843                   ELEFT  = ETTOT - ESPENT
 844                   EELEFT = ELEFT + DELEFT
 845                   UMO2   = AIT2 + EZERO**2 - PFRMSQ + 2.D+00 * ( EZERO
 846      &                   * ELPP - PLPP * PFRSCA )
 847                   DELTU2 = UMIN2 - UMO2
 848 *  |  |  |  |  +-------------------------------------------------------*
 849 *  |  |  |  |  |  We need more invariant mass for the collision!!!
 850                   IF ( DELTU2 .GT. 0.D+00 ) THEN
 851 *  |  |  |  |  |  +----------------------------------------------------*
 852 *  |  |  |  |  |  |
 853                      IF ( IUMO .GT. 4 ) THEN
 854                         GO TO 653
 855 *  |  |  |  |  |  |
 856 *  |  +-<|-<|-<|-<|--<--<--<
 857                      END IF
 858 *  |  |  |  |  |  |
 859 *  |  |  |  |  |  +----------------------------------------------------*
 860
 861 *  |  |  |  |  |  +----------------------------------------------------*
 862 *  |  |  |  |  |  |
 863                      IF ( LDELTX .AND. TPNVK .GT. 0.D+00 ) THEN
 864                         AHELP  = ESPENT / ELEFT
 865                         FFRAC  = ( EEX - TVEPSI ) / TPNVK
 866                         BHELP  = - ELPP - ELPP / PLPP * ( PFRSCA + AHELP
 867      &                         * ( PFRSCA - PPLUS ) )
 868                         DELTAE = 0.5D+00 * DELTU2 / ( BHELP + FFRAC *
 869      &                           AMSTAR * AHELP )
 870                         TMPDEL = 1.01D+00 * DELTAE
 871                         DELTAE = MIN ( TMPDEL, TPNVK )
 872                         DELTEX = FFRAC * DELTAE
 873                         TKOR = 1.D+00 - DELTAE / TPNVK
 874                         TPK  = TPK * TKOR
 875                         TNK  = TNK * TKOR
 876                         TV   = TV  * TKOR
 877                         TVGRE0 = TV
 878                         TPNVK  = TPNVK - DELTAE
 879 *  |  |  |  |  |  |
 880 *  |  |  |  |  |  +----------------------------------------------------*
 881 *  |  |  |  |  |  |
 882                      ELSE IF ( TPNVK .GT. 0.D+00 ) THEN
 883                         DELTAE = 0.5D+00 * DELTU2 / ( EZERO - PFRSCA *
 884      &                           ELPP / PLPP )
 885                         TMPDEL = 1.01D+00 * DELTAE
 886                         DELTAE = MIN ( TMPDEL, TPNVK )
 887                         DELTEX = 0.D+00
 888                         TKOR = 1.D+00 - DELTAE / TPNVK
 889                         TPK  = TPK * TKOR
 890                         TNK  = TNK * TKOR
 891                         TV   = TV  * TKOR
 892                         TVGRE0 = TV
 893                         TPNVK  = TPNVK - DELTAE
 894 *  |  |  |  |  |  |
 895 *  |  |  |  |  |  +----------------------------------------------------*
 896 *  |  |  |  |  |  |
 897                      ELSE IF ( LDELTX ) THEN
 898                         DELTAE = 0.D+00
 899                         DELTEX = 0.5D+00 * DELTU2 * ELEFT / ( ESPENT *
 900      &                           AMSTAR )
 901                         TMPDEL = 1.01D+00 * DELTEX
 902                         DELTEX = MIN ( TMPDEL, EEX - TVEPSI )
 903 *  |  |  |  |  |  |
 904 *  |  |  |  |  |  +----------------------------------------------------*
 905 *  |  |  |  |  |  |
 906                      ELSE
 907                         DELTAE = 0.D+00
 908                         DELTEX = 0.D+00
 909                      END IF
 910 *  |  |  |  |  |  |
 911 *  |  |  |  |  |  +----------------------------------------------------*
 912 *                    ELPP   = ELPP  + DELTAE, AMMIT
 913                      TMPAMM = AMMIT + 0.01D+00
 914                      ELPP   = MAX ( ELPP  + DELTAE, TMPAMM )
 915                      DEX    = DEX   - DELTEX
 916                      GO TO 2500
 917 *  |  |  |  |  |
 918 *  |  |  |  +-<|--<--<--<--<--< go to compute the new invariant mass!!
 919 *  |  |  |  |  |
 920 *  |  |  |  |  +-------------------------------------------------------*
 921 *  |  |  |  |  |
 922                   ELSE
 923                      TVEUZ = EEX
 924                      GO TO 510
 925 *  |  |  |  |  |
 926 *  |  |  |  |  |-->-->-->-->--> invariant mass larger than the projec-
 927 *  |  |  |  |  |                tile mass plus target nucleon mass, go
 928 *  |  |  |  |  |                to the event simulation !!
 929                   END IF
 930 *  |  |  |  |  |
 931 *  |  |  |  |  +-------------------------------------------------------*
 932 *  |  |  |  |
 933 *  |  |  |  +----------------------------------------------------------*
 934             END IF
 935 *  |  |  |
 936 *  |  |  +-------------------------------------------------------------*
 937             IF ( IBAR (IT) .GE. 0 ) GO TO 653
 938 *  |  |
 939 *  |  +--<--<--<--<--< go to resampling Tv, Tnk, Tpk, Elpp < T3 and no
 940 *  |  |                antinucleon projectile
 941             PLPPSQ = ( ELPP - AMMIT ) * ( ELPP + AMMIT )
 942 *  |  |  +-------------------------------------------------------------*
 943 *  |  |  |
 944             IF ( PLPPSQ .GT. 0.D+00 ) THEN
 945                PLPP = SQRT ( PLPPSQ )
 946 *  |  |  |
 947 *  |  |  +-------------------------------------------------------------*
 948 *  |  |  |
 949             ELSE
 950                PLPP = 0.D+00
 951             END IF
 952 *  |  |  |
 953 *  |  |  +-------------------------------------------------------------*
 954             PPLUS  = PLAB - PLPP
 955             EPLUS  = ELAB - ELPP
 956             EEX    = TVEUZ
 957             AMSTAR = AMNRES + EEX
 958             EKRES  = SQRT ( AMSTAR**2 + PPLUS**2 + PFRMSQ - 2.D+00
 959      &             * PFRSCA * PPLUS ) - AMSTAR
 960             EKRECL = EKRES - EPLUS
 961 *  |  |  +-------------------------------------------------------------*
 962 *  |  |  |  Check the recoil energy:
 963             IF ( EKRECL .LT. EKREC0 ) THEN
 964                EKRECL = EKREC0
 965             END IF
 966 *  |  |  |
 967 *  |  |  +-------------------------------------------------------------*
 968             EZERO  = EFRM  - V0WELL (ITJ) - EKRECL - EBNDNG (ITJ)
 969             ESPENT = ELPP  + EZERO
 970             ELEFT  = ETTOT - ESPENT
 971             UMO2   = AIT2 + EZERO**2 - PFRMSQ + 2.D+00 * ( EZERO
 972      &             * ELPP - PLPP * PFRSCA )
 973             DELTU2 = UMIN2 - UMO2
 974             IF ( DELTU2 .LT. 0.D+00 ) GO TO 510
 975 *  |  |
 976 *  |  |-->-->-->-->--> invariant mass larger than the projec-
 977 *  |  |                tile mass plus target nucleon mass, go
 978 *  |  |                to the event simulation !!
 979 *  |  |  If we are here we are dealing with an antibaryon and the total
 980 *  |  |  available energy in the centre of mass frame is less than the
 981 *  |  |  projectile mass plus the target nucelon mass
 982 C--------------------------------------------------
 983 C*** IF ELPP GREATERTHEN THE MASS OF IT, GO TO MOMENTUM LIMITATION FOR H
 984 C*** ELSE IF NO EFFECTIVE  ENERGY WAS CALCULATED DO THE SAME
 985 C*** ELSE IF IT IS NO ANTINUCLEON AND ELPP LT MASS OF IT, START AGAIN WI
 986 C***   EFFECTIVE COLLISION ENERGY CALCULATION
 987 C*** ELSE NEW ENERGY DEFINITION FOR IT AND NEW MASS DEFINITION FOR BARYO
 988 C*** AND ANTIBARYONS
 989 C--------------------------------------------------
 990 *  |  |  +-------------------------------------------------------------*
 991 *  |  |  |
 992             IF ( EMD .LE. 0.D+00 ) THEN
 993                ELPP = T3
 994                PLPP = 0.D+00
 995                GO TO 1510
 996 *  |  |  |
 997 *  |  |  |-->-->-->-->--> Skip mass redefinition if Emd is negative or
 998 *  |  |  | zero. In my opinion now Emd is always = 0 for non stopping
 999 *  |  |  | antibaryons, because of the if inside the Emd stuff and I am
1000 *  |  |  | not sure we have really to arrive here under the minimum
1001 *  |  |  | invariant mass...
1002             END IF
1003 *  |  |  |
1004 *  |  |  +-------------------------------------------------------------*
1005 C
1006 C--------------------------------------------------
1007 C*** CASE OF ANNIHILATION,PSEUDO MASS DEFINITION
1008 C--------------------------------------------------
1009             PLPP = 0.D+00
1010             UMO2 = ( ELPP + EZERO )**2 - PFRMSQ
1011             T3   = 0.5D+00 * SQRT ( UMO2 )
1012             IUMO = 0
1013             INEG = 0
1014 *  |  |  +-------------------------------------------------------------*
1015 *  |  |  |
1016 5050        CONTINUE
1017                T3   = 0.9D+00 * T3
1018                T3SQ = T3 * T3
1019                PLPPSQ = ELPP**2 - T3SQ
1020 *  |  |  |  +----------------------------------------------------------*
1021 *  |  |  |  |
1022                IF ( PLPPSQ .GT. 0.D+00 ) THEN
1023                   IUMO = IUMO + 1
1024                   PLPP = SQRT ( PLPPSQ )
1025 *  |  |  |  |
1026 *  |  |  |  +----------------------------------------------------------*
1027 *  |  |  |  |
1028                ELSE
1029                   INEG = INEG + 1
1030 *  |  |  |  |  +-------------------------------------------------------*
1031 *  |  |  |  |  |
1032                   IF ( INEG .GT. 100 ) THEN
1033                      WRITE (LUNOUT,*)
1034      &               ' Nucriv: impossible to get the pseudo-mass!!',
1035      &               ' Plpp always negative, Plpp,Elpp', PLPP, ELPP
1036                      WRITE (LUNERR,*)
1037      &               ' Nucriv: impossible to get the pseudo-mass!!',
1038      &               ' Plpp always negative, Plpp,Elpp', PLPP, ELPP
1039                      V0WELL (ITJ) = V0WELL (ITJ) - V0EXTR
1040                      ANOW  = ANOW  + IBAR (IT)
1041                      KTARN = KTARN - IBAR (IT) + ICH (IT)
1042                      ZNOW  = ZNOW  + ICH  (IT)
1043                      KTARP = KTARP - ICH  (IT)
1044                      GO TO 405
1045                   END IF
1046 *  |  |  |  |  |
1047 *  |  |  |  |  +-------------------------------------------------------*
1048                   GO TO 5050
1049                END IF
1050 *  |  |  |  |
1051 *  |  |  |  +----------------------------------------------------------*
1052                PPLUS  = PLAB - PLPP
1053                EPLUS  = ELAB - ELPP
1054                EEX    = TVEUZ
1055                AMSTAR = AMNRES + EEX
1056                EKRES  = SQRT ( AMSTAR**2 + PPLUS**2 + PFRMSQ - 2.D+00
1057      &                * PFRSCA * PPLUS ) - AMSTAR
1058                EKRECL = EKRES - EPLUS
1059 *  |  |  |  +----------------------------------------------------------*
1060 *  |  |  |  |  Check the recoil energy:
1061                IF ( EKRECL .LT. EKREC0 ) THEN
1062                   EKRECL = EKREC0
1063                END IF
1064 *  |  |  |  |
1065 *  |  |  |  +----------------------------------------------------------*
1066                EZERO  = EFRM - V0WELL (ITJ) - EKRECL - EBNDNG (ITJ)
1067                UMO2   = T3SQ + EZERO**2 - PFRMSQ + 2.D+00 * ( EZERO
1068      &                * ELPP - PLPP * PFRSCA )
1069                DELTU2 = 2.D+00 * T3 - UMO2
1070 *  |  |  |  +----------------------------------------------------------*
1071 *  |  |  |  |
1072                IF ( IUMO .GT. 10 ) THEN
1073                   WRITE (LUNERR,*)
1074      &            ' Nucriv: impossible to get the pseudo-mass!!',
1075      &              IUMO,DELTU2
1076                   V0WELL (ITJ) = V0WELL (ITJ) - V0EXTR
1077                   ANOW  = ANOW  + IBAR (IT)
1078                   KTARN = KTARN - IBAR (IT) + ICH (IT)
1079                   ZNOW  = ZNOW  + ICH  (IT)
1080                   KTARP = KTARP - ICH  (IT)
1081                   GO TO 405
1082                END IF
1083 *  |  |  |  |
1084 *  |  |  |  +----------------------------------------------------------*
1085                IF ( DELTU2 .GT. 0.D+00 ) GO TO 5050
1086 *  |  |  |
1087 *  |  |  |--<--<--<--<--< invariant mass too small !!
1088 *  |  |  +-------------------------------------------------------------*
1089 *  |  |  +-------------------------------------------------------------*
1090 *  |  |  |   Check the pseudo-mass value (note that the following check
1091 *  |  |  |   becomes useless (except for I653 > 10)
1092             IF ( T3 .LT. 0.3D+00 * AM (1) .AND. I653 .LT. 10 ) THEN
1093                I653 = MAX ( I653 + 1, 4 )
1094                ELIM = 0.33D+00 * AM (1) * ELPP / T3
1095                T3 = AM (IT)
1096                GO TO 6530
1097 *  |  |  |--<--<--<--<--< pseudo-mass is too small !!
1098             END IF
1099 *  |  |  |
1100 *  |  |  +-------------------------------------------------------------*
1101 *  |  |  +-------------------------------------------------------------*
1102 *  |  |  | However let a warning message if it is met!!!
1103             IF ( T3 .LT. 0.14D+00 .AND. INXPI .LT. 0 .AND. NXPI .GT. 0 )
1104      &           THEN
1105                NXPI  = 1
1106                V0WELL (ITJ) = V0WELL (ITJ) - V0EXTR
1107                ANOW  = ANOW  + IBAR (IT)
1108                KTARN = KTARN - IBAR (IT) + ICH (IT)
1109                ZNOW  = ZNOW  + ICH  (IT)
1110                KTARP = KTARP - ICH  (IT)
1111                GO TO 405
1112 *  |  |  |
1113 *  |  |  |-->-->-->-->--> go to the cascade simulation
1114             END IF
1115 *  |  |  |
1116 *  |  |  +-------------------------------------------------------------*
1117
1118 *  |  |  +-------------------------------------------------------------*
1119 *  |  |  |  Enter this if only with anti-nucleons!!!!!
1120             IF ( IBAR (ITTTT) .LT. 0 ) THEN
1121                ETDD  = 0.D0
1122 *  |  |  |  Etkor
1123                ETKOR = TPNVK + AMMIT
1124 1800           CONTINUE
1125 *  |  |  |  +----------------------------------------------------------*
1126 *  |  |  |  | Check if it was a stopping particle or not !!
1127                IF ( NXPI .LE. 0 ) THEN
1128                   ET = ELAB
1129 C
1130 C--------------------------------------------------
1131 C*** THE PARTICLE IT (=AN) IS NOT STOPPED IN THE NUCLEUS
1132 C--------------------------------------------------
1133                   ETREST = ET - ETKOR
1134                   IF ( ETREST .LT. 0.D0 ) ETDD = ETREST
1135                   TKOR = 1.D0 + ETDD / ( TPNVK + 1.D-10 )
1136 *  |  |  |  |  +-------------------------------------------------------*
1137 *  |  |  |  |  |   Treat it as a stopped particle
1138                   IF ( TKOR .LT. 1.D-6 ) THEN
1139                      NXPI = 1
1140                      GO TO 1800
1141 *  |  |  |  |  |
1142 *  |  |  |  |  +-------------------------------------------------------*
1143 *  |  |  |  |  |        update tv, tnk, tpk to account for etdd < 0
1144                   ELSE
1145                      TV  = TV  * TKOR
1146                      TPK = TPK * TKOR
1147                      TNK = TNK * TKOR
1148                      TVGRE0 = TV
1149                      TPNVK = TPK + TNK + TV
1150                      ELPP= ET - TPNVK
1151                      GO TO 510
1152                   END IF
1153 *  |  |  |  |  |
1154 *  |  |  |  |  +-------------------------------------------------------*
1155 *  |  |  |  |
1156 *  |  |  |  +----------------------------------------------------------*
1157 *  |  |  |  |       It is a stopping antibaryon !!
1158                ELSE
1159                   ETDD  = 0.D0
1160 C
1161 C--------------------------------------------------
1162 C*** THE PARTICLE IT (=AN) IS STOPPED IN THE NUCLEUS
1163 C--------------------------------------------------
1164                   ET = ELAB + AMMTA
1165                   ETREST = ET - ETKOR
1166                   IF ( ETREST .LT. 0.D0 ) ETDD = ETREST
1167                   TKOR = 1.D+00 + ETDD / TPNVK
1168                   IF ( TKOR .LT. 1.D-6 ) GO TO 653
1169 *  |  |  |  |
1170 *  |  +-<|-<|--<--<--< go to resampling Tv, Tnk, Tpk
1171                   TV  = TV  * TKOR
1172                   TPK = TPK * TKOR
1173                   TNK = TNK * TKOR
1174                   TVGRE0 = TV
1175                   TPNVK = TPK + TNK + TV
1176                   ELPP = ELAB - TPNVK
1177                   PLPP = SQRT ( ELPP**2 - T3SQ )
1178                END IF
1179 *  |  |  |  |
1180 *  |  |  |  +----------------------------------------------------------*
1181             END IF
1182 *  |  |  |
1183 *  |  |  +-------------------------------------------------------------*
1184 *  |  | end of the "653" loop for bad energy correction in annihilation
1185 *  |  +----------------------------------------------------------------*
1186 C
1187 C--------------------------------------------------
1188 C*** NEW MASS DEFINITION OF BARYONS AND ANTIBARYONS
1189 C*** (ONLY FOR REMAINING EFFECTIVE CM-ENERGIES,LOWER THAN THE TOTAL
1190 C*** 1.88GEV-THRESHOLD CM-ENERGY IN ANNIHILATION)
1191 C--------------------------------------------------
1192 *  |  +----------------------------------------------------------------*
1193 *  |  |  First pass through mass redefinition
1194          IF ( .NOT. LVMASS ) THEN
1195             LVMASS = .TRUE.
1196 *  |  |  +-------------------------------------------------------------*
1197 *  |  |  |  Loop over masses to be redefined
1198             DO 1511 IAM=1,15
1199                JAM        = IAMC (IAM)
1200                AMHH (IAM) = AM   (JAM)
1201                IF ( IBAR (JAM) .NE. 0 ) AM (JAM) = MIN  ( T3, AM (JAM) )
1202 1511        CONTINUE
1203 *  |  |  |
1204 *  |  |  +-------------------------------------------------------------*
1205 *  |  |
1206 *  |  +----------------------------------------------------------------*
1207 *  |  |  Masses already redefined
1208          ELSE
1209 *  |  |  +-------------------------------------------------------------*
1210 *  |  |  |  Loop over masses to be redefined
1211             DO 15111 IAM=1,15
1212                JAM        = IAMC (IAM)
1213                IF ( IBAR (JAM) .NE. 0 ) AM (JAM) = MIN ( T3, AMHH (IAM))
1214 15111       CONTINUE
1215 *  |  |  |
1216 *  |  |  +-------------------------------------------------------------*
1217          END IF
1218 *  |  |
1219 *  |  +----------------------------------------------------------------*
1220 *  | We will arrive here also with T3 = AM (IT) and ELPP < T3 from
1221 *  | the if on EMD =< 0, it is not clear if it is correct!!
1222 1510     CONTINUE
1223 *  |  +----------------------------------------------------------------*
1224 *  |  |                             Loop for Ferhad/hadrin call
1225 510      CONTINUE
1226 C
1227 C--------------------------------------------------
1228 C*** MOMENTUM LIMITATION FOR HADRIN IN FERMI-MOM.-VERSION
1229 C--------------------------------------------------
1230             PLABOU=15.D0
1231 *  |  |  +-------------------------------------------------------------*
1232 *  |  |  |  This check seems to me to be useless, since the maximum
1233 *  |  |  |  allowed momentum is 10 GeV/c .... anyway get a warning
1234 *  |  |  |  message if this condition is met
1235             IF ( PLPP .GE. PLABOU ) THEN
1236                WRITE ( LUNERR, * )' Nucriv: momentum limitation met!!',
1237      &                              PLPP, PLABOU
1238                APL  = PLPP - PLABOU
1239                ELPP = SQRT (PLABOU**2 + AM(IT)**2)
1240                TPK  = TPK + APL*.3D0
1241                TNK  = TNK + APL*.3D0
1242                TV   = TV  + APL*.4D0
1243                IFORBI = IFORBI + 1
1244                IF ( IFORBI .GT. 10 ) GO TO 405
1245             END IF
1246 *  |  |  |
1247 *  |  |  +-------------------------------------------------------------*
1248             EELAB = ELPP
1249             PPLAB = PLPP
1250             IHACA = IHACA + 1
1251             ITFRHD = IT
1252             IF ( PPLAB .LT. 1.D-04 ) THEN
1253                WRITE (LUNERR,*)
1254      &              ' Nucriv: iperco=0,eelab,pplab,it',
1255      &                EELAB,PPLAB,IT
1256                WRITE (LUNERR,*)'      IHACA,LVMASS,ELAB,PLAB',
1257      &                                IHACA,LVMASS,ELAB,PLAB
1258             END IF
1259             CALL FERHAV ( IT, EELAB, PPLAB, CX, CY, CZ )
1260 C
1261 C--------------------------------------------------
1262 C***
1263 C***
1264 C*** HADRON NUCLEON COLLISION SIMULATION
1265 C*** EELAB=LABENERGY OF H,PPLAB=ABSOLUT MOMENTUM OF H, CX,CY,CZ DIRECTIO
1266 C*** COSINES OF H,IN THE LAB SYSTEM, ITTA TARGET NUCLEON INDEX,
1267 C*** ANUC,ZNUC NUCLEON AND PROTON NUMBERS
1268 C***
1269 C*** REGARDING FERMI MOMENTUM OF THE NUCLEON
1270 C***
1271 C--------------------------------------------------
1272 C--------------------------------------------------
1273 C***
1274 C***
1275 C*** IHACA TIMES CALLED H-N-COLL., IF IN FINAL STATE NO PARTICLE,CALL
1276 C*** AGAIN,
1277 C*** IF IN ANNIHILATION OF A STOPPING NUCLEON,   IN THE FINAL
1278 C*** STATE IS A BARYON, CALL AGAIN H-N-COLL.
1279 C*** (NOT MORE THAN 30 TIMES, ELSE TAKE THIS STATE)
1280 C***
1281 C***
1282 C--------------------------------------------------
1283 *  |  |  +-------------------------------------------------------------*
1284 *  |  |  | Check if we have reached the maximum numbers of calls (which
1285 *  |  |  | has been reduced to 15 .... 30 seems to be a bit too large
1286 *  |  |  | but there were two "IHACA = IHACA + 1" ....
1287             IF ( IHACA .GT. 15 ) THEN
1288               GO TO 88
1289 *  |  |  |-->-->-->-->--> give up !!!
1290             END IF
1291 *  |  |  |
1292 *  |  |  +-------------------------------------------------------------*
1293 *  |  | Ir is the number of secondaries produced in the Ferhad/Hadrin
1294 *  |  | call
1295 *  |  |  +-------------------------------------------------------------*
1296 *  |  |  |
1297             IF (IR .EQ. 0) THEN
1298                V0WELL (ITJ) = V0WELL (ITJ) + V0EXTR
1299                GO TO 510
1300 *  |  |  |
1301 *  +-<|-<|-<--<--<--<--< go to resampling the ferhad/hadrin interaction
1302             END IF
1303 *  |  |  |
1304 *  |  |  +-------------------------------------------------------------*
1305 *  |  | Check that for annihilated antibaryons we have no (anti)baryon
1306 *  |  | in the final state: if yes resample !
1307 *  |  |  +-------------------------------------------------------------*
1308 *  |  |  |
1309             IF ( INXPI .LT. 0 .AND. NXPI .GT. 0 ) THEN
1310                IBAROT = 0
1311                DO 520 IRZ = 1, IR
1312                   IBAROT = IBAROT + ABS (IBAR(ITR(IRZ)))
1313   520          CONTINUE
1314                IF ( IBAROT .GT. 0 ) THEN
1315                   V0WELL (ITJ) = V0WELL (ITJ) + V0EXTR
1316                   GO TO 510
1317                END IF
1318 *  |  |  |
1319 *  +-<|-<|-<--<--<--<--< go to resampling the ferhad/hadrin interaction
1320             END IF
1321 *  |  |  |
1322 *  |  |  +-------------------------------------------------------------*
1323 C
1324 C--------------------------------------------------
1325 C*** IF NO ANNIHILATION INTO MESONS FOR AN+N
1326 C--------------------------------------------------
1327 *  |  | Of course ir=1 means that the only outgoing particle is still
1328 *  |  | the projectile !! If we have an incoming antibaryon
1329 *  |  | we need at least 2 particles in the final state (why????, if
1330 *  |  | we have no annihilation for example a charge exchange recation
1331 *  |  | can give ir=2)
1332 *  |  |  +-------------------------------------------------------------*
1333 *  |  |  |
1334             IF ( IR .LT. 2 .AND. INXPI .LT. 0 ) THEN
1335                V0WELL (ITJ) = V0WELL (ITJ) + V0EXTR
1336                GO TO 510
1337 *  |  |  |
1338 *  +-<|-<|-<--<--<--<--< go to resampling the ferhad/hadrin interaction
1339             END IF
1340 *  |  |  |
1341 *  |  |  +-------------------------------------------------------------*
1342 88       CONTINUE
1343 *  |  |                                End loop for Ferhad/hadrin call
1344 *  |  +----------------------------------------------------------------*
1345
1346 *  |  +----------------------------------------------------------------*
1347 *  |  |  Redefinition of masses, if they were changed for annihilation
1348 *  |  |  beyond 1.88 GeV - threshold
1349          IF ( LVMASS ) THEN
1350 *  |  |  +-------------------------------------------------------------*
1351 *  |  |  |
1352             DO 1512 IAM=1,15
1353                AM (IAMC(IAM)) = AMHH (IAM)
1354 1512        CONTINUE
1355 *  |  |  |
1356 *  |  |  +-------------------------------------------------------------*
1357          END IF
1358 *  |  |
1359 *  |  +----------------------------------------------------------------*
1360 C
1361 C--------------------------------------------------
1362 C*** CORRECTION OF CASCADE AND EXCITATION ENERGIES FOR FERMI-MOMENTUM
1363 C*** IF ENERGY CONSERVATION ALLOWS NO H-N-COLLISION, GO TO THE START
1364 C*** POINT OF EVENT SIMUL.
1365 C*** ELSE ENERGY CORRECTION BY VARIABLE EKIKOR
1366 C--------------------------------------------------
1367 C--------------------------------------------------
1368 2800     CONTINUE
1369          ECHCK  = ETTOT
1370          PXCHCK = PXTTOT
1371          PYCHCK = PYTTOT
1372          PZCHCK = PZTTOT
1373 *  |  +----------------------------------------------------------------*
1374 *  |  |
1375          DO 3000 IPART = 1, IR
1376             ECHCK  = ECHCK  - EL (IPART)
1377             PXCHCK = PXCHCK - PL (IPART) * CXR (IPART)
1378             PYCHCK = PYCHCK - PL (IPART) * CYR (IPART)
1379             PZCHCK = PZCHCK - PL (IPART) * CZR (IPART)
1380 3000     CONTINUE
1381 *  |  |
1382 *  |  +----------------------------------------------------------------*
1383          UMO2 = (ECHCK + DELEFT)**2 - PXCHCK**2 - PYCHCK**2 - PZCHCK**2
1384 *  |  +----------------------------------------------------------------*
1385 *  |  |    Now moved to atomic masses !!!
1386          IF ( IR .GT. 1 ) THEN
1387             AMMRE2 = AMMRES**2
1388 *  |  |
1389 *  |  +----------------------------------------------------------------*
1390 *  |  |    It is supposed that the only emitted particle is still the
1391 *  |  |    projectile
1392          ELSE
1393             AMMRES = AMMTAR
1394             AMNRES = AMNTAR
1395             AMMRE2 = AMMTAR**2
1396          END IF
1397 *  |  |
1398 *  |  +----------------------------------------------------------------*
1399 *  |  +----------------------------------------------------------------*
1400 *  |  |
1401          IF ( UMO2 .LT. AMMRE2 ) THEN
1402             I1651 = I1651 + 1
1403             IF ( I1651 .GT. 4 .OR. IHACA .GT. 10 ) THEN
1404                LRESMP = .TRUE.
1405                RETURN
1406             END IF
1407             V0WELL (ITJ) = V0WELL (ITJ) + V0EXTR
1408             GO TO 1651
1409 *  |  |
1410 *  +-<|-<--<--<--<--<--< go to the beginning of the event simulation!
1411          END IF
1412 *  |  |
1413 *  |  +----------------------------------------------------------------*
1414 *  |              End of the event simulation: finally we got it
1415 *  +-------------------------------------------------------------------*
1416 *   Update Tv, Tnk, Tpk taking into account Ekikor: simply ekikor is
1417 *   added (with its sign) leaving unchanged the ratios
1418 *   Again this stuff is now obsolete, we make the same from the energy
1419 *   balance
1420 *  +-------------------------------------------------------------------*
1421 *  |
1422       IF ( IPERCO .LT. 0 ) THEN
1423          ELRES  = ECHCK - AMNRES
1424          ELRES0 = TVEUZ + TPNVK
1425          AA65 = ELRES / ELRES0
1426 *  | Tpnvk records the sum of the energy available for cascade protons,
1427 *  | neutrons and excitation
1428 *  |  +----------------------------------------------------------------*
1429 *  |  |
1430          IF ( TPNVK .LE. 0.D+00 ) THEN
1431             TPK = 0.4D+00 * ELRES
1432             TNK = 0.5D+00 * ELRES
1433             TV  = 0.1D+00 * ELRES
1434             TVGRE0 = 0.5D+00 * TV
1435             TVEUZ  = 0.5D+00 * TV
1436 *  |  |
1437 *  |  +----------------------------------------------------------------*
1438 *  |  |
1439          ELSE
1440             TPK = TPK * AA65
1441             TNK = TNK * AA65
1442             TVGRE0 = TVGRE0 * AA65
1443             TVEUZ  = TVEUZ  * AA65
1444             TV = TVGRE0 + TVEUZ
1445          END IF
1446 *  |  |
1447 *  |  +----------------------------------------------------------------*
1448 *  |
1449 *  +-------------------------------------------------------------------*
1450 *  |                           peripheral collision, no cascade
1451       ELSE
1452          ELRES0 = TVEUZ
1453          ELRES  = ECHCK - AMNRES
1454          AA65 = ELRES / ELRES0
1455          TPK = 0.D+00
1456          TNK = 0.D+00
1457          TVGRE0 = 0.D+00
1458          TVEUZ  = TVEUZ * AA65
1459          TV = TVEUZ
1460       END IF
1461 *  |
1462 *  +-------------------------------------------------------------------*
1463 C
1464 C--------------------------------------------------
1465 C
1466 C*** REDEFINE THE ENERGIES OF SAMPLED PARTICLES IN
1467 C*** CASE OF CHANGED MASSES
1468 C
1469 C--------------------------------------------------
1470 *  +-------------------------------------------------------------------*
1471 *  |
1472       IF ( LVMASS ) THEN
1473          MESON=0
1474 *  |  +----------------------------------------------------------------*
1475 *  |  |
1476          DO 2510 IRZ=1,IR
1477             ITRZ   = MIN (ITR(IRZ),30)
1478             IBAROT = JAMC(ITRZ)
1479 *  |  |  +-------------------------------------------------------------*
1480 *  |  |  |   Check if we need to redefine the energy
1481             IF ( ABS (IBAROT) .GE. 1 ) THEN
1482 *  |  |  |   This energy redefinition can be very dangerous for any
1483 *  |  |  |   energy balance, since we retain momentum and not energy!!
1484                ECHCK = ECHCK + EL (IRZ)
1485                EL (IRZ) = SQRT ( PL (IRZ)**2 + AM (ITRZ)**2 )
1486                ECHCK = ECHCK - EL (IRZ)
1487 *  |  |  |
1488 *  |  |  +-------------------------------------------------------------*
1489 *  |  |  |
1490             ELSE
1491                 MESON = MESON + 1
1492             END IF
1493 *  |  |  |
1494 *  |  |  +-------------------------------------------------------------*
1495 2510     CONTINUE
1496 *  |  |
1497 *  |  +----------------------------------------------------------------*
1498       END IF
1499 *  |
1500 *  +-------------------------------------------------------------------*
1501
1502 *  +-------------------------------------------------------------------*
1503 *  |                           No particle in the secondary stack
1504       IF ( IR .EQ. 0 ) THEN
1505          ANOW  = ANOW  + IBAR (IT)
1506          KTARN = KTARN - IBAR (IT) + ICH (IT)
1507          ZNOW  = ZNOW  + ICH  (IT)
1508          KTARP = KTARP - ICH  (IT)
1509          GO TO 405
1510 *  |
1511 *  +-->-->-->-->--> go to the cascade simulation
1512 *  |
1513 *  +-------------------------------------------------------------------*
1514 *  |  Only one particle in the final state, it should be the projectile
1515 *  |  since it is assumed that no nucleon number correction occurs
1516 *  |  ( see below ): give random angle like for cascade particles to the
1517 *  |  single secondary from Hadrin
1518       ELSE IF ( IR .EQ. 1 ) THEN
1519          AMMRES = AMMTAR
1520          AMNRES = AMNTAR
1521          TKI = EL(1) - AM (ITR(1))
1522          IF(TKI .LE. 1.D-20 ) GO TO 1221
1523          ADE=0.12D0*(1.D0+0.003D0*BBTAR)/TKI
1524          DEX=EXP(-1.57D0*1.57D0/ADE)
1525          AN1=(1.D0-DEX)*ADE*.5D0
1526          AN2=DEX*1.57D0
1527          AN=AN1+AN2
1528          AN1=AN1/AN
1529          CALL GRNDM(RNDM,1)
1530          IF (RNDM(1) .GT. AN1) GO TO 1222
1531  1223    CONTINUE
1532 C
1533 C--------------------------------------------------
1534 C*** TETA ANGLE DETERMINATION (SINGLE PART.CASE), = DE
1535 C--------------------------------------------------
1536          CALL GRNDM(RNDM,1)
1537          DE=SQRT(-ADE*LOG(1.D0-RNDM(1)*(1.D0-DEX)))
1538          IF(DE.GT.1.57D0) GO TO 1223
1539          GOTO 1224
1540  1222    CONTINUE
1541          CALL GRNDM(RNDM,1)
1542          DE=-RNDM(1)
1543          DE=ATAN2(SQRT(ONEONE-DE**2),DE)
1544  1224    CONTINUE
1545          CALL COSI(SFE,CFE)
1546 C
1547 C--------------------------------------------------
1548 C*** COS PHI, SIN PHI DETERMINATION(S.P.CASE)
1549 C--------------------------------------------------
1550          SID=SIN(DE)
1551          COD=COS(DE)
1552          CALL TTRANS (CXR(1), CYR(1), CZR(1), COD, SID, SFE, CFE,
1553      &                CCXR, CCYR, CCZR)
1554 C
1555 C--------------------------------------------------
1556 C*** TURNING OF ANGLE S INTO THE LABSYSTEM (S.P.CASE)
1557 C--------------------------------------------------
1558          CXR(1)=CCXR
1559          CYR(1)=CCYR
1560          CZR(1)=CCZR
1561  1221    CONTINUE
1562 *  |
1563 *  +-------------------------------------------------------------------+
1564 *  |                           Two or more particles produced
1565       ELSE
1566          I1 = ITJ - 1
1567          A1 = I1
1568          ANOW = ANOW - 1.D0
1569          ZNOW = ZNOW - 1.D0 + A1
1570          KTARP = KTARP + 1 - I1
1571          KTARN = KTARN + I1
1572       END IF
1573 *  |
1574 *  +-------------------------------------------------------------------+
1575       IRN   = NP0
1576 *  +-------------------------------------------------------------------*
1577 *  |
1578       DO 10 I = 1,IR
1579 C
1580 C--------------------------------------------------
1581 C*** STORAGE OF H-N-COLL. PRODUCTS INTO NUCLEUS FINAL STATE C./FINUC/
1582 C*** CUT OFF FOR KINETIC ENERGIES LESS THEN ETHR, PUSH THEM INTO EX.EN.
1583 C*** TV
1584 C--------------------------------------------------
1585          T1 = EL  (I)
1586          I1 = ITR (I)
1587 *  |  +----------------------------------------------------------------*
1588 *  |  |
1589          IF ( T1 .LE. ETHR + AM (I1) ) THEN
1590             ANOW  = ANOW  + IBAR (I1)
1591             KTARN = KTARN - IBAR (I1) + ICH (I1)
1592             ZNOW  = ZNOW  + ICH (I1)
1593             KTARP = KTARP - ICH  (I1)
1594             AMMRES = ANOW * AMUAMU + 1.D-03 * FKENER ( ANOW, ZNOW )
1595             AMNRES = AMMRES - ZNOW * AMELEC + ELBNDE ( NINT (ZNOW) )
1596             EAVAIL = ECHCK - AMNRES + T1
1597             TPNVK  = TNK + TPK + TV
1598 *  |  |  +-------------------------------------------------------------*
1599 *  |  |  |
1600             IF ( TPNVK .GT. 0.D+00 ) THEN
1601                TKOR   = EAVAIL / TPNVK
1602                TPK = TPK * TKOR
1603                TNK = TNK * TKOR
1604                TV  = TV  * TKOR
1605                TVEUZ  = TVEUZ  * TKOR
1606                TVGRE0 = TVGRE0 * TKOR
1607 *  |  |  |
1608 *  |  |  +-------------------------------------------------------------*
1609 *  |  |  |
1610             ELSE
1611                TPK = 0.4444444444444444D+00 * EAVAIL
1612                TNK = 0.5555555555555556D+00 * EAVAIL
1613             END IF
1614 *  |  |  |
1615 *  |  |  +-------------------------------------------------------------*
1616 *  |  |
1617 *  |  +----------------------------------------------------------------*
1618 *  |  |                                  particle acceptable
1619          ELSE
1620             IRN = IRN+1
1621 C
1622 C--------------------------------------------------
1623 C*** STORE THE PARTICLE CHARACTERISTICS INTO C./FINUC/, THEN TAKE NEXT
1624 C*** PARTICLE
1625 C--------------------------------------------------
1626             ITRN(IRN) = I1
1627             CXRN(IRN) = CXR(I)
1628             CYRN(IRN) = CYR(I)
1629             CZRN(IRN) = CZR(I)
1630             ELR(IRN)  = EL(I)
1631             PLR(IRN)  = PL(I)
1632             IBAT = IBAT+1
1633 *  |  |   updating the conservation counters in common balanc
1634             ENUCR  = ENUCR  + ELR(IRN)
1635             PXNUCR = PXNUCR + PLR(IRN)*CXRN(IRN)
1636             PYNUCR = PYNUCR + PLR(IRN)*CYRN(IRN)
1637             PZNUCR = PZNUCR + PLR(IRN)*CZRN(IRN)
1638             ICNUCR = ICNUCR + ICH(ITRN(IRN))
1639             IBNUCR = IBNUCR + IBAR(ITRN(IRN))
1640             IRNTOT = IRNTOT + 1
1641          END IF
1642 *  |  |
1643 *  |  +----------------------------------------------------------------*
1644   10  CONTINUE
1645 *  |
1646 *  +-------------------------------------------------------------------*
1647
1648 *  +-------------------------------------------------------------------*
1649 *  |
1650       IF ( IRN - NP0 - IGREYP - IGREYN .EQ. 1 ) THEN
1651 *  |  The following balance is ok if and only if both the incident and
1652 *  |  the emitted particles were nucleons (proton or neutron), or the
1653 *  |  emitted particle is equal to the incident one
1654 *  |  +----------------------------------------------------------------*
1655 *  |  |
1656          IF ( ITTTT .NE. ITRN (NP0+1) ) THEN
1657             LPRONU = ITTTT .EQ. 1 .OR. ITTTT .EQ. 8
1658             LEMINU = ITRN (NP0+1) .EQ. 1 .OR. ITRN (NP0+1) .EQ. 8
1659 *  |  |  +-------------------------------------------------------------*
1660 *  |  |  |
1661             IF ( .NOT. LPRONU .OR. .NOT. LEMINU ) THEN
1662                ICNOW = NINT (ZNOW)
1663                IBNOW = NINT (ANOW)
1664                ICCC = ICHTAR + ICH  (ITTTT) - ICNUCR - ICINTR
1665                IBBB = IBTAR  + IBAR (ITTTT) - IBNUCR - IBINTR
1666 *  |  |  |  +----------------------------------------------------------*
1667 *  |  |  |  |
1668                IF ( IBBB .NE. IBNOW .AND. LEMINU .AND. IBAR (ITTTT)
1669      &              .EQ. 0 ) THEN
1670 *  |  |  |  |  +-------------------------------------------------------*
1671 *  |  |  |  |  |
1672                   IF ( ICCC .GT. ICNOW ) THEN
1673                      KTARP = KTARP - 1
1674                      KTARN = KTARN + 2
1675                      ANOW  = ANOW - 1.D+00
1676                      ZNOW  = ZNOW + 1.D+00
1677 *  |  |  |  |  |
1678 *  |  |  |  |  +-------------------------------------------------------*
1679 *  |  |  |  |  |
1680                   ELSE IF (ICCC .LT. ICNOW) THEN
1681                      KTARP = KTARP + 1
1682                      KTARN = KTARN
1683                      ANOW  = ANOW - 1.D+00
1684                      ZNOW  = ZNOW - 1.D+00
1685 *  |  |  |  |  |
1686 *  |  |  |  |  +-------------------------------------------------------*
1687 *  |  |  |  |  |
1688                   ELSE
1689                      KTARN = KTARN + 1
1690                      ANOW  = ANOW - 1.D+00
1691                   END IF
1692 *  |  |  |  |  |
1693 *  |  |  |  |  +-------------------------------------------------------*
1694 *  |  |  |  |
1695 *  |  |  |  +----------------------------------------------------------*
1696 *  |  |  |  |
1697                ELSE IF ( ICCC .NE. ICNOW .AND. LEMINU .AND. IBAR (ITTTT)
1698      &                   .EQ. 0 ) THEN
1699 *  |  |  |  |  +-------------------------------------------------------*
1700 *  |  |  |  |  |
1701                   IF ( ICCC .GT. ICNOW ) THEN
1702                      KTARP = KTARP - 1
1703                      KTARN = KTARN + 1
1704                      ZNOW  = ZNOW + 1.D+00
1705 *  |  |  |  |  |
1706 *  |  |  |  |  +-------------------------------------------------------*
1707 *  |  |  |  |  |
1708                   ELSE
1709                      KTARP = KTARP + 1
1710                      KTARN = KTARN - 1
1711                      ZNOW  = ZNOW - 1.D+00
1712                   END IF
1713 *  |  |  |  |  |
1714 *  |  |  |  |  +-------------------------------------------------------*
1715                END IF
1716 *  |  |  |  |
1717 *  |  |  |  +----------------------------------------------------------*
1718             END IF
1719 *  |  |  |
1720 *  |  |  +-------------------------------------------------------------*
1721          END IF
1722 *  |  |
1723 *  |  +----------------------------------------------------------------*
1724       END IF
1725 *  |
1726 *  +-------------------------------------------------------------------*
1727       GO TO 406
1728 *  +-->-->-->-->--> go to the cascade simulation!!!
1729 C
1730 C--------------------------------------------------
1731 C*** IF ITTTT IS A STOPPED NEGATIVE MESON,(BUT ALSO FOR ANNIHILATION)
1732 C*** DISTRIBUTE THE AVAILABLE TOTAL ENERGY OF THE PARTICLE TO THE CASCAD
1733 C*** AND EXCITATION
1734 C--------------------------------------------------
1735 *  +-------------------------------------------------------------------*
1736 *  | Zero secondary part. case (h-n-coll.), stopping mesons or anyway
1737 *  | all situations not allowing a call to Hadrin
1738 405   CONTINUE
1739          V0WELL (ITJ) = V0OLD
1740          IRN = NP0
1741          AMMRES = AMUAMU * ANOW + 1.D-03 * FKENER ( ANOW, ZNOW )
1742          AMNRES = AMMRES - ZNOW * AMELEC + ELBNDE ( NINT (ZNOW) )
1743          TPNVK  = TPK + TNK + TV
1744          EAVAIL = ETTOT - AMNRES - TPNVK
1745 *  |  +----------------------------------------------------------------*
1746 *  |  |
1747          IF ( TPNVK .GT. 0.D+00 ) THEN
1748             TKOR = 1.D+00 + EAVAIL / TPNVK
1749             TPK  = TPK * TKOR
1750             TNK  = TNK * TKOR
1751             TV     = TV * TKOR
1752             TVGRE0 = TV
1753             TVEUZ  = 0.D+00
1754 *  |  |
1755 *  |  +----------------------------------------------------------------*
1756 *  |  |
1757          ELSE
1758             TPK  = 0.4D+00 * EAVAIL
1759             TNK  = 0.5D+00 * EAVAIL
1760             TVGRE0 = 0.1D+00 * EAVAIL
1761             TV     = TVGRE0
1762             TVEUZ  = 0.D+00
1763          END IF
1764 *  |  |
1765 *  |  +----------------------------------------------------------------*
1766          ELPP = 0.0D+00
1767 *  |
1768 *  +-------------------------------------------------------------------*
1769 *  Now the cascade simulation!!!
1770 406   CONTINUE
1771       EINCP = 0.D+00
1772       EINCN = 0.D+00
1773       IU  = IRN - NP0
1774 C
1775 C--------------------------------------------------
1776 C
1777 C   SIMULATION OF CASCADE NEUTRONS
1778 C
1779 C--------------------------------------------------
1780 C
1781 C--------------------------------------------------
1782 C*** NUCLEON NUMBER AND CASCADE ENERGY(VERSUS CUT OFF)-TEST
1783 C*** GIVE THE ENERGY TO  THE EXCITATION , IF TEST DEMANDS THIS
1784 C--------------------------------------------------
1785 *  +-------------------------------------------------------------------*
1786 *  |
1787       IF ((TPK.LE.TPKTE).OR.(ZNOW.LT.0.5D0)) THEN
1788          TV=TV+TPK
1789          TPK = 0.D+00
1790          LPCASC = .FALSE.
1791 *  |
1792 *  +-------------------------------------------------------------------*
1793 *  |
1794       ELSE
1795          LPCASC = .TRUE.
1796       END IF
1797 *  |
1798 *  +-------------------------------------------------------------------*
1799 *  +-------------------------------------------------------------------*
1800 *  |
1801       IF ((TNK.LE.TNKTE).OR.(ANOW-ZNOW.LT.0.5D0)) THEN
1802          TV=TV+TNK
1803          TNK = 0.D+00
1804          LNCASC = .FALSE.
1805 *  |
1806 *  +-------------------------------------------------------------------*
1807 *  |
1808       ELSE
1809          LNCASC = .TRUE.
1810       END IF
1811 *  |
1812 *  +-------------------------------------------------------------------*
1813       KNREJE = 0
1814       KPREJE = 0
1815 *  +-------------------------------------------------------------------*
1816 *  |    Cascade nucleons selection !!!!
1817 1900  CONTINUE
1818 *  |  +----------------------------------------------------------------*
1819 *  |  |
1820          IF ( LPCASC ) THEN
1821             ZSAMP = ZNOW
1822             IF ( ZSAMP .LE. 0.5D+00 ) THEN
1823                LPCASC = .FALSE.
1824             END IF
1825 *  |  |
1826 *  |  +----------------------------------------------------------------*
1827 *  |  |
1828          ELSE
1829             ZSAMP = 0.D+00
1830          END IF
1831 *  |  |
1832 *  |  +----------------------------------------------------------------*
1833 *  |  +----------------------------------------------------------------*
1834 *  |  |
1835          IF ( LNCASC ) THEN
1836             ANSMP = ANOW - ZNOW
1837             IF ( ANSMP .LE. 0.5D+00 ) THEN
1838                LNCASC = .FALSE.
1839             END IF
1840 *  |  |
1841 *  |  +----------------------------------------------------------------*
1842 *  |  |
1843          ELSE
1844             ANSMP = 0.D+00
1845          END IF
1846 *  |  |
1847 *  |  +----------------------------------------------------------------*
1848          IF ( ( .NOT. LPCASC ) .AND. ( .NOT. LNCASC ) ) GO TO 2000
1849          IF ( NINT (ANOW) .LE. 0 ) THEN
1850             WRITE(LUNERR,*)' Nucriv: <<<< ANOW < 0 >>>>',
1851      &      ANOW,ZNOW,LPCASC,LNCASC
1852             LPCASC = .FALSE.
1853             LNCASC = .FALSE.
1854             GO TO 2000
1855          END IF
1856          CALL GRNDM(RNDM,1)
1857          IF ( RNDM (1) .LT. ANSMP / ( ANSMP + ZSAMP ) ) THEN
1858             GO TO 20
1859          ELSE
1860             GO TO 30
1861          END IF
1862    20    CONTINUE
1863          T2=AM(8)
1864 C
1865 C--------------------------------------------------
1866 C*** NEUTRON NUMBER AND NEUTRON CASCADE ENERGY TEST
1867 C--------------------------------------------------
1868          IF ( ANSMP .LT. 0.5D+00 .OR. TNK .LE. TNKTE ) THEN
1869             LNCASC = .FALSE.
1870             TV  = TV + TNK
1871             TNK = 0.D+00
1872             GO TO 1900
1873          END IF
1874 C
1875 C--------------------------------------------------
1876 C*** SAMPLING OF KINETIC ENERGY TN AND TETA ANGLE DN OF THE EMITTED NEU
1877 C--------------------------------------------------
1878          CALL RBKEKV ( 2, EXSOP, TOEFF, BBTAR, TKIN, TSTRCK,
1879      &                 PSTRCK, ANOW, TKRECL, COD, SID )
1880          TN  = TKIN
1881          TNK = TNK-TN
1882 C
1883 C--------------------------------------------------
1884 C*** TEST: IF THE REMAINING ENERGY AND CORRESP.NUCLEON NUMBER ALLOWS FUR
1885 C*** THER CORRESP. NUCLEON EMISSION, ELSE GIVE THE  THE ENERGY OR TAKE
1886 C*** THE MISSING FROM THE EXCITATION
1887 C*** IF FOR LAST CHOSEN NUCLEON TN>REMAINING REST, TAKE IN 90% THE ENER
1888 C*** DIFFERENCE FROM EXCIT.ENERGY, IN 10% OR IF TV WOULD BE LESS THEN 0.
1889 C***  ADD TN TO TV (It does the opposite!!! A. Ferrari)
1890 C--------------------------------------------------
1891          IF ( ANSMP .LT. 1.5D+00 .OR. TNK .LE. ETHR  ) THEN
1892             LNCASC = .FALSE.
1893             CALL GRNDM(RNDM,1)
1894             VT = RNDM (1)
1895             THK=TV+TNK
1896             IF (VT .GT. 0.9D+00 .OR. THK .LT. 0.D+00) THEN
1897                TV  = TV + TNK + TN
1898                TNK = 0.D+00
1899                GO TO 1900
1900             END IF
1901             TNK=0.D0
1902             TV=THK
1903          END IF
1904 C
1905 C--------------------------------------------------
1906 C*** CORRESP. NUCLEON NUMBER COUNTING
1907 C--------------------------------------------------
1908          IF ( TSTRCK .LT. ETHR ) THEN
1909             TV = TV + TN
1910             GO TO 1900
1911 *  |  |
1912 *  +-<|--<--<--<--<--< Try to select another nucleon!!
1913          END IF
1914          ANOW   = ANOW - 1.D+00
1915          FRSURV = ANOW / BBTAR
1916 C
1917 C--------------------------------------------------
1918 C*** TAKE THE EMITTED CORRESP. NUCLEON INTO THE FINAL PARTICLE TABLE
1919 C*** C./FINUC/ AFTER CHOICE OF PHI AND TURNING INTO THE LAB SYSTEM
1920 C--------------------------------------------------
1921          IRN = IRN + 1
1922          CALL COSI ( SFE, CFE )
1923 *  |  +----------------------------------------------------------------*
1924 *  |  |
1925          IF ( ANOW .LT. 30.D+00 .OR. FRSURV .LT. 0.66D+00 ) THEN
1926             PXLAST = PXTTOT - PXINTR - PXNUCR
1927             PYLAST = PYTTOT - PYINTR - PYNUCR
1928             PZLAST = PZTTOT - PZINTR - PZNUCR
1929             PPLAS2 = PXLAST**2 + PYLAST**2 + PZLAST**2
1930 *  |  |  +-------------------------------------------------------------*
1931 *  |  |  |
1932             IF ( PPLAS2 .GT. 1.D-6 * PTTOT**2 ) THEN
1933                PPLAST = SQRT ( PPLAS2 )
1934                CXXINC = PXLAST / PPLAST
1935                CYYINC = PYLAST / PPLAST
1936                CZZINC = PZLAST / PPLAST
1937                CDMIN = MIN ( ONEONE, PSTRCK / PPLAST ) * ( 1.D+00 -
1938      &                       0.8D+00 * FRSURV )
1939                COD = 1.D+00 - 0.5D+00 * ( 1.D+00 - COD ) * ( 1.D+00 -
1940      &               CDMIN )
1941                COD = MIN ( COD, ONEONE )
1942                SID = SQRT ( 1.D+00 - COD**2 )
1943 *  |  |  |
1944 *  |  |  +-------------------------------------------------------------*
1945 *  |  |  |
1946             ELSE
1947                CXXINC = CX
1948                CYYINC = CY
1949                CZZINC = CZ
1950             END IF
1951 *  |  |  |
1952 *  |  |  +-------------------------------------------------------------*
1953             CALL TTRANS ( CXXINC, CYYINC, CZZINC, COD, SID, SFE, CFE,
1954      &                    CXRN (IRN), CYRN (IRN), CZRN (IRN) )
1955 *  |  |  +-------------------------------------------------------------*
1956 *  |  |  |
1957             IF ( NINT (ANOW) .GT. 0 ) THEN
1958                AMNRES = AMUAMU * ANOW - ZNOW * AMELEC + 1.D-03 * FKENER
1959      &                ( ANOW, ZNOW ) + ELBNDE ( NINT ( ZNOW ) )
1960 *  |  |  |
1961 *  |  |  +-------------------------------------------------------------*
1962 *  |  |  |
1963             ELSE
1964                AMNRES = 0.D+00
1965             END IF
1966 *  |  |  |
1967 *  |  |  +-------------------------------------------------------------*
1968             AMNRE2 = AMNRES  * AMNRES
1969             ELEFT  = ETTOT - ENUCR - EINTR - TSTRCK - T2
1970 *  |  | Update TV !!!
1971             TV     = ELEFT - TPK - TNK - AMNRES
1972             IUMO   = 0
1973 *  |  |  +-------------------------------------------------------------*
1974 *  |  |  |
1975 2020        CONTINUE
1976                PXLEFT = PXLAST - PSTRCK * CXRN (IRN)
1977                PYLEFT = PYLAST - PSTRCK * CYRN (IRN)
1978                PZLEFT = PZLAST - PSTRCK * CZRN (IRN)
1979                PPLEF2 = PXLEFT**2 + PYLEFT**2 + PZLEFT**2
1980                UMO2   = ELEFT**2 - PPLEF2
1981                DELTU2 = AMNRE2 - UMO2
1982 *  |  |  |  +----------------------------------------------------------*
1983 *  |  |  |  |
1984                IF ( DELTU2 .GT. 0.D+00 ) THEN
1985 *  |  |  |  |  +-------------------------------------------------------*
1986 *  |  |  |  |  |
1987                   IF ( IUMO .EQ. 0 .AND. PSTRCK .GE. PPLAST ) THEN
1988                      CXRN (IRN) = CXXINC
1989                      CYRN (IRN) = CYYINC
1990                      CZRN (IRN) = CZZINC
1991                      IUMO   = IUMO + 1
1992                      GO TO 2020
1993                   END IF
1994 *  |  |  |  |  |
1995 *  |  |  |  |  +-------------------------------------------------------*
1996                   PLDOTU = PXLEFT * CXRN (IRN) + PYLEFT * CYRN (IRN) +
1997      &                     PZLEFT * CZRN (IRN)
1998                   DELTAE = 0.51D+00 * DELTU2 / ( ELEFT - ( T2 + TSTRCK )
1999      &                   * PLDOTU / PSTRCK )
2000 *  |  |  |  |  +-------------------------------------------------------*
2001 *  |  |  |  |  |
2002                   IF ( DELTAE .GE. TSTRCK .OR. IUMO .GT. 3 ) THEN
2003                      IRN   = IRN - 1
2004                      KNREJE = KNREJE + 1
2005                      ANOW  = ANOW + 1.D+00
2006 *  |  |  |  |  |  +----------------------------------------------------*
2007 *  |  |  |  |  |  |
2008                      IF ( KNREJE .GT. 3 ) THEN
2009                         LNCASC = .FALSE.
2010 *  |  |  |  |  |  |  +-------------------------------------------------*
2011 *  |  |  |  |  |  |  |
2012                         IF ( LPCASC ) THEN
2013                            TPK = TPK + TNK + TN
2014 *  |  |  |  |  |  |  |
2015 *  |  |  |  |  |  |  +-------------------------------------------------*
2016 *  |  |  |  |  |  |  |
2017                         ELSE
2018                            TV = TV + TNK + TN
2019                         END IF
2020 *  |  |  |  |  |  |  |
2021 *  |  |  |  |  |  |  +-------------------------------------------------*
2022                         TNK = 0.D+00
2023 *  |  |  |  |  |  |
2024 *  |  |  |  |  |  +----------------------------------------------------*
2025 *  |  |  |  |  |  |
2026                      ELSE
2027                         TNK = TNK + TN
2028                      END IF
2029 *  |  |  |  |  |  |
2030 *  |  |  |  |  |  +----------------------------------------------------*
2031                      GO TO 1900
2032 *  |  |  |  |  |
2033 *  +-<|-<|-<|-<|--<--<--<--<--<
2034 *  |  |  |  |  |
2035 *  |  |  |  |  +-------------------------------------------------------*
2036 *  |  |  |  |  |
2037                   ELSE
2038                      TSTRCK = TSTRCK - DELTAE
2039                      PSTRCK = SQRT ( TSTRCK * ( TSTRCK + 2.D+00 * T2 ) )
2040                      TKRECL = TKRECL + DELTAE
2041                      ELEFT  = ELEFT  + DELTAE
2042                      GO TO 2020
2043 *  |  |  |  |
2044 *  |  |  +-<|-<|--<--<--<--<
2045                   END IF
2046 *  |  |  |  |  |
2047 *  |  |  |  |  +-------------------------------------------------------*
2048                END IF
2049 *  |  |  |  |
2050 *  |  |  |  +----------------------------------------------------------*
2051 *  |  |  |
2052 *  |  |  +-------------------------------------------------------------*
2053 *  |  |
2054 *  |  +----------------------------------------------------------------*
2055 *  |  |
2056          ELSE
2057             CALL TTRANS ( CX, CY, CZ, COD, SID, SFE, CFE, CXRN (IRN),
2058      &                    CYRN (IRN), CZRN (IRN) )
2059          END IF
2060 *  |  |
2061 *  |  +----------------------------------------------------------------*
2062          IGREYN = IGREYN + 1
2063          EINCN  = EINCN + TKIN
2064          ITRN(IRN) = 8
2065 *
2066 *  T2 is the mass energy of the neutron (see before)
2067 *
2068          TVGRE0 = TVGRE0
2069          TVGREY = TVGREY + TKIN + TKRECL - TSTRCK - EBNDNG (2)
2070          EEX = TKIN + TKRECL - TSTRCK - EBNDNG (2)
2071          TV  = TV + EEX
2072          ELR (IRN)  = TSTRCK + T2
2073          PLR(IRN) = PSTRCK
2074 *  |   updating the conservation counters in common balanc
2075          EINTR  = EINTR + ELR(IRN)
2076          PXINTR = PXINTR + PLR(IRN)*CXRN(IRN)
2077          PYINTR = PYINTR + PLR(IRN)*CYRN(IRN)
2078          PZINTR = PZINTR + PLR(IRN)*CZRN(IRN)
2079          ICINTR = ICINTR + ICH(ITRN(IRN))
2080          IBINTR = IBINTR + IBAR(ITRN(IRN))
2081          IRNTOT = IRNTOT + 1
2082       GO TO 1900
2083 *  |
2084 *  +--<--<--<--<--< Try to select another cascade nucleon!!!!
2085
2086 *  +-------------------------------------------------------------------*
2087 *  |    Cascade protons selection!!!!!
2088    30 CONTINUE
2089          IF ( ZNOW .LT. 0.5D+00 .OR. TPK .LE. TPKTE ) THEN
2090             LPCASC = .FALSE.
2091             TV  = TV + TPK
2092             TPK = 0.D+00
2093             GO TO 1900
2094          END IF
2095 C
2096 C--------------------------------------------------
2097 C
2098 C   SIMULATION OF CASCADE PROTONS
2099 C
2100 C--------------------------------------------------
2101 C*** PROTON NUMBER  AND PROTON CASCADE ENERGY TEST
2102 C
2103 C
2104 C--------------------------------------------------
2105 C*** SAMPLING OF KINETIC ENERGY TN AND TETA ANGLE DN OF THE EMITTED P
2106 C--------------------------------------------------
2107 C--------------------------------------------------
2108          T2 = AM(1)
2109          CALL RBKEKV ( 1, EXSOP, TOEFF, BBTAR, TKIN, TSTRCK,
2110      &                 PSTRCK, ANOW, TKRECL, COD, SID )
2111          TN = TKIN
2112          TPK = TPK - TN
2113 C
2114 C--------------------------------------------------
2115 C*** TEST: IF THE REMAINING ENERGY AND CORRESP.NUCLEON NUMBER ALLOWS FOR
2116 C*** THE CORRESP. NUCLEON EMISSION, ELSE GIVE THE ENERGY OR TAKE
2117 C*** THE MISSING FROM THE EXCITATION
2118 C*** IF FOR LAST CHOSEN NUCLEON TN>REMAINING REST, TAKE IN 90% THE ENERG
2119 C*** DIFFERENCE FROM EXCIT.ENERGY, IN 10% OR IF TV WOULD BE LESS THEN 0.
2120 C***  ADD TN TO TV
2121 C--------------------------------------------------
2122          IF ( ZNOW .LT. 1.5D+00 .OR. TPK .LE. ETHR  ) THEN
2123             LPCASC = .FALSE.
2124             CALL GRNDM(RNDM,1)
2125             VT  = RNDM(1)
2126             THK = TV + TPK
2127             IF ( VT .GT. 0.9D+00 .OR. THK .LT. 0.D+00 ) THEN
2128                TV=TV+TPK+TN
2129                TPK = 0.D+00
2130                GO TO 1900
2131             END IF
2132             TPK = 0.D0
2133             TV  = THK
2134          END IF
2135 C
2136 C--------------------------------------------------
2137 C*** CORRESP. NUCLEON NUMBER COUNTING
2138 C--------------------------------------------------
2139          IF (TN .LT. ETHR) THEN
2140             TV = TV + TN
2141             GO TO 1900
2142 *  |  |
2143 *  +-<|--<--<--<--<--< try to select another nucleon
2144          END IF
2145          ANOW = ANOW - 1.D0
2146          ZNOW = ZNOW - 1.D0
2147          FRSURV = ANOW / BBTAR
2148 C
2149 C--------------------------------------------------
2150 C*** TAKE THE EMITTED CORRESP. NUCLEON INTO THE FINAL PARTICLE TABLE
2151 C*** C./FINUC/ AFTER CHOICE OF PHI AND TURNING INTO THE LAB SYSTEM
2152 C--------------------------------------------------
2153          IRN = IRN + 1
2154          CALL COSI ( SFE, CFE )
2155 *  |  +----------------------------------------------------------------*
2156 *  |  |
2157          IF ( ANOW .LT. 30.D+00 .OR. FRSURV .LT. 0.66D+00 ) THEN
2158             PXLAST = PXTTOT - PXINTR - PXNUCR
2159             PYLAST = PYTTOT - PYINTR - PYNUCR
2160             PZLAST = PZTTOT - PZINTR - PZNUCR
2161             PPLAS2 = PXLAST**2 + PYLAST**2 + PZLAST**2
2162 *  |  |  +-------------------------------------------------------------*
2163 *  |  |  |
2164             IF ( PPLAS2 .GT. 1.D-6 * PTTOT**2 ) THEN
2165                PPLAST = SQRT ( PPLAS2 )
2166                CXXINC = PXLAST / PPLAST
2167                CYYINC = PYLAST / PPLAST
2168                CZZINC = PZLAST / PPLAST
2169                CDMIN = MIN ( ONEONE, PSTRCK / PPLAST ) * ( 1.D+00 -
2170      &                       0.8D+00 * FRSURV )
2171                COD = 1.D+00 - 0.5D+00 * ( 1.D+00 - COD ) * ( 1.D+00 -
2172      &               CDMIN )
2173                COD = MIN ( COD, ONEONE )
2174                SID = SQRT ( 1.D+00 - COD**2 )
2175 *  |  |  |
2176 *  |  |  +-------------------------------------------------------------*
2177 *  |  |  |
2178             ELSE
2179                CXXINC = CX
2180                CYYINC = CY
2181                CZZINC = CZ
2182             END IF
2183 *  |  |  |
2184 *  |  |  +-------------------------------------------------------------*
2185             CALL TTRANS ( CXXINC, CYYINC, CZZINC, COD, SID, SFE, CFE,
2186      &                    CXRN (IRN), CYRN (IRN), CZRN (IRN) )
2187 *  |  |  +-------------------------------------------------------------*
2188 *  |  |  |
2189             IF ( NINT (ANOW) .GT. 0 ) THEN
2190                AMNRES = AMUAMU * ANOW - ZNOW * AMELEC + 1.D-03 * FKENER
2191      &                ( ANOW, ZNOW ) + ELBNDE ( NINT ( ZNOW ) )
2192 *  |  |  |
2193 *  |  |  +-------------------------------------------------------------*
2194 *  |  |  |
2195             ELSE
2196                AMNRES = 0.D+00
2197             END IF
2198 *  |  |  |
2199 *  |  |  +-------------------------------------------------------------*
2200             AMNRE2 = AMNRES  * AMNRES
2201             ELEFT  = ETTOT  - ENUCR  - EINTR  - TSTRCK - T2
2202             IUMO   = 0
2203 *  |  |  +-------------------------------------------------------------*
2204 *  |  |  |
2205 3030        CONTINUE
2206                PXLEFT = PXLAST - PSTRCK * CXRN (IRN)
2207                PYLEFT = PYLAST - PSTRCK * CYRN (IRN)
2208                PZLEFT = PZLAST - PSTRCK * CZRN (IRN)
2209                PPLEF2 = PXLEFT**2 + PYLEFT**2 + PZLEFT**2
2210                UMO2   = ELEFT**2 - PPLEF2
2211                DELTU2 = AMNRE2 - UMO2
2212 *  |  |  |  +----------------------------------------------------------*
2213 *  |  |  |  |
2214                IF ( DELTU2 .GT. 0.D+00 ) THEN
2215 *  |  |  |  |  +-------------------------------------------------------*
2216 *  |  |  |  |  |
2217                   IF ( IUMO .EQ. 0 .AND. PSTRCK .GE. PPLAST ) THEN
2218                      CXRN (IRN) = CXXINC
2219                      CYRN (IRN) = CYYINC
2220                      CZRN (IRN) = CZZINC
2221                      IUMO = IUMO + 1
2222                      GO TO 3030
2223                   END IF
2224 *  |  |  |  |  |
2225 *  |  |  |  |  +-------------------------------------------------------*
2226                   IUMO   = IUMO + 1
2227                   PLDOTU = PXLEFT * CXRN (IRN) + PYLEFT * CYRN (IRN) +
2228      &                     PZLEFT * CZRN (IRN)
2229                   DELTAE = 0.51D+00 * DELTU2 / ( ELEFT - ( T2
2230      &                   + TSTRCK ) * PLDOTU / PSTRCK )
2231 *  |  |  |  |  +-------------------------------------------------------*
2232 *  |  |  |  |  |
2233                   IF ( DELTAE .GE. TSTRCK .OR. IUMO .GT. 3 ) THEN
2234                      IRN   = IRN - 1
2235                      KPREJE = KPREJE + 1
2236                      ANOW  = ANOW + 1.D+00
2237                      ZNOW  = ZNOW + 1.D+00
2238 *  |  |  |  |  |  +----------------------------------------------------*
2239 *  |  |  |  |  |  |
2240                      IF ( KPREJE .GT. 3 ) THEN
2241                         LPCASC = .FALSE.
2242 *  |  |  |  |  |  |  +-------------------------------------------------*
2243 *  |  |  |  |  |  |  |
2244                         IF ( LNCASC ) THEN
2245                            TNK = TPK + TNK + TN
2246 *  |  |  |  |  |  |  |
2247 *  |  |  |  |  |  |  +-------------------------------------------------*
2248 *  |  |  |  |  |  |  |
2249                         ELSE
2250                            TV = TV + TPK + TN
2251                         END IF
2252 *  |  |  |  |  |  |  |
2253 *  |  |  |  |  |  |  +-------------------------------------------------*
2254                         TPK = 0.D+00
2255 *  |  |  |  |  |  |
2256 *  |  |  |  |  |  +----------------------------------------------------*
2257 *  |  |  |  |  |  |
2258                      ELSE
2259                         TPK = TPK + TN
2260                      END IF
2261 *  |  |  |  |  |  |
2262 *  |  |  |  |  |  +----------------------------------------------------*
2263                      GO TO 1900
2264 *  |  |  |  |  |
2265 *  +-<|-<|-<|-<|--<--<--<--<--<
2266 *  |  |  |  |  |
2267 *  |  |  |  |  +-------------------------------------------------------*
2268 *  |  |  |  |  |
2269                   ELSE
2270                      TSTRCK = TSTRCK - DELTAE
2271                      PSTRCK = SQRT ( TSTRCK * ( TSTRCK + 2.D+00 * T2 ) )
2272                      TKRECL = TKRECL + DELTAE
2273                      ELEFT  = ELEFT  + DELTAE
2274                      GO TO 3030
2275 *  |  |  |  |  |
2276 *  |  |  +-<|-<|--<--<--<--<
2277                   END IF
2278 *  |  |  |  |  |
2279 *  |  |  |  |  +-------------------------------------------------------*
2280                END IF
2281 *  |  |  |  |
2282 *  |  |  |  +----------------------------------------------------------*
2283 *  |  |  |
2284 *  |  |  +-------------------------------------------------------------*
2285 *  |  |
2286 *  |  +----------------------------------------------------------------*
2287 *  |  |
2288          ELSE
2289             CALL TTRANS ( CX, CY, CZ, COD, SID, SFE, CFE, CXRN (IRN),
2290      &                    CYRN (IRN), CZRN (IRN) )
2291          END IF
2292 *  |  |
2293 *  |  +----------------------------------------------------------------*
2294          IGREYP = IGREYP + 1
2295          EINCP  = EINCP + TKIN
2296 C
2297          ITRN (IRN) = 1
2298 *
2299 *  T2 is the mass energy of the proton (see before)
2300 *
2301 *        TVGRE0 = TVGRE0 + EXSOP (1)
2302          TVGRE0 = TVGRE0
2303          TVGREY = TVGREY + TKIN + TKRECL - TSTRCK - EBNDNG (1)
2304          EEX = TKIN + TKRECL - TSTRCK - EBNDNG (1)
2305          TV  = TV + EEX
2306          ELR (IRN) = TSTRCK + T2
2307          PLR (IRN) = PSTRCK
2308 *  |   updating the conservation counters in common balanc
2309          EINTR  = EINTR + ELR(IRN)
2310          PXINTR = PXINTR + PLR(IRN)*CXRN(IRN)
2311          PYINTR = PYINTR + PLR(IRN)*CYRN(IRN)
2312          PZINTR = PZINTR + PLR(IRN)*CZRN(IRN)
2313          ICINTR = ICINTR + ICH(ITRN(IRN))
2314          IBINTR = IBINTR + IBAR(ITRN(IRN))
2315          IRNTOT = IRNTOT + 1
2316       GO TO 1900
2317 *  |
2318 *  +--<--<--<--<--< Try to select another cascade nucleon!!!!
2319 C
2320 C--------------------------------------------------
2321 C*** END OF CASCADE NUCLEON EMISSION
2322 C*** HANDLING OF THE CASES: NO OR ONE PARTICLE IN H-A-FINAL STATE
2323 C*** FINAL CALCULATION OF TV
2324 C--------------------------------------------------
2325 2000  CONTINUE
2326 *  +-------------------------------------------------------------------*
2327 *  |   Now working with atomic masses !
2328       IF ( NINT (ANOW) .GT. 0 ) THEN
2329          AMMRES = AMUAMU * ANOW + 1.D-03 * FKENER ( ANOW, ZNOW )
2330          AMNRES = AMMRES - ZNOW * AMELEC + ELBNDE ( NINT (ZNOW) )
2331          DELEFT = AMMTAR - AMNTAR - ( ICHTAR - ZNOW ) * AMELEC
2332 *  |
2333 *  +-------------------------------------------------------------------*
2334 *  |
2335       ELSE
2336          AMMRES = 0.D+00
2337          AMNRES = 0.D+00
2338          DELEFT = 0.D+00
2339       END IF
2340 *  |
2341 *  +-------------------------------------------------------------------*
2342 *  +-------------------------------------------------------------------*
2343 *  |                        Check the number of emitted particles
2344       IF ( IRN - NP0 - 1 ) 4001,4002,4003
2345 4001  CONTINUE
2346 *  |  No particle emitted !!!
2347 *  |  Now working with atomic masses
2348          TV = ETTOT - AMNRES
2349       GO TO 4013
2350 *  |
2351 *  +-------------------------------------------------------------------*
2352 *  |  One particle emitted
2353 4002  CONTINUE
2354 *  |  +----------------------------------------------------------------*
2355 *  |  |
2356          IF ( TV .LE. 0.001D+00 ) THEN
2357             TV  = TVEUZ0
2358          END IF
2359 *  |  |
2360 *  |  +----------------------------------------------------------------*
2361 *  |   Now working with atomic masses
2362          UMO2   = ( ETTOT + DELEFT )**2 - PTTOT**2
2363          IUMO   = 0
2364 *  |  +----------------------------------------------------------------*
2365 *  |  |
2366 4040     CONTINUE
2367             AMSTAR = AMMRES + TV
2368             AMEMIT = AM ( ITRN (NP0+1) )
2369             UMIN2  = ( AMEMIT + AMSTAR )**2
2370 *  |  |  +-------------------------------------------------------------*
2371 *  |  |  |
2372             IF ( UMO2 .LE. UMIN2 ) THEN
2373                DELTAE = 0.51D+00 * ( UMIN2 - UMO2 ) / ( AMEMIT +
2374      &                  AMSTAR )
2375                IUMO = IUMO + 1
2376 *  |  |  |  +----------------------------------------------------------*
2377 *  |  |  |  |
2378                IF ( TV .GE. DELTAE .AND. IUMO .LE. 3 ) THEN
2379                   TV = TV - DELTAE
2380                   GO TO 4040
2381 *  |  |  |  |
2382 *  |  +-<|-<|--<--<--<--<--<
2383 *  |  |  |  |
2384 *  |  |  |  +----------------------------------------------------------*
2385 *  |  |  |  |
2386                ELSE
2387                   WRITE ( LUNERR, * ) ' Nucriv: single particle emiss',
2388      &                  'ion not possible for missing invariant mass!!',
2389      &                     NINT (ANOW), NINT (ZNOW), UMIN2, UMO2,
2390      &                     ITRN (NP0+1), IUMO
2391                   ANOW = BBTAR
2392                   ZNOW = ZZTAR
2393                   KTARP = 0
2394                   KTARN = 0
2395                   IGREYP = 0
2396                   IGREYN = 0
2397                   ENUCR  = 0.D+00
2398                   PXNUCR = 0.D+00
2399                   PYNUCR = 0.D+00
2400                   PZNUCR = 0.D+00
2401                   EINTR   = 0.D+00
2402                   PXINTR  = 0.D+00
2403                   PYINTR  = 0.D+00
2404                   PZINTR  = 0.D+00
2405                   TVGREY  = 0.D+00
2406                   UMO     = SQRT ( UMO2 )
2407                   TVGRE0  = UMO - AMMTAR
2408                   TV  = TVGRE0
2409                   IRN = NP0
2410                   RETURN
2411                END IF
2412 *  |  |  |  |
2413 *  |  |  |  +----------------------------------------------------------*
2414             END IF
2415 *  |  |  |
2416 *  |  |  +-------------------------------------------------------------*
2417 *  |  |
2418 *  |  +----------------------------------------------------------------*
2419          UMO    = SQRT ( UMO2 )
2420          GAMCM  = ( ETTOT + DELEFT ) / UMO
2421          ETACM  = PTTOT / UMO
2422          ETAX   = PXTTOT / UMO
2423          ETAY   = PYTTOT / UMO
2424          ETAZ   = PZTTOT / UMO
2425          PXNUCR = PLR (NP0+1) * CXRN (NP0+1)
2426          PYNUCR = PLR (NP0+1) * CYRN (NP0+1)
2427          PZNUCR = PLR (NP0+1) * CZRN (NP0+1)
2428          ETAPCM = - ETAX * PXNUCR - ETAY * PYNUCR - ETAZ * PZNUCR
2429          PHELP  = ETAPCM / ( GAMCM + 1.D+00 ) + ELR (NP0+1)
2430 *  |  Old direction cosines of the emitted particle in CMS
2431          PXINTR = PXNUCR - ETAX * PHELP
2432          PYINTR = PYNUCR - ETAY * PHELP
2433          PZINTR = PZNUCR - ETAZ * PHELP
2434          PCMS   = SQRT ( PXINTR**2 + PYINTR**2 + PZINTR**2 )
2435          CXXINC = PXINTR / PCMS
2436          CYYINC = PYINTR / PCMS
2437          CZZINC = PZINTR / PCMS
2438          ECMSST = 0.5D+00 * ( UMO2 + AMSTAR**2 - AMEMIT**2 ) / UMO
2439          ECMSEM = UMO - ECMSST
2440          PCMS   = SQRT ( ECMSEM**2 - AMEMIT**2 )
2441 *  |   Now save the old direction
2442          PXINTR = PCMS * CXXINC
2443          PYINTR = PCMS * CYYINC
2444          PZINTR = PCMS * CZZINC
2445          ETAPCM = ETAX * PXINTR + ETAY * PYINTR + ETAZ * PZINTR
2446          ELR (NP0+1) = GAMCM * ECMSEM + ETAPCM
2447          PHELP   = ETAPCM / ( GAMCM + 1.D+00 ) + ECMSEM
2448          PXNUCR  = PXINTR + ETAX * PHELP
2449          PYNUCR  = PYINTR + ETAY * PHELP
2450          PZNUCR  = PZINTR + ETAZ * PHELP
2451          PLR (NP0+1) = SQRT ( PXNUCR**2 + PYNUCR**2 + PZNUCR**2 )
2452          ENUCR   = ELR  (NP0+1)
2453          CXRN (NP0+1) = PXNUCR / PLR (NP0+1)
2454          CYRN (NP0+1) = PYNUCR / PLR (NP0+1)
2455          CZRN (NP0+1) = PZNUCR / PLR (NP0+1)
2456          KTARP   = KTARP + IGREYP
2457          KTARN   = KTARN + IGREYN
2458          IGREYP  = 0
2459          IGREYN  = 0
2460          EINTR   = 0.D+00
2461          PXINTR  = 0.D+00
2462          PYINTR  = 0.D+00
2463          PZINTR  = 0.D+00
2464          TVGRE0  = 0.D+00
2465          TVGREY  = 0.D+00
2466          TVEUZ   = TV
2467       GO TO 4013
2468 *  |
2469 *  +-------------------------------------------------------------------*
2470
2471 *  +-------------------------------------------------------------------*
2472 *  |  More than one particle emitted!!
2473 4003  CONTINUE
2474 *  |  Now working with atomic masses !
2475          TV = ETTOT + DELEFT - AMMRES - ENUCR - EINTR
2476 *  |  +----------------------------------------------------------------*
2477 *  |  |
2478          IF ( TV .LE. 0.D+00 ) THEN
2479             WRITE ( LUNERR, * )' *** Nucrin failure: Tv < 0!!',TV
2480          END IF
2481 *  |  |
2482 *  |  +----------------------------------------------------------------*
2483  4013 CONTINUE
2484 *  |
2485 *  +-------------------------------------------------------------------*
2486 *  +-------------------------------------------------------------------*
2487 *  |  Check if a neutron/proton has been left: if so add to the
2488 *  |  emitted particle list, conserving energyu (but not momentum)
2489       IF ( NINT (ANOW) .EQ. 1 ) THEN
2490          IF ( TV .LT. 0.D+00 ) THEN
2491             LRESMP = .TRUE.
2492             RETURN
2493          END IF
2494          PXRES = PXTTOT - PXINTR - PXNUCR
2495          PYRES = PYTTOT - PYINTR - PYNUCR
2496          PZRES = PZTTOT - PZINTR - PZNUCR
2497          PTRES = SQRT ( PXRES**2 + PYRES**2 + PZRES**2 )
2498          IRN = IRN + 1
2499          ANOW  = ANOW  - 1.D+00
2500          IF ( NINT (ZNOW) .GT. 0.D+00 ) THEN
2501             ITRN (IRN) = 1
2502             ZNOW  = ZNOW  - 1.D+00
2503             KTARP = KTARP + 1
2504             IGREYP = IGREYP + 1
2505             TSTRCK = ETTOT - AM (ITRN(IRN)) - ENUCR - EINTR
2506             EINCP  = EINCP + TSTRCK
2507          ELSE
2508             ITRN (IRN) = 1
2509             KTARN = KTARN + 1
2510             IGREYN = IGREYN + 1
2511             TSTRCK = ETTOT - AM (ITRN(IRN)) - ENUCR - EINTR
2512             EINCN  = EINCN + TSTRCK
2513          END IF
2514          TV  = 0.D+00
2515          ELR (IRN) = TSTRCK + AM (ITRN(IRN))
2516          PLR (IRN) = SQRT ( TSTRCK * ( TSTRCK + 2.D+00 * AM(ITRN(IRN)) )
2517      &             )
2518          PTRES = PLR (IRN) / MAX ( PTRES, ANGLGB )
2519          CXRN (IRN) = PXRES * PTRES
2520          CYRN (IRN) = PYRES * PTRES
2521          CZRN (IRN) = PZRES * PTRES
2522 *  |   updating the conservation counters in common balanc
2523          EINTR  = EINTR + ELR(IRN)
2524          PXINTR = PXINTR + PLR(IRN)*CXRN(IRN)
2525          PYINTR = PYINTR + PLR(IRN)*CYRN(IRN)
2526          PZINTR = PZINTR + PLR(IRN)*CZRN(IRN)
2527          ICINTR = ICINTR + ICH(ITRN(IRN))
2528          IBINTR = IBINTR + IBAR(ITRN(IRN))
2529          IRNTOT = IRNTOT + 1
2530       END IF
2531 *  |
2532 *  +-------------------------------------------------------------------*
2533       RETURN
2534       END