GEANT321/fluka/nucevv.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 NUCEVV.FOR
  13 *COPY NUCEVV
  14 *
  15 *=== nucevv ===========================================================*
  16 *
  17       SUBROUTINE NUCEVV ( NNHAD, KPROJ, PPPROJ, EKPROJ, TXI, TYI, TZI )
  18
  19 #include "geant321/dblprc.inc"
  20 #include "geant321/dimpar.inc"
  21 #include "geant321/iounit.inc"
  22 *----------------------------------------------------------------------*
  23 *     Nucevt90: new version by A. Ferrari INFN-Milan and CERN-SPS
  24 *
  25 *       Besides code cleaning, some changes have been made to extract
  26 *       the quantities needed for a correct energy, momentum, electric
  27 *       charge and baryonic charge conservation. This quantities are put
  28 *       in /balanc/ common
  29 *----------------------------------------------------------------------*
  30 *
  31 C
  32 C**************************************************************
  33 C        /NUCPAR/
  34 C        PARTICLES CREATED IN NUCEVT
  35 C           PXNU,PYNU,PZNU = X-, Y- AND Z-COMPONENTS OF THE
  36 C           LAB MOMENTUM OF THE SECONDARY. COORDINATE SYSTEM
  37 C           DEFINED BY THE PRIMARY PARTICLE.
  38 C           HEPNU    =  TOTAL ENERGY IN THE LAB
  39 C           AMNU     =  REST MASS
  40 C           ICHNU    =  CHARGE
  41 C           IBARNU   =  BARYONIC NUMBER
  42 C           ANNU     =  NAME OF THE PARTICLE
  43 C           NRENU    =  TYPE NUMBER OF THE PARTICLE
  44 C**************************************************************
  45 C
  46 #include "geant321/balanc.inc"
  47 #include "geant321/corinc.inc"
  48 #include "geant321/depnuc.inc"
  49 #include "geant321/hadpar.inc"
  50 #include "geant321/inpdat2.inc"
  51 #include "geant321/nucdat.inc"
  52 #include "geant321/nucpar.inc"
  53 #include "geant321/parevt.inc"
  54 #include "geant321/part.inc"
  55 #include "geant321/resnuc.inc"
  56       COMMON /FKPRIN/ IPRI, INIT
  57 *  The following dimension statement is now obsolete
  58 *     DIMENSION XQT(20), XDQT(20), IFQT(3,20)
  59       REAL RNDM(1)
  60       LOGICAL LDSAVE
  61       DIMENSION PTHRSH (39)
  62       SAVE PTHRSH
  63       DATA PTHRSH / 16*5.D+00,2*2.5D+00,5.D+00,3*2.5D+00,8*5.D+00,
  64      &              9*2.5D+00 /
  65 *
  66 *
  67 *
  68 *----------------------------------------------------------------------*
  69 *         Eproj  = total energy of the projectile
  70 *         Amproj = mass energy of the projectile
  71 *         Ekproj = kinetic energy of the projectile
  72 *         Pproj  = momentum of the projectile
  73 *         Eeproj = original total energy of the projectile
  74 *         Ppproj = original momentum of the projectile
  75 *         Ibproj = barionic charge of the projectile
  76 *         Icproj = electric charge of the projectile
  77 *         Nnhad  = number of particles produced in Nucevt
  78 *         Nevt   = number of high energy collisions
  79 *         Atemp  = local variable with the mass number of the residual
  80 *                  nucleus : it is updated at any interaction and alrea-
  81 *                  dy includes contributions from cascade particles
  82 *                  (even though actually they are produced after the
  83 *                   high energy interactions: since there is no real
  84 *                   correlation except for the one with Nsea, this does
  85 *                   not represent a problem)
  86 *         Ztemp  = same as Atemp for the atomic number
  87 *----------------------------------------------------------------------*
  88 *
  89       TXX = TXI
  90       TYY = TYI
  91       TZZ = TZI
  92       PPROJ = PPPROJ
  93       AMPROJ= AM (KPROJ)
  94       EPROJ = EKPROJ + AMPROJ
  95       EEPROJ= EPROJ
  96 C
  97       IBPROJ= IBAR (KPROJ)
  98       ICPROJ= ICH  (KPROJ)
  99 *  Set Atemp and Ztemp to their initial values
 100       ATEMP = IBTAR  - IGREYP - IGREYN
 101       ZTEMP = ICHTAR - IGREYP
 102       NNHAD = 0
 103       NEVT  = 0
 104 * Now Nsea is sampled inside Corrin and passed through common /corrinc
 105       IF ( NSEA .EQ. 0 ) GO TO 1000
 106       IF ( INIT .EQ. 1 ) WRITE(LUNOUT,1)NNHAD,KPROJ,PPROJ,
 107      *                   AMPROJ,EPROJ,ANUAV,NSEA
 108    1  FORMAT (1X,2I5,4F12.5,I10)
 109 *
 110 *  We have now sampled Nsea, the number of quark-antiquark pairs
 111 *  interacting besides the valence interaction (total interactions
 112 *  = Nsea+1
 113 *
 114 *or      IF (INIT.EQ.1) WRITE(LUNOUT,4)XO,(XSEA(I),XASEA(I),I=1,NSEA)
 115 *or   4  FORMAT (10F12.6)
 116 *  +-------------------------------------------------------------------*
 117 *  |  Sample effective meson nucleus collisions
 118 *  |  Em, im, ekm, pm refer to the meson quantities
 119       DO 6 I = 1, NSEA
 120          EM = EKPROJ * ( XSEA(I) + XASEA(I) )
 121 *  |  +----------------------------------------------------------------*
 122 *  |  |  Selection of the ddbar or qqbar composition
 123          CALL GRNDM(RNDM,1)
 124          IF ( RNDM (1) .LT. 0.5D+00 ) THEN
 125             IM = 23
 126 *  |  |
 127 *  |  +----------------------------------------------------------------*
 128 *  |  |
 129          ELSE
 130             IM = 26
 131          END IF
 132 *  |  |
 133 *  |  +----------------------------------------------------------------*
 134          AMPIO  = AM (IM)
 135          EKM = EM - AMPIO
 136 *  |  +----------------------------------------------------------------*
 137 *  |  |  Energy and/or momentum is too low!!!
 138          IF ( EKM .LT. ETHSEA .OR. PPROJ .LT. PTHRSH (KPTOIP(KPROJ)) )
 139      &      THEN
 140 *  |  |  +-------------------------------------------------------------*
 141 *  |  |  |
 142             IF ( I .LT. NSEA ) THEN
 143                II = I + 1
 144                XSEA(II) = XSEA(II)  + XSEA(I)
 145                XASEA(II)= XASEA(II) + XASEA(I)
 146 *  |  |  |
 147 *  |  |  +-------------------------------------------------------------*
 148 *  |  |  |
 149             ELSE
 150             END IF
 151 *  |  |  |
 152 *  |  |  +-------------------------------------------------------------*
 153             KT = IJTARG (I)
 154 *  |  |  +-------------------------------------------------------------*
 155 *  |  |  | Kt is the index of the target nucleon (1=proton,8=neutron)
 156             IF ( KT .EQ. 1 ) THEN
 157                ZNOW  = ZNOW + 1.D+00
 158                KTARP = KTARP - 1
 159                ZNCOLL = ZNCOLL - 1.D+00
 160 *  |  |  |
 161 *  |  |  +-------------------------------------------------------------*
 162 *  |  |  |
 163             ELSE
 164                KTARN = KTARN - 1
 165             END IF
 166 *  |  |  |
 167 *  |  |  +-------------------------------------------------------------*
 168             ANOW   = ANOW   + 1.D+00
 169             ANCOLL = ANCOLL - 1.D+00
 170             GO TO 6
 171          END IF
 172 *  |  |
 173 *  |  +----------------------------------------------------------------*
 174          AMPIO2 = AMPIO * AMPIO
 175          EPROLD = EKPROJ + AMPROJ
 176          PPROLD = PPROJ
 177 *  | After selection of the meson energy, Ekproj is updated
 178          EKPROJ = EKPROJ - EM
 179          EPROJ  = EKPROJ + AMPROJ
 180 *  |  +----------------------------------------------------------------*
 181 *  |  | All ekproj energy is transferred to em if ekproj < 4 GeV
 182 *  |  | Now modified since it causes troubles to energy conservation by
 183 *  |  | A. Ferrari: force them to have only a valence interaction
 184          IF ( EKPROJ .LT. 3.99D0 ) THEN
 185             EKPROJ = EKPROJ + EM
 186 *  |  |  +-------------------------------------------------------------*
 187 *  |  |  |
 188             DO 666 IN = I, NSEA
 189                KT = IJTARG (IN)
 190 *  |  |  |  +----------------------------------------------------------*
 191 *  |  |  |  | Kt is the index of the target nucleon (1=proton,8=neutron)
 192                IF ( KT .EQ. 1 ) THEN
 193                   ZNOW  = ZNOW + 1.D0
 194                   KTARP = KTARP - 1
 195 *  |  |  |  |
 196 *  |  |  |  +----------------------------------------------------------*
 197 *  |  |  |  |
 198                ELSE
 199                   KTARN = KTARN - 1
 200                END IF
 201 *  |  |  |  |
 202 *  |  |  |  +----------------------------------------------------------*
 203                ANOW  = ANOW + 1.D0
 204  666        CONTINUE
 205 *  |  |  |
 206 *  |  |  +-------------------------------------------------------------*
 207             GO TO 66
 208          END IF
 209 *  |  |
 210 *  |  +----------------------------------------------------------------*
 211 *  |  Now the wounded nucleus is selected inside Corevt
 212          KT    = IJTARG (I)
 213          ATEMP = ATEMP - 1.D+00
 214          IF ( KT .EQ. 1 ) ZTEMP = ZTEMP - 1.D+00
 215 *  |  +---------------------------------------------------------------*
 216 *  |  |
 217          IF ( NINT ( ATEMP ) .EQ. 1 ) THEN
 218             LLASTN = .TRUE.
 219             KTLAST = KT
 220             KTINC  = IJTARG ( NSEA + 1 )
 221             PM = SQRT ( EM**2 - AMPIO2 )
 222             AMLAST = AM (KTLAST)
 223             AMINC  = AM (KTINC)
 224             UMIN2  = ( AMINC + AMLAST )**2
 225             ITJ = MIN ( KTINC, 2 )
 226             DELTAE = V0WELL (ITJ) - EFRMAV (ITJ)
 227 *  |  |  +------------------------------------------------------------*
 228 *  |  |  |   Selection of the invariant mass of the two last nucleon
 229 *  |  |  |   system
 230 2020        CONTINUE
 231                EKPROJ = EKPROJ - DELTAE
 232                EFRM   = EFRM  - DELTAE
 233                ELEFT  = ETTOT - EINTR - EUZ - EKPROJ - AMPROJ - EM
 234                PPROJ  = SQRT ( EKPROJ * ( EKPROJ + 2.D+00 * AMPROJ ) )
 235                PXLEFT = PXTTOT - PXINTR - PUX - ( PPROJ + PM ) * TXX
 236                PYLEFT = PYTTOT - PYINTR - PUY - ( PPROJ + PM ) * TYY
 237                PZLEFT = PZTTOT - PZINTR - PUZ - ( PPROJ + PM ) * TZZ
 238                UMO2 = ELEFT**2 - PXLEFT**2 - PYLEFT**2 - PZLEFT**2
 239 *  |  |  |  +----------------------------------------------------------*
 240 *  |  |  |  |
 241                IF ( UMO2 .LE. UMIN2 ) THEN
 242                   PPDTPL = PXLEFT * TXX + PYLEFT * TYY + PZLEFT * TZZ
 243                   DELTAE = 0.51D+00 * ( UMIN2 - UMO2 ) / ( ELEFT -
 244      &                     PPDTPL * ( EKPROJ + AMPROJ ) / PPROJ )
 245 *  |  |  |  |  +-------------------------------------------------------*
 246 *  |  |  |  |  |   Skip the sea interaction if deltae < 0
 247                   IF ( DELTAE .LT. 0.D+00 ) THEN
 248                      EKPROJ = EKPROJ + EM
 249                      ATEMP = ATEMP + 1.D+00
 250                      IF ( KT .EQ. 1 ) ZTEMP = ZTEMP + 1.D+00
 251 *  |  |  |  |  |  +----------------------------------------------------*
 252 *  |  |  |  |  |  |
 253                      DO 777 IN = I, NSEA
 254                         KT = IJTARG (IN)
 255 *  |  |  |  |  |  |  +-------------------------------------------------*
 256 *  |  |  |  |  |  |  |  Kt is the index of the target nucleon
 257 *  |  |  |  |  |  |  | (1=proton,8=neutron)
 258                         IF ( KT .EQ. 1 ) THEN
 259                            ZNOW  = ZNOW + 1.D0
 260                            KTARP = KTARP - 1
 261 *  |  |  |  |  |  |  |
 262 *  |  |  |  |  |  |  +-------------------------------------------------*
 263 *  |  |  |  |  |  |  |
 264                         ELSE
 265                            KTARN = KTARN - 1
 266                         END IF
 267 *  |  |  |  |  |  |  |
 268 *  |  |  |  |  |  |  +-------------------------------------------------*
 269                         ANOW  = ANOW + 1.D0
 270  777                 CONTINUE
 271 *  |  |  |  |  |  |
 272 *  |  |  |  |  |  +----------------------------------------------------*
 273                      PPROJ  = SQRT ( EKPROJ * ( EKPROJ + 2.D+00
 274      &                      * AMPROJ ) )
 275                      GO TO 66
 276 *  |  |  |  |  |
 277 *  |  |  |  |  +-------------------------------------------------------*
 278 *  |  |  |  |  |
 279                   ELSE IF ( DELTAE .GE. EKPROJ ) THEN
 280                      WRITE ( LUNERR,* )' Nucevv: sea call impossible ',
 281      &                                 ' to get Umin2, go to resampling'
 282                      LRESMP = .TRUE.
 283                      RETURN
 284                   END IF
 285 *  |  |  |  |  |
 286 *  |  |  |  |  +-------------------------------------------------------*
 287                   GO TO 2020
 288 *  |  |  |
 289 *  |  |  +-<|--<--<--<--< We need more invariant mass!!
 290                END IF
 291 *  |  |  |  |
 292 *  |  |  |  +----------------------------------------------------------*
 293 *  |  |  |
 294 *  |  |  +-------------------------------------------------------------*
 295 *  |  |  Now we will divide Eleft and Pleft between the two
 296 *  |  |  left nucleons!
 297             UMO = SQRT ( UMO2 )
 298             ELACMS = 0.5D+00 * ( UMO2 + AMLAST**2 - AMINC**2 )
 299      &             / UMO
 300             EINCMS = UMO - ELACMS
 301             PCMS   = SQRT ( ELACMS**2 - AMLAST**2 )
 302             GAMCM = ELEFT  / UMO
 303             ETAX  = PXLEFT / UMO
 304             ETAY  = PYLEFT / UMO
 305             ETAZ  = PZLEFT / UMO
 306             CALL RACO ( CXXINC, CYYINC, CZZINC )
 307             PCMSX = PCMS * CXXINC
 308             PCMSY = PCMS * CYYINC
 309             PCMSZ = PCMS * CZZINC
 310 *  |  |  Now go back from the CMS frame to the lab frame!!!
 311 *  |  |  First the "inc" nucleon:
 312             ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
 313             EKINC  = GAMCM * EINCMS + ETAPCM - AMINC
 314             PHELP  = ETAPCM / (GAMCM + 1.D+00) + EINCMS
 315             PXXINC = PCMSX + ETAX * PHELP
 316             PYYINC = PCMSY + ETAY * PHELP
 317             PZZINC = PCMSZ + ETAZ * PHELP
 318 *  |  |  Now the "last" nucleon
 319             EKLAST = GAMCM * ELACMS - ETAPCM - AMLAST
 320             PHELP  = - ETAPCM / (GAMCM + 1.D+00) + ELACMS
 321             PXLAST = - PCMSX + ETAX * PHELP
 322             PYLAST = - PCMSY + ETAY * PHELP
 323             PZLAST = - PCMSZ + ETAZ * PHELP
 324             TVEUZ  = 0.D+00
 325             TVGRE0 = 0.D+00
 326             EINCT  = EINCP + EINCN
 327             IF ( EINCT .GT. 0.D+00 ) THEN
 328                EINCP  = EINCP - EINCP / EINCT * TVGREY
 329                EINCN  = EINCN - EINCN / EINCT * TVGREY
 330             END IF
 331             TVGREY = 0.D+00
 332             PSEA = PM + PPROJ - PPROLD
 333             LDSAVE = LDIFFR (KPTOIP(IM))
 334             LDIFFR (KPTOIP(IM)) = .FALSE.
 335 *  |  |  The last argument for ferevv is set to 0 to signal that it
 336 *  |  |  is not a true valence call
 337             IVFLAG = 0
 338             CALL FEREVV ( NHAD, IM, KT, PM, EKM, TXX, TYY, TZZ, ATEMP,
 339      &                    ZTEMP, IVFLAG )
 340             LDIFFR (KPTOIP(IM)) = LDSAVE
 341             IF ( LRESMP ) RETURN
 342 *  |  |
 343 *  |  +---------------------------------------------------------------*
 344 *  |  |
 345          ELSE
 346             TXXOLD = TXX
 347             TYYOLD = TYY
 348             TZZOLD = TZZ
 349             PPROJ  = SQRT ( EKPROJ * ( EKPROJ + 2.D+00 * AMPROJ ) )
 350             IF ( PPROJ .GT. PPROLD ) THEN
 351                PPROJ = PPROLD * EPROJ / EPROLD
 352             END IF
 353             CALL FERSEA ( NHAD, IM, KT, EKM, TXX, TYY, TZZ, ATEMP,
 354      &                    ZTEMP, KPROJ, EPROJ, PPROJ, EPROLD, PPROLD )
 355             IF ( LRESMP ) RETURN
 356          END IF
 357 *  |  |
 358 *  |  +----------------------------------------------------------------*
 359          NEVT = NEVT + 1
 360 *  |  +----------------------------------------------------------------*
 361 *  |  |
 362          IF ( NNHAD + NHAD .GT. MXPNUC ) THEN
 363             WRITE (LUNOUT,1101) NNHAD
 364 1101        FORMAT('  NHAD IN NUCEVT TOO BIG CHANGE DIMENSION',
 365      *             '  NNHAD=',I10)
 366             STOP
 367          END IF
 368 *  |  |
 369 *  |  +----------------------------------------------------------------*
 370
 371 *  |  +----------------------------------------------------------------*
 372 *  |  |          Looping over the produced particles!!!
 373          DO 7 JJ = 1, NHAD
 374             NNHAD = NNHAD + 1
 375             J     = NNHAD
 376             PXNU(J)  = PXH(JJ)
 377             PYNU(J)  = PYH(JJ)
 378             PZNU(J)  = PZH(JJ)
 379             HEPNU(J) = HEPH(JJ)
 380             AMNU(J)  = AMH(JJ)
 381             IBARNU(J)= IBARH(JJ)
 382             ICHNU(J) = ICHH(JJ)
 383             NRENU(J) = NREH(JJ)
 384             ANNU(J)  = ANH(JJ)
 385 *  |  |          Updating energy and momentum accumulators
 386             PUX = PUX+PXNU(J)
 387             PUY = PUY+PYNU(J)
 388             PUZ = PUZ+PZNU(J)
 389             EUZ = EUZ+HEPNU(J)
 390             ICU = ICU+ICHNU(J)
 391             IBU = IBU+IBARNU(J)
 392    7     CONTINUE
 393 *  |  |
 394 *  |  +----------------------------------------------------------------*
 395    6  CONTINUE
 396 *  |
 397 *  +-------------------------------------------------------------------*
 398   66  CONTINUE
 399 *
 400 *  Computing the actual pproj
 401 *
 402       EPROJ = EKPROJ + AMPROJ
 403 1000  CONTINUE
 404 *  +-------------------------------------------------------------------+
 405 *  |  Now modified because of troubles to energy conservation by
 406 *  |  A. Ferrari: if the incident projectile is a proton or a neutron
 407 *  |  maybe it was stopped in the nucleus, else it is added to the stack,
 408 *  |  other particles are also added to the secondary stack.
 409 *  |  May be is not physical, but better than nothing!!! A different
 410 *  |  solution maybe to let meson and other baryons decay at rest or
 411 *  |  also force them to have only a valence interaction or also
 412 *  |  convert a charged meson plus a residual neutron/proton to a
 413 *  |  proton/neutron and add the extra energy to the sea meson and
 414 *  |  so on        ????????????????????
 415 *  |  Now the condition Pproj .ge. .4 is always fulfilled!!!!!
 416       IF (PPROJ .LT. 4.D0) THEN
 417 *  |  +---------------------------------------------------------------*
 418 *  |  |
 419          DKPR30 = KPROJ-30
 420          IF ( PPROJ - 1.5D+00 * SIGN ( ONEONE, DKPR30 )
 421      &        .LT. 4.D0 ) THEN
 422             WRITE ( LUNERR,* )' Nucevt: Pproj < 4 GeV/c!!!', PPROJ,
 423      &                          KPROJ
 424          END IF
 425 *  |  |
 426 *  |  +---------------------------------------------------------------*
 427       END IF
 428 *  |
 429 *  +------------------------------------------------------------------*
 430 *   Now the wounded nucleus is selected inside Corevt
 431       KT = IJTARG ( NSEA + 1 )
 432 *   Update the Nsea value
 433       NSEA = NEVT
 434 *or      IF (INIT.EQ.1) WRITE(LUNOUT,162)NEVT,NHAD,IM,KPROJ,KT,PM,
 435 *or     &EM,EKM,PM,EKPROJ,PPROJ,EPROJ
 436 *
 437 *  Now ferevt performs the interaction: no longer a meson but the actual
 438 *  projectile
 439 *
 440       ATEMP = ATEMP - 1.D+00
 441       IF ( KT .EQ. 1 ) ZTEMP = ZTEMP - 1.D+00
 442 *  +-------------------------------------------------------------------*
 443 *  |
 444       IF ( NINT ( ATEMP ) .EQ. 1 ) THEN
 445          LLASTN = .TRUE.
 446          KTLAST = KT
 447          IF ( ZNOW .GT. 0.D+00 ) THEN
 448             KTINC  = 1
 449          ELSE
 450             KTINC  = 8
 451          END IF
 452          AMLAST = AM (KTLAST)
 453          AMINC  = AM (KTINC)
 454          UMIN2  = ( AMINC + AMLAST )**2
 455          ITJ = MIN ( KTINC, 2 )
 456 *  |  +----------------------------------------------------------------*
 457 *  |  |   Selection of the invariant mass of the two last nucleon
 458 *  |  |   system
 459 2040     CONTINUE
 460             ELEFT  = ETTOT - EINTR - EUZ - EKPROJ - AMPROJ
 461             PXLEFT = PXTTOT - PXINTR - PUX - PPROJ * TXX
 462             PYLEFT = PYTTOT - PYINTR - PUY - PPROJ * TYY
 463             PZLEFT = PZTTOT - PZINTR - PUZ - PPROJ * TZZ
 464             UMO2 = ELEFT**2 - PXLEFT**2 - PYLEFT**2 - PZLEFT**2
 465 *  |  |  +-------------------------------------------------------------*
 466 *  |  |  |
 467             IF ( UMO2 .LE. UMIN2 ) THEN
 468                PPDTPL = PXLEFT * TXX + PYLEFT * TYY + PZLEFT * TZZ
 469                DELTAE = 0.51D+00 * ( UMIN2 - UMO2 ) / ( ELEFT -
 470      &                  PPDTPL * ( EKPROJ + AMPROJ ) / PPROJ )
 471 *  |  |  |  +----------------------------------------------------------*
 472 *  |  |  |  |
 473                IF ( DELTAE .GE. EKPROJ .OR. DELTAE .LT. 0.D+00 ) THEN
 474                   WRITE ( LUNERR,* )' Nucevv: valence call impossible ',
 475      &                              ' to get Umin2, go to resampling'
 476                   LRESMP = .TRUE.
 477                   RETURN
 478                END IF
 479 *  |  |  |  |
 480 *  |  |  |  +----------------------------------------------------------*
 481                EKPROJ = EKPROJ - DELTAE
 482                PLA    = PPROJ
 483                PPROJ  = SQRT ( EKPROJ * ( EKPROJ + 2.D+00 * AMPROJ ) )
 484                PLA    = PPROJ - PLA
 485                EFRM   = EFRM  - DELTAE
 486                PXFRM  = PXFRM + TXX * PLA
 487                PYFRM  = PYFRM + TYY * PLA
 488                PZFRM  = PZFRM + TZZ * PLA
 489                GO TO 2040
 490 *  |  |
 491 *  |  +-<|--<--<--<--< We need more invariant mass!!
 492             END IF
 493 *  |  |  |
 494 *  |  |  +-------------------------------------------------------------*
 495 *  |  |
 496 *  |  +----------------------------------------------------------------*
 497 *  |  Now we will divide Eleft and Pleft between the two
 498 *  |  left nucleons!
 499          UMO = SQRT ( UMO2 )
 500          ELACMS = 0.5D+00 * ( UMO2 + AMLAST**2 - AMINC**2 ) / UMO
 501          EINCMS = UMO - ELACMS
 502          PCMS   = SQRT ( ELACMS**2 - AMLAST**2 )
 503          GAMCM = ELEFT  / UMO
 504          ETAX  = PXLEFT / UMO
 505          ETAY  = PYLEFT / UMO
 506          ETAZ  = PZLEFT / UMO
 507          CALL RACO ( CXXINC, CYYINC, CZZINC )
 508          PCMSX = PCMS * CXXINC
 509          PCMSY = PCMS * CYYINC
 510          PCMSZ = PCMS * CZZINC
 511 *  |  Now go back from the CMS frame to the lab frame!!!
 512 *  |  First the "inc" nucleon:
 513          ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
 514          EKINC  = GAMCM * EINCMS + ETAPCM - AMINC
 515          PHELP  = ETAPCM / (GAMCM + 1.D+00) + EINCMS
 516          PXXINC = PCMSX + ETAX * PHELP
 517          PYYINC = PCMSY + ETAY * PHELP
 518          PZZINC = PCMSZ + ETAZ * PHELP
 519 *  |  Now the "last" nucleon
 520          EKLAST = GAMCM * ELACMS - ETAPCM - AMLAST
 521          PHELP  = - ETAPCM / (GAMCM + 1.D+00) + ELACMS
 522          PXLAST = - PCMSX + ETAX * PHELP
 523          PYLAST = - PCMSY + ETAY * PHELP
 524          PZLAST = - PCMSZ + ETAZ * PHELP
 525          TVEUZ  = 0.D+00
 526          TVGRE0 = 0.D+00
 527          EINCT  = EINCP + EINCN
 528          IF ( EINCT .GT. 0.D+00 ) THEN
 529             EINCP  = EINCP - EINCP / EINCT * TVGREY
 530             EINCN  = EINCN - EINCN / EINCT * TVGREY
 531          END IF
 532          TVGREY = 0.D+00
 533       END IF
 534 *  |
 535 *  +-------------------------------------------------------------------*
 536 *  The last argument for ferevv is set to 1 to signal that it
 537 *  is a true valence call
 538       IVFLAG = 1
 539       CALL FEREVV ( NHAD, KPROJ, KT, PPROJ, EKPROJ, TXX, TYY, TZZ,
 540      &              ATEMP, ZTEMP, IVFLAG )
 541       IF ( LRESMP ) RETURN
 542       NEVT = NEVT + 1
 543 *  +-------------------------------------------------------------------*
 544 *  |
 545       IF ( NNHAD + NHAD .GT. MXPNUC ) THEN
 546          WRITE (LUNOUT,1101) NNHAD
 547          STOP
 548       END IF
 549 *  |
 550 *  +-------------------------------------------------------------------*
 551
 552 *  +-------------------------------------------------------------------*
 553 *  |       Looping over the produced particles
 554       DO 8 JJ=1,NHAD
 555          NNHAD = NNHAD+1
 556          J = NNHAD
 557          PXNU(J)  = PXH(JJ)
 558          PYNU(J)  = PYH(JJ)
 559          PZNU(J)  = PZH(JJ)
 560          HEPNU(J) = HEPH(JJ)
 561          AMNU(J)  = AMH(JJ)
 562          IBARNU(J)= IBARH(JJ)
 563          ICHNU(J) = ICHH(JJ)
 564          NRENU(J) = NREH(JJ)
 565          ANNU(J)  = ANH(JJ)
 566 *  |             Updating energy and momentum accumulators
 567          PUX = PUX+PXNU(J)
 568          PUY = PUY+PYNU(J)
 569          PUZ = PUZ+PZNU(J)
 570          EUZ = EUZ+HEPNU(J)
 571          ICU = ICU+ICHNU(J)
 572          IBU = IBU+IBARNU(J)
 573    8  CONTINUE
 574 *  |
 575 *  +-------------------------------------------------------------------*
 576  888  CONTINUE
 577 *     IF (PPROJ .LT. 4.D0) INIT=1
 578 **** print and test energy conservation
 579 *
 580       ETOT  = EEPROJ + KTARP * ( AMNUCL (1) - EBNDNG (1) )
 581      &      + KTARN * ( AMNUCL (2) - EBNDNG (2) ) + EFRM
 582       ICTOT = ICH  (KPROJ) + KTARP
 583       IBTOT = IBAR (KPROJ) + KTARP + KTARN
 584       EMIN  = 1.D-10 * ETOT
 585       PMIN  = 1.D-10 * PPPROJ
 586 *  +-------------------------------------------------------------------*
 587 *  |
 588       IF (ICTOT .NE. ICU .OR. IBTOT .NE. IBU) THEN
 589          LRESMP = .TRUE.
 590          WRITE(LUNERR,*)' NUCEVT-FEREVT CHARGE CONSERVATION FAILURE:',
 591      &                  ' ICU,ICTOT,IBU,IBTOT',ICU,ICTOT,IBU,IBTOT
 592       END IF
 593 *  |
 594 *  +-------------------------------------------------------------------*
 595       PXCHCK = PPPROJ * TXI + PSEA * TXX
 596       PYCHCK = PPPROJ * TYI + PSEA * TYY
 597       PZCHCK = PPPROJ * TZI + PSEA * TZZ
 598 *  +-------------------------------------------------------------------*
 599 *  |
 600       IF ( ABS ( PUX - PXFRM - PXCHCK ) .GT. PMIN ) THEN
 601          LRESMP = .TRUE.
 602          WRITE(LUNERR,*)' NUCEVT-FEREVT PX CONSERVATION FAILURE:',
 603      &                  ' PUX,PXFRM+PXCHCK',PUX,PXFRM+PXCHCK
 604       END IF
 605 *  |
 606 *  +-------------------------------------------------------------------*
 607 *  +-------------------------------------------------------------------*
 608 *  |
 609       IF ( ABS ( PUY - PYFRM - PYCHCK ) .GT. PMIN ) THEN
 610          WRITE(LUNERR,*)' NUCEVT-FEREVT PY CONSERVATION FAILURE:',
 611      &                  ' PUY,PYFRM+PYCHCK*TYY',PUY,PYFRM+PYCHCK
 612          LRESMP = .TRUE.
 613       END IF
 614 *  |
 615 *  +-------------------------------------------------------------------*
 616 *  +-------------------------------------------------------------------*
 617 *  |
 618       IF ( ABS ( PUZ - PZFRM - PZCHCK ) .GT. PMIN ) THEN
 619          LRESMP = .TRUE.
 620          WRITE(LUNERR,*)' NUCEVT-FEREVT PZ CONSERVATION FAILURE:',
 621      &                  ' PUZ,PZFRM+PZCHCK',PUZ,PZFRM+PZCHCK
 622       END IF
 623 *  |
 624 *  +-------------------------------------------------------------------*
 625 *  +-------------------------------------------------------------------*
 626 *  |
 627       IF ( ABS ( ETOT - EUZ ) .GT. EMIN ) THEN
 628          LRESMP = .TRUE.
 629          WRITE(LUNERR,*)' NUCEVT-FEREVT E CONSERVATION FAILURE:',
 630      &                  ' EUZ,ETOT',EUZ,ETOT
 631       END IF
 632 *  |
 633 *  +-------------------------------------------------------------------*
 634       EUZ0 = EUZ-NEVT*AM(1)
 635       IF (INIT.NE.1) GO TO 11
 636       WRITE(LUNOUT,12)NNHAD,KPROJ,PPPROJ,EPROJ,PUX,PUY,PUZ,EUZ0,
 637      *ICU,IBU,NNHAD
 638   12  FORMAT (1X,2I5,6F12.6,3I5)
 639 *  +-------------------------------------------------------------------*
 640 *  |
 641       DO 13 I=1,NNHAD
 642          WRITE(LUNOUT,14)I,NRENU(I),ICHNU(I),IBARNU(I),ANNU(I),PXNU(I),
 643      *                   PYNU(I),PZNU(I),HEPNU(I),AMNU(I)
 644   14     FORMAT (1X,4I5,A8,6F12.6)
 645   13  CONTINUE
 646 *  |
 647 *  +-------------------------------------------------------------------*
 648       INIT=0
 649   11  CONTINUE
 650       RETURN
 651       END