GEANT321/fluka/eventv.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:19:55  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.42  by  S.Giani
  11 *-- Author :
  12 *$ CREATE EVENTV.FOR
  13 *COPY EVENTV
  14 *
  15 *=== eventv ===========================================================*
  16 *
  17       SUBROUTINE EVENTV ( IJIJ, POO, EKE, TXI, TYI, TZI, WE, IMAT )
  18
  19 #include "geant321/dblprc.inc"
  20 #include "geant321/dimpar.inc"
  21 #include "geant321/iounit.inc"
  22 *
  23 *----------------------------------------------------------------------*
  24 *                                                                      *
  25 *  Created  on  7 september 1989      by         Alfredo Ferrari       *
  26 *                                                 INFN - Milan         *
  27 *                                                                      *
  28 *  Last change  on  15 april 1993     by         Alfredo Ferrari       *
  29 *                                                                      *
  30 *  This routine is a strongly modified and improved version of the     *
  31 *  original Eventq routine of the Fluka package: it requires modified  *
  32 *  versions of many other routines                                     *
  33 *                                                                      *
  34 *  Eventv90: new version by A. Ferrari, INFN-Milan and CERN-SPS, 7-9-89*
  35 *       main modifications:                                            *
  36 *                  - new treatment for hydrogen (following Moehring's  *
  37 *                                                suggestion)           *
  38 *                  - new commons and coding for energy, momentum,      *
  39 *                    electric and baryonic charge conservation         *
  40 *                  - new correlations both in the "low" (p < 5 GeV/c)  *
  41 *                    and in the "high" energy models                   *
  42 *                  - complete revision of kinematics for both models   *
  43 *                  - tentative treatment for Lambdas, Sigmas with the  *
  44 *                    "high" energy model down to 2.5 GeV/c to overcome *
  45 *                    the impossibility of Hadrin to treat these parti- *
  46 *                    cles                                              *
  47 *                  - tentative introduction of Xsi0, Xsi-, Omega and   *
  48 *                    their antiparticles and of Sigmabars's            *
  49 *                  - complete new treatment of cascade nucleons        *
  50 *                  - accurate treatment of binding energy effects,     *
  51 *                    using atomic mass tables and the proper isotopic  *
  52 *                    composition of elements                           *
  53 *                  - possibility to use an evaporation module (derived *
  54 *                    from the Hetc-Kfa one) after the interaction      *
  55 *                  - possibility to use a nuclear gamma deexcitation   *
  56 *                    module (written by P.Sala, INFN-Milan) to produce *
  57 *                    deexcitation gammas after the evaporation step    *
  58 *                                                                      *
  59 *----------------------------------------------------------------------*
  60 *
  61       PARAMETER ( COS45 = 0.7071067811865475D+00 )
  62       PARAMETER ( SIN30 = 0.5D+00 )
  63       PARAMETER ( COS30 = 0.8660254037844387D+00 )
  64       PARAMETER ( CSNPRN = 50.D+00 * CSNNRM )
  65 #include "geant321/balanc.inc"
  66 #include "geant321/corinc.inc"
  67 #include "geant321/depnuc.inc"
  68 #include "geant321/fheavy.inc"
  69 #include "geant321/finlsp3.inc"
  70 #include "geant321/finuc.inc"
  71 #include "geant321/hadflg.inc"
  72 #include "geant321/hadpar.inc"
  73 #include "geant321/isotop.inc"
  74 #include "geant321/mapa.inc"
  75 #include "geant321/nucdat.inc"
  76 #include "geant321/nucpar.inc"
  77 #include "geant321/part.inc"
  78 #include "geant321/parevt.inc"
  79 #include "geant321/resnuc.inc"
  80       COMMON / FKCASF / PKFRMI, COSTH, PKIN
  81       COMMON /FKEVNT/ LNUCRI, LHADRI
  82       LOGICAL LNUCRI, LHADRI, LINCCK, LACCEP, LCHTYP
  83       REAL RNDM(1)
  84 *  Note Pthrsh is also in Ferevv!!
  85       DIMENSION SOPPP(2),EXSOP(2),PTHRSH(39),IJNUCR(39)
  86 * Use this "save" if common dbgdbg is commented
  87       SAVE PTHRSH, IPRI, IEVT, IJNUCR
  88       DATA PTHRSH / 16*5.D+00,2*2.5D+00,5.D+00,3*2.5D+00,8*5.D+00,
  89      &              9*2.5D+00 /
  90       DATA IJNUCR / 16*1,2*0,1,3*0,8*1,9*0 /
  91       DATA IPRI / 0 /
  92       DATA IEVT / 0 /
  93 *
  94 *  +-------------------------------------------------------------------*
  95 *  |  Heavy ions are treated in eventa: now switched off!!
  96       IF (IJIJ .EQ. 30) THEN
  97          STOP 'HEAVY ION'
  98       END IF
  99 *  |
 100 *  +-------------------------------------------------------------------*
 101 *   First of all: get the "part" index of the incoming "paprop" ijij
 102       IJ = IPTOKP (IJIJ)
 103       IEVT = IEVT + 1
 104       AMCHCK = 0.5D+00 * ( POO * POO - EKE * EKE ) / EKE
 105       AHELP  =  ABS ( AMCHCK - AM ( IJ ) )
 106 *  +-------------------------------------------------------------------*
 107 *  |
 108       IF ( AHELP .GT. ANGLGB * AM ( IJ ) ) THEN
 109          POO = SQRT ( EKE * ( EKE + 2.D+00 * AM (IJ) ) )
 110       END IF
 111 *  |
 112 *  +-------------------------------------------------------------------*
 113 *
 114 *  Variable initialization for conservation checks
 115 *
 116       PTTOT  = POO
 117       PXTTOT = PTTOT*TXI
 118       PYTTOT = PTTOT*TYI
 119       PZTTOT = PTTOT*TZI
 120       ICHTAR = NINT(ZTAR(IMAT))
 121       ICHTOT = ICHTAR + ICH(IJ)
 122 *  +-------------------------------------------------------------------*
 123 *  |    Choice of the mass number of the target nucleus: use the input
 124 *  |    one if any
 125       IF ( MSSNUM (IMAT) .GT. 0 ) THEN
 126          IBTAR = MSSNUM (IMAT)
 127 *  |
 128 *  +-------------------------------------------------------------------*
 129 *  |    Choice of the mass number of the target nucleus according
 130 *  |    to the natural isotopic composition of element ichtar
 131       ELSE
 132          CALL GRNDM(RNDM,1)
 133          RNDISO = RNDM (1)
 134 *  |  +----------------------------------------------------------------*
 135 *  |  |    Loop on the stable isotopes
 136          DO 25 IS = ISONDX (1,ICHTAR), ISONDX (2,ICHTAR) - 1
 137             RNDISO = RNDISO - ABUISO (IS)
 138             IF ( RNDISO .LT. 0.D+00 ) GO TO 30
 139    25    CONTINUE
 140 *  |  |
 141 *  |  +----------------------------------------------------------------*
 142          IS = ISONDX (2,ICHTAR)
 143    30    CONTINUE
 144          IBTAR = ISOMNM (IS)
 145       END IF
 146 *  |
 147 *  +-------------------------------------------------------------------*
 148       IBTOT  = IBTAR + IBAR(IJ)
 149       BBTAR  = IBTAR
 150       ZZTAR  = ICHTAR
 151       ZTO103 = ZZTAR**0.3333333333333333D+00
 152 *  +-------------------------------------------------------------------*
 153 *  |
 154       IF ( IBTAR .EQ. 1 ) THEN
 155          ITTA   = 8 - 7 * ICHTAR
 156          ETTOT  = AM (ITTA) + EKE + AM(IJ)
 157          AMNTAR = AM (ITTA)
 158 *  |
 159 *  +-------------------------------------------------------------------*
 160 *  | The following should be done with the proper mass of the nuclide
 161       ELSE
 162 *        AMNTAR = BBTAR * AMUC12
 163          AMMTAR = BBTAR * AMUAMU + 0.001D+00 * FKENER ( BBTAR, ZZTAR )
 164          AMNTAR = AMMTAR - ZZTAR * AMELEC + ELBNDE ( ICHTAR )
 165          ETTOT  = AMNTAR + EKE + AM (IJ)
 166       END IF
 167 *  |
 168 *  +-------------------------------------------------------------------*
 169 *  +-------------------------------------------------------------------*
 170 *  |  Kaon-short and kaon-long are transformed into a kaon0 or a
 171 *  |  kaon0-bar
 172       IF (IJ .EQ. 12 .OR. IJ. EQ. 19) THEN
 173          IJ = 24
 174          CALL GRNDM(RNDM,1)
 175          IF (RNDM(1) .GT. 0.5D0) IJ = 25
 176 *  |
 177 *  +-------------------------------------------------------------------*
 178 *  |  Pi0 quark configurations are selected according to 50% uubar and
 179 *  |  50% ddbar
 180       ELSE IF ( IJ .EQ. 23 .OR. IJ .EQ. 26 ) THEN
 181 *  |  +----------------------------------------------------------------*
 182 *  |  |
 183          CALL GRNDM(RNDM,1)
 184          IF ( RNDM (1) .LT. 0.5D+00 ) THEN
 185             IJ = 23
 186 *  |  |
 187 *  |  +----------------------------------------------------------------*
 188 *  |  |
 189          ELSE
 190             IJ = 26
 191          END IF
 192 *  |  |
 193 *  |  +----------------------------------------------------------------*
 194       END IF
 195 *  |
 196 *  +-------------------------------------------------------------------*
 197       ETEPS  = MAX ( 1.D-10 * EKE, 1.D-10 )
 198   50  CONTINUE
 199 *-->-->-->-->--> Come here for resampling
 200 * Reset all the accumulators
 201       NP     = NP0
 202       NPHEAV = 0
 203       LRNFSS = .FALSE.
 204       LEVDIF = .FALSE.
 205       LHADRI = .FALSE.
 206       ICHTOT = ICHTAR + ICH(IJ)
 207       TV = 0.D+00
 208       TVRECL = 0.D+00
 209       TVHEAV = 0.D+00
 210       TVCMS  = 0.D+00
 211       EOTEST = ETTOT
 212       LRESMP = .FALSE.
 213       LLASTN = .FALSE.
 214       LLAST1 = .FALSE.
 215       EUZ = 0.D+00
 216       PUX = 0.D+00
 217       PUY = 0.D+00
 218       PUZ = 0.D+00
 219       TVEUZ = 0.D+00
 220       ICU = 0
 221       IBU = 0
 222       EINTR  = 0.D+00
 223       PXINTR = 0.D+00
 224       PYINTR = 0.D+00
 225       PZINTR = 0.D+00
 226       ICINTR = 0
 227       IBINTR = 0
 228       ENUCR  = 0.D+00
 229       PXNUCR = 0.D+00
 230       PYNUCR = 0.D+00
 231       PZNUCR = 0.D+00
 232       ICNUCR = 0
 233       IBNUCR = 0
 234       EFRM  = 0.D+00
 235       PXFRM = 0.D+00
 236       PYFRM = 0.D+00
 237       PZFRM = 0.D+00
 238       PSEA  = 0.D+00
 239       TVGRE0  = 0.D+00
 240       TVGREY  = 0.D+00
 241       IGREYP  = 0
 242       IGREYN  = 0
 243       KTARP = 0
 244       KTARN = 0
 245       IEVAPL = 0
 246       IEVAPH = 0
 247       IEVPRO = 0
 248       IEVNEU = 0
 249       IEVDEU = 0
 250       IEVTRI = 0
 251       IEV3HE = 0
 252       IEV4HE = 0
 253       IDEEXG = 0
 254       ANOW   = BBTAR
 255       ZNOW   = ZZTAR
 256       DSOPP  = 0.D+00
 257       NNHAD  = 0
 258       RN1GSC = -1.D+00
 259       RN2GSC = -1.D+00
 260       IF ( POO .GE. PTHRSH (IJIJ) ) GO TO 200
 261 * Below 5 GeV/c use Nucrin:
 262  100  CONTINUE
 263       ANCOLL = 1.D+00
 264 *  +-------------------------------------------------------------------*
 265 *  |  Check for Hydrogen
 266       IF ( IBTAR .NE. 1  ) THEN
 267 *  |  +----------------------------------------------------------------*
 268 *  |  |  Steering for Peanut: very temporary
 269 *  |  |  Protons and Neutrons below 260 MeV:
 270          IF ( EKE .LT. 0.260D+00 .AND. ( IJ .EQ. 1 .OR. IJ .EQ. 8 )
 271      &        .AND. IBTAR .GT. 4 ) THEN
 272             CALL PEANUT ( IJ, EKE, POO, TXI, TYI, TZI, WE )
 273             IF ( LRESMP ) GO TO 50
 274             RETURN
 275 *  |  |
 276 *  |  +----------------------------------------------------------------*
 277 *  |  |  Stopping pi-s:
 278          ELSE IF ( IJ .EQ. 14 .AND. EKE .LT. 2.D+00 * GAMMIN .AND. IBTAR
 279      &             .GT. 4 ) THEN
 280             CALL PEANUT ( IJ, EKE, POO, TXI, TYI, TZI, WE )
 281             IF ( LRESMP ) GO TO 50
 282             RETURN
 283          END IF
 284 *  |  |
 285 *  |  +----------------------------------------------------------------*
 286 *  |  Now at first we sample the number of interactions
 287          CALL CORRIN ( ZZTAR, BBTAR, IJ, POO, EKE )
 288       END IF
 289       LNUCRI = .TRUE.
 290       NP = NP0
 291       ELAB  = EKE + AM(IJ)
 292 *  Redefinition of particle types: the recovery of proper charge and
 293 *  strangeness conservation after the interaction is not yet implemented
 294       KPROJ = IJ
 295       LCHTYP = .FALSE.
 296 *  +-------------------------------------------------------------------*
 297 *  |
 298       IF ( KPROJ .LE. 30 ) THEN
 299          GO TO ( 2017, 2018, 2999, 2020, 2021, 2022, 2999, 2999, 2999,
 300      &           2026 ), KPROJ - 16
 301          GO TO 2999
 302 *  |
 303 *  +-------------------------------------------------------------------*
 304 *  |
 305       ELSE
 306          GO TO ( 2031, 2032, 2033, 2034, 2035, 2036, 2037, 2038, 2039 ),
 307      &      KPTOIP (KPROJ) - 30
 308          GO TO 2999
 309       END IF
 310 *  |
 311 *  +-------------------------------------------------------------------*
 312 *  +-------------------------------------------------------------------*
 313 *  |  Lambdas are considered as neutrons: after the interaction
 314 *  |  a possible way to recover strangeness and charge conservation
 315 *  |  is to flip a pi- --> K-, or an eventual neutron into
 316 *  |  a Sigma0/Lambda,
 317 *  |  (we must flip a d quark into an s quark), or a pi0 --> K0bar,
 318 *  |  or a p --> Sigma+
 319  2017 CONTINUE
 320          LCHTYP = .TRUE.
 321          KPROJ = 8
 322 *  |   Adjusting mass and momentum for lambdas (now n)
 323          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 324          PTTOT  = POO
 325          PXTTOT = PTTOT*TXI
 326          PYTTOT = PTTOT*TYI
 327          PZTTOT = PTTOT*TZI
 328       GO TO 2999
 329 *  |
 330 *  +-------------------------------------------------------------------*
 331 *  |  Alambdas are considered as antineutrons : after the interaction
 332 *  |  a possible way to recover strangeness and charge conservation
 333 *  |  is to flip a pi+ --> K+, or an eventual neutronbar into a
 334 *  |  Lambdabar, (we must flip a dbar quark into an sbar quark),
 335 *  |  or a pi0 --> K0, or a pbar --> Sigma+bar (ASigma-) ( there
 336 *  |  are also possible ways through eta, rho and omega mesons )
 337  2018 CONTINUE
 338          LCHTYP = .TRUE.
 339          KPROJ = 9
 340 *  |   Adjusting mass and momentum for alambdas (now an)
 341          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 342          PTTOT  = POO
 343          PXTTOT = PTTOT*TXI
 344          PYTTOT = PTTOT*TYI
 345          PZTTOT = PTTOT*TZI
 346       GO TO 2999
 347 *  |
 348 *  +-------------------------------------------------------------------*
 349 *  |  Sigma-s are considered as neutrons : after the interaction
 350 *  |  a possible way to recover strangeness and charge conservation
 351 *  |  is to flip a pi+ --> K0bar, or an eventual neutron into a Sigma-,
 352 *  |  (we must flip a u quark into an s quark), or a pi0 --> K-,
 353 *  |  or a p --> Sigma0/Lambda
 354  2020 CONTINUE
 355          LCHTYP = .TRUE.
 356          KPROJ = 8
 357 *  |   The following card must be commented following the strangeness
 358 *  |   and charge conservation recovery implementation
 359          ICHTOT = ICHTOT + 1
 360 *  |   Adjusting mass and momentum for Sigma-s (now n)
 361          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 362          PTTOT  = POO
 363          PXTTOT = PTTOT*TXI
 364          PYTTOT = PTTOT*TYI
 365          PZTTOT = PTTOT*TZI
 366       GO TO 2999
 367 *  |
 368 *  +-------------------------------------------------------------------*
 369 *  |  Sigma+s are considered as protons : after the interaction
 370 *  |  a possible way to recover strangeness and charge conservation
 371 *  |  is to flip a pi- --> K-, or an eventual proton into a Sigma+,
 372 *  |  (we must flip a d quark into an s quark), or a pi0 --> K0bar,
 373 *  |  or a n --> Sigma0/Lambda
 374  2021 CONTINUE
 375          LCHTYP = .TRUE.
 376          KPROJ = 1
 377 *  |   Adjusting mass and momentum for Sigma+s (now p)
 378          PTTOT  = POO
 379          PXTTOT = PTTOT*TXI
 380          PYTTOT = PTTOT*TYI
 381          PZTTOT = PTTOT*TZI
 382       GO TO 2999
 383 *  |
 384 *  +-------------------------------------------------------------------*
 385 *  |  Sigma0s are considered as neutrons: after the interaction
 386 *  |  a possible way to recover strangeness and charge conservation
 387 *  |  is to flip a pi- --> K-, or an eventual neutron into
 388 *  |  a Sigma0/Lambda,
 389 *  |  (we must flip a d quark into an s quark), or a pi0 --> K0bar,
 390 *  |  or a p --> Sigma+
 391  2022 CONTINUE
 392          LCHTYP = .TRUE.
 393          KPROJ = 8
 394 *  |   Adjusting mass and momentum for Sigma0s (now n)
 395          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 396          PTTOT  = POO
 397          PXTTOT = PTTOT*TXI
 398          PYTTOT = PTTOT*TYI
 399          PZTTOT = PTTOT*TZI
 400       GO TO 2999
 401 *  |
 402 *  +-------------------------------------------------------------------*
 403 *  |  Pi0 must be 23 for Nucrin
 404  2026 CONTINUE
 405          KPROJ = 23
 406       GO TO 2999
 407 *  |
 408 *  +-------------------------------------------------------------------*
 409 *  |  ASigma-s are considered as pbars: after the interaction
 410 *  |  a possible way to recover strangeness and charge conservation
 411 *  |  is to flip a pi+ --> K+, or an eventual pbar into a ASigma-,
 412 *  |  (we must flip a dbar quark into an sbar quark), or a pi0 --> K0,
 413 *  |  or a nbar --> ASigma0/ALambda
 414  2031 CONTINUE
 415          LCHTYP = .TRUE.
 416          KPROJ = 2
 417 *  |   Adjusting mass and momentum for ASigma-s (now pbar)
 418          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 419          PTTOT  = POO
 420          PXTTOT = PTTOT*TXI
 421          PYTTOT = PTTOT*TYI
 422          PZTTOT = PTTOT*TZI
 423       GO TO 2999
 424 *  |
 425 *  +-------------------------------------------------------------------*
 426 *  |  ASigma0s are considered as nbar: after the interaction
 427 *  |  a possible way to recover strangeness and charge conservation
 428 *  |  is to flip a pi+ --> K+, or an eventual nbar into
 429 *  |  a ASigma0/ALambda,
 430 *  |  (we must flip a dbar quark into an sbar quark), or a pi0 --> K0,
 431 *  |  or a pbar --> ASigma-
 432  2032 CONTINUE
 433          LCHTYP = .TRUE.
 434          KPROJ = 9
 435 *  |   Adjusting mass and momentum for ASigma0s (now nbar)
 436          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 437          PTTOT  = POO
 438          PXTTOT = PTTOT*TXI
 439          PYTTOT = PTTOT*TYI
 440          PZTTOT = PTTOT*TZI
 441       GO TO 2999
 442 *  |
 443 *  +-------------------------------------------------------------------*
 444 *  |  ASigma+s are considered as nbar : after the interaction
 445 *  |  a possible way to recover strangeness and charge conservation
 446 *  |  is to flip a pi- --> K0, or an eventual nbar into a ASigma+,
 447 *  |  (we must flip a ubar quark into an sbar quark), or a pi0 --> K+,
 448 *  |  or a pbar --> ASigma0/ALambda
 449  2033 CONTINUE
 450          LCHTYP = .TRUE.
 451          KPROJ = 9
 452 *  |   The following card must be commented following the strangeness
 453 *  |   and charge conservation recovery implementation
 454          ICHTOT = ICHTOT - 1
 455 *  |   Adjusting mass and momentum for ASigma+s (now nbar)
 456          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 457          PTTOT  = POO
 458          PXTTOT = PTTOT*TXI
 459          PYTTOT = PTTOT*TYI
 460          PZTTOT = PTTOT*TZI
 461       GO TO 2999
 462 *  |
 463 *  +-------------------------------------------------------------------*
 464 *  |  Xsi0s are considered as neutrons: after the interaction
 465 *  |  a possible way to recover strangeness and charge conservation
 466 *  |  is to flip an eventual neutron into a Xsi0,
 467 *  |  or a proton --> Sigma+ and a pi- --> K- (or a pi0 --> K0bar) etc.
 468 *  |  (we must flip two d quarks into s quarks)
 469  2034 CONTINUE
 470          LCHTYP = .TRUE.
 471          KPROJ = 8
 472 *  |   Adjusting mass and momentum for Xsi0s (now n)
 473          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 474          PTTOT  = POO
 475          PXTTOT = PTTOT*TXI
 476          PYTTOT = PTTOT*TYI
 477          PZTTOT = PTTOT*TZI
 478       GO TO 2999
 479 *  |
 480 *  +-------------------------------------------------------------------*
 481 *  |  AXsi0s are considered as nbars: after the interaction
 482 *  |  a possible way to recover strangeness and charge conservation
 483 *  |  is to flip an eventual nbar into a AXsi0,
 484 *  |  or a pbar --> ASigma- and a pi+ --> K+ (or a pi0 --> K0) etc.
 485 *  |  (we must flip two dbar quarks into sbar quarks)
 486  2035 CONTINUE
 487          LCHTYP = .TRUE.
 488          KPROJ = 9
 489 *  |   Adjusting mass and momentum for AXsi0s (now nbar)
 490          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 491          PTTOT  = POO
 492          PXTTOT = PTTOT*TXI
 493          PYTTOT = PTTOT*TYI
 494          PZTTOT = PTTOT*TZI
 495       GO TO 2999
 496 *  |
 497 *  +-------------------------------------------------------------------*
 498 *  |  Xsi-s are considered as protons: after the interaction
 499 *  |  a possible way to recover strangeness and charge conservation
 500 *  |  is to flip an eventual proton into a Xsi-,
 501 *  |  or a neutron --> Sigma- and a pi+ --> K0bar (or a pi0 --> K-) etc.
 502 *  |  (we must flip two u quarks into s quarks)
 503  2036 CONTINUE
 504          LCHTYP = .TRUE.
 505          KPROJ = 1
 506 *  |   The following card must be commented following the strangeness
 507 *  |   and charge conservation recovery implementation
 508          ICHTOT = ICHTOT + 2
 509 *  |   Adjusting mass and momentum for Xsi-s (now p)
 510          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 511          PTTOT  = POO
 512          PXTTOT = PTTOT*TXI
 513          PYTTOT = PTTOT*TYI
 514          PZTTOT = PTTOT*TZI
 515       GO TO 2999
 516 *  |
 517 *  +-------------------------------------------------------------------*
 518 *  |  AXsi+s are considered as pbars: after the interaction
 519 *  |  a possible way to recover strangeness and charge conservation
 520 *  |  is to flip an eventual pbar into a AXsi+,
 521 *  |  or a nbar --> ASigma+ and a pi- --> K0 (or a pi0 --> K+) etc.
 522 *  |  (we must flip two ubar quarks into sbar quarks)
 523  2037 CONTINUE
 524          LCHTYP = .TRUE.
 525          KPROJ = 2
 526 *  |   The following card must be commented following the strangeness
 527 *  |   and charge conservation recovery implementation
 528          ICHTOT = ICHTOT - 2
 529 *  |   Adjusting mass and momentum for AXsi+s (now pbar)
 530          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 531          PTTOT  = POO
 532          PXTTOT = PTTOT*TXI
 533          PYTTOT = PTTOT*TYI
 534          PZTTOT = PTTOT*TZI
 535       GO TO 2999
 536 *  |
 537 *  +-------------------------------------------------------------------*
 538 *  |  Omega-s are considered as protons: after the interaction there
 539 *  |  are many possible ways to recover strangeness and charge
 540 *  |  conservation (we must flip a d quark into an s quark and two u
 541 *  |  quarks into s quarks).
 542  2038 CONTINUE
 543          LCHTYP = .TRUE.
 544          KPROJ = 1
 545 *  |   The following card must be commented following the strangeness
 546 *  |   and charge conservation recovery implementation
 547          ICHTOT = ICHTOT + 2
 548 *  |   Adjusting mass and momentum for Omega-s (now p)
 549          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 550          PTTOT  = POO
 551          PXTTOT = PTTOT*TXI
 552          PYTTOT = PTTOT*TYI
 553          PZTTOT = PTTOT*TZI
 554       GO TO 2999
 555 *  |
 556 *  +-------------------------------------------------------------------*
 557 *  |  Omegabar+s are considered as pbars: after the interaction there
 558 *  |  are many possible ways to recover strangeness and charge
 559 *  |  conservation (we must flip a dbar quark into an sbar quark and
 560 *  |  two ubar quarks into sbar quarks).
 561  2039 CONTINUE
 562          LCHTYP = .TRUE.
 563          KPROJ = 2
 564 *  |   The following card must be commented following the strangeness
 565 *  |   and charge conservation recovery implementation
 566          ICHTOT = ICHTOT - 2
 567 *  |   Adjusting mass and momentum for AOmega+s (now pbar)
 568          POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
 569          PTTOT  = POO
 570          PXTTOT = PTTOT*TXI
 571          PYTTOT = PTTOT*TYI
 572          PZTTOT = PTTOT*TZI
 573       GO TO 2999
 574 *  |
 575 *  +-------------------------------------------------------------------*
 576  2999 CONTINUE
 577 *  +-------------------------------------------------------------------*
 578 *  |                                   Check for hydrogen!!!!
 579       IF ( IBTAR .EQ. 1 ) THEN
 580          TV = 0.D+00
 581          NP = NP0
 582          ITTA = 8 - 7 * ICHTAR
 583          IF ( KPROJ .EQ. 14 .AND. EKE .LT. 2.D+00 * GAMMIN ) THEN
 584             CALL PMPRAB ( KPROJ, EKE, POO, TXI, TYI, TZI, WE )
 585             RETURN
 586          END IF
 587 *  |   Set a flag to avoid elastic collision in Hadrin
 588          IELFLG = +1
 589 *  |   Set a flag to fully exploit charge exchange
 590          ICXFLG = 0
 591          CALL HADRIV ( KPROJ, POO, ELAB, TXI, TYI, TZI, ITTA )
 592          IELFLG = -1
 593          ICXFLG = -1
 594 *  |  +----------------------------------------------------------------*
 595 *  |  |
 596          IF ( IH .LE. 1 ) THEN
 597 *  |  |  +-------------------------------------------------------------*
 598 *  |  |  |
 599             IF ( ITH (1) .EQ. KPROJ ) THEN
 600                IH  = IH + 1
 601                ITH (2) = 1
 602                CXH (2) = TXI
 603                CYH (2) = TYI
 604                CZH (2) = TZI
 605                PLH (2) = 0.D+00
 606                ELH (2) = AM (1)
 607 *  |  |  |
 608 *  |  |  +-------------------------------------------------------------*
 609 *  |  |  |
 610             ELSE IF ( ITH (1) .EQ. 1 ) THEN
 611                IH  = IH + 1
 612                ITH (2) = KPROJ
 613                CXH (2) = TXI
 614                CYH (2) = TYI
 615                CZH (2) = TZI
 616                PLH (2) = 0.D+00
 617                ELH (2) = AM (KPROJ)
 618 *  |  |  |
 619 *  |  |  +-------------------------------------------------------------*
 620 *  |  |  |
 621             ELSE
 622             END IF
 623 *  |  |  |
 624 *  |  |  +-------------------------------------------------------------*
 625           END IF
 626 *  |  |
 627 *  |  +----------------------------------------------------------------*
 628
 629 *  |  +----------------------------------------------------------------*
 630 *  |  |                  Looping over the particles produced in hadrin
 631          DO 110 J=1,IH
 632             NP = NP + 1
 633             TKI(NP) = ELH(J) - AM(ITH(J))
 634 *  |  |  +-------------------------------------------------------------*
 635 *  |  |  |
 636             IF ( TKI(NP) .GE. 0.D+00 ) THEN
 637                KPART(NP) = ITH(J)
 638                CXR(NP)   = CXH(J)
 639                CYR(NP)   = CYH(J)
 640                CZR(NP)   = CZH(J)
 641                PLR(NP)   = PLH(J)
 642                WEI(NP)   = WE
 643 *  |  |  |  Updating conservation counters
 644                ENUCR  = ENUCR + ELH(J)
 645                PXNUCR = PXNUCR + PLR(NP)*CXR(NP)
 646                PYNUCR = PYNUCR + PLR(NP)*CYR(NP)
 647                PZNUCR = PZNUCR + PLR(NP)*CZR(NP)
 648                ICNUCR = ICNUCR + ICH(KPART(NP))
 649                IBNUCR = IBNUCR + IBAR(KPART(NP))
 650 *  |  |  |
 651 *  |  |  +-------------------------------------------------------------*
 652 *  |  |  |
 653             ELSE IF ( ABS ( TKI (NP) ) .LE. ANGLGB ) THEN
 654 *  |  |  |  Updating conservation counters
 655                ENUCR  = ENUCR  + ELH(J)
 656                ICNUCR = ICNUCR + ICH(ITH(J))
 657                IBNUCR = IBNUCR + IBAR(ITH(J))
 658                NP = NP - 1
 659 *  |  |  |
 660 *  |  |  +-------------------------------------------------------------*
 661 *  |  |  |
 662             ELSE
 663                NP = NP - 1
 664             END IF
 665 *  |  |  |
 666 *  |  |  +-------------------------------------------------------------*
 667  110     CONTINUE
 668 *  |  |
 669 *  |  +----------------------------------------------------------------*
 670          KTARP = 1
 671          ANOW = 0.D+00
 672          ZNOW = 0.D+00
 673          ICRES = 0
 674          IBRES = 0
 675          ERES  = 0.D+00
 676          EOTEST = EOTEST - ENUCR
 677 *  |  +----------------------------------------------------------------*
 678 *  |  |
 679          IF ( ABS (EOTEST) .GT. ETEPS ) THEN
 680             WRITE (LUNERR,*)' Eventq: eotest failure with Hadrin',
 681      &                        EOTEST
 682             LRESMP = .TRUE.
 683             GO TO 50
 684 *  |
 685 *  +-<|--<--<--<--<--< go to resampling
 686          END IF
 687 *  |  |
 688 *  |  +----------------------------------------------------------------*
 689          RETURN
 690       END IF
 691 *  |                                       end hydrogen check
 692 *  +-------------------------------------------------------------------*
 693 *
 694 *  Now calling Nucrin to produce secondary particles, which are put
 695 *  into the /Finuc/ common
 696 *
 697       CALL NUCRIV ( KPROJ, ELAB, TXI, TYI, TZI, BBTAR,
 698      &              ZZTAR, RHO(IMAT) )
 699       IF ( LRESMP ) GO TO 50
 700 *--<--<--<--<--< go to resampling if something was wrong
 701       TVGRE0 = TV - TVGREY
 702       TVEUZ  = 0.D+00
 703       TVTENT = TV
 704       TV = 0.D+00
 705       IF (NP .GT. MXP) GO TO 400
 706 *  +-------------------------------------------------------------------*
 707 *  |                     Looping over the particles produced in nucrin
 708 *  |                     No need for Nucrin produced particles to
 709 *  |                     change the numbering
 710       DO 101 I=NP0+1,NP
 711          I1 = KPART(I)
 712          WEI(I) = WE
 713          TKI(I) = TKI(I) - AM(I1)
 714          TXYZ   = CXR (I)**2 + CYR (I)**2 + CZR (I)**2
 715          IF ( ABS (TXYZ-1.D0) .GT. CSNPRN ) THEN
 716             WRITE ( LUNERR,* )
 717      &      ' **** Bad normalized cosines from Nucriv!!!',I,TXYZ,CXR(I),
 718      &        CYR(I),CZR(I)
 719             TXYZ = TXYZ - 1.D+00
 720             TXYZ = 1.D+00 - 0.5D+00 * TXYZ + 0.375D+00 * TXYZ * TXYZ
 721             CXR (I) = CXR (I) * TXYZ
 722             CYR (I) = CYR (I) * TXYZ
 723             CZR (I) = CZR (I) * TXYZ
 724             TXYZ   = CXR (I)**2 + CYR (I)**2 + CZR (I)**2
 725             IF ( ABS (TXYZ-1.D0) .GT. CSNPRN ) THEN
 726                LRESMP = .TRUE.
 727                GO TO 50
 728             END IF
 729          END IF
 730          IF (TKI(I) .GE. 0.D+00) GO TO 101
 731             TKI(I) = 0.D+00
 732             PLR(I) = 0.D+00
 733  101  CONTINUE
 734 *  |
 735 *  +-------------------------------------------------------------------*
 736       GO TO 500
 737 C
 738 C********************************************************************
 739 C     GENERATE FIRST THE LOW ENERGY NUCLEONS FROM THE INTRANUCLEAR
 740 C     CASCADE - CHECK IF ACCEPTABLE.
 741 C     NOTE THAT AINEL AND RAKEKA ARE USING THE OLD PARTICLE NUMBERING
 742 C********************************************************************
 743 C
 744  200  CONTINUE
 745       NP = NP0
 746       LNUCRI = .FALSE.
 747 *  +-------------------------------------------------------------------*
 748 *  |                                   Check for hydrogen!!!!
 749       IF ( IBTAR .EQ. 1 ) THEN
 750          ANCOLL = 1.D+00
 751          NHAD = 0
 752          KPROJ = IJ
 753          KTARG = 1
 754          EPROJ = EKE + AM (KPROJ)
 755          UMO   = SQRT ( AM (1)**2 + AM (KPROJ)**2 + 2.D+00 * AM (1) *
 756      &                  EPROJ )
 757 *  |  +----------------------------------------------------------------*
 758 *  |  | Now diffractive events are taken into account for projectile
 759 *  |  | -proton collisions!!! A. Ferrari and J. Ranft, 25-3-90
 760          IF ( LDIFFR (KPTOIP(KPROJ)) ) THEN
 761             CALL GRNDM(RNDM,1)
 762             IF ( RNDM (1) .LT. FRDIFF ) THEN
 763                IF ( POO .GT. PTHDFF ) THEN
 764                   LEVDIF = .TRUE.
 765                   CALL DIFEVV ( NHAD, KPROJ, KTARG, POO, EPROJ, UMO )
 766                ELSE
 767                   CALL HADEVV ( NHAD, KPROJ, KTARG, POO, EPROJ, UMO )
 768                END IF
 769             ELSE
 770                CALL HADEVV ( NHAD, KPROJ, KTARG, POO, EPROJ, UMO )
 771             END IF
 772 *  |  |
 773 *  |  +----------------------------------------------------------------*
 774 *  |  |
 775          ELSE
 776             CALL HADEVV ( NHAD, KPROJ, KTARG, POO, EPROJ, UMO )
 777          END IF
 778 *  |  |
 779 *  |  +----------------------------------------------------------------*
 780 *  |  +----------------------------------------------------------------*
 781 *  |  |                   Loop over particles produced in hadevt/difevt
 782          DO 207 J=1,NHAD
 783             NP  = NP + 1
 784             CFE = 1.D+00
 785             SFE = 0.D+00
 786 *  |  |   Get the "paprop" index
 787             KPART(NP) = KPTOIP (NREH(J))
 788             IF ( KPART (NP) .EQ. 0 ) THEN
 789                WRITE (LUNERR,*)
 790      &         ' **** Charmed particle produced in Hadevv/Difevv',
 791      &         NREH(J),HEPH(J),AMH(J),' ****'
 792                KPART (NP) = NREH (J)
 793             END IF
 794             PLR(NP) = SQRT ( PXH(J)**2 + PYH(J)**2 + PZH(J)**2 )
 795             PTH = SQRT( PXH(J)**2 + PYH(J)**2 )
 796             CDE = PZH(J) / PLR(NP)
 797             SDE = PTH    / PLR(NP)
 798             IF (SDE .GE. ANGLGB) THEN
 799                CFE = PXH(J) / PTH
 800                SFE = PYH(J) / PTH
 801             END IF
 802             CALL TTRANS ( TXI, TYI, TZI, CDE, SDE, SFE, CFE,
 803      &                    CXR(NP), CYR(NP), CZR(NP) )
 804             TKI(NP) = HEPH(J) - AMH(J)
 805             WEI(NP) = WE
 806 *  |  |  +-------------------------------------------------------------*
 807 *  |  |  |                   Updating conservation counters
 808             IF ( TKI(NP) .GT. 0.D+00 ) THEN
 809                EUZ = EUZ + HEPH(J)
 810                PUX = PUX + PLR(NP)*CXR(NP)
 811                PUY = PUY + PLR(NP)*CYR(NP)
 812                PUZ = PUZ + PLR(NP)*CZR(NP)
 813                ICU = ICU + ICH(NREH(J))
 814                IBU = IBU + IBAR(NREH(J))
 815 *  |  |  |
 816 *  |  |  +-------------------------------------------------------------*
 817 *  |  |  |
 818             ELSE
 819                NP = NP-1
 820             END IF
 821 *  |  |  |
 822 *  |  |  +-------------------------------------------------------------*
 823  207     CONTINUE
 824 *  |  |                                end loop
 825 *  |  +----------------------------------------------------------------*
 826          ANOW = 0.D+00
 827          ZNOW = 0.D+00
 828          KTARP = 1
 829          ICRES = 0
 830          IBRES = 0
 831          ERES  = 0.D+00
 832          EOTEST = EOTEST - EUZ
 833 *  |  +----------------------------------------------------------------*
 834 *  |  |
 835          IF ( ABS (EOTEST) .GT. ETEPS ) THEN
 836             WRITE (LUNERR,*)' Eventq: eotest failure with',
 837      &                      ' Hadevt/Diffevt',EOTEST
 838          END IF
 839 *  |  |
 840 *  |  +----------------------------------------------------------------*
 841          RETURN
 842       END IF
 843 *  |                                       end hydrogen check
 844 *  +-------------------------------------------------------------------*
 845 *  Now at first we sample the number of interactions
 846       CALL COREVT ( ZZTAR, BBTAR, IJ, POO, EKE )
 847 *
 848 *  First decide how much energy will go into the intranuclear cascade
 849 *
 850       EKIN = EKE
 851 *  Here starts the cascade particle production module!!!
 852       SOPPP (1) = EINCP
 853       SOPPP (2) = EINCN
 854       ARECL  = MAX ( ANOW - 1.D+00, 1.D+00 )
 855       KREJE  = 0
 856       NGREYT = NGREYP + NGREYN
 857       TVTENT = ANCOLL * ( AV0WEL - AEFRMA ) + TVGRE0
 858       TMP015 = 0.15D+00
 859       EEXTRA = 2.D+00* MIN ( TWOTWO * ( EKUPNU (1) + EKUPNU (2) )
 860      &                     , TMP015 *
 861      &       ( EKE - EINCP - EINCN - TVTENT ), EINCP + EINCN )
 862       EEXTRA = MAX ( EEXTRA, EKUPNU (1), EKUPNU (2) )
 863       LINCCK = .FALSE.
 864       PXCASC = 0.D+00
 865       PYCASC = 0.D+00
 866       PZCASC = 0.D+00
 867       PTXINT = 0.D+00
 868       PTYINT = 0.D+00
 869       PPINTR = 0.D+00
 870       EKRECL = 0.D+00
 871       IEXTRN = 0
 872       IEXTRP = 0
 873       LACCEP = .FALSE.
 874       IGREYT = 0
 875       IF ( NGREYT .EQ. 0 ) GO TO 206
 876 *  +-------------------------------------------------------------------*
 877 *  | Now looping until we reach the correct number of cascade
 878 *  | nucleons
 879  205  CONTINUE
 880          IGREYT = IGREYP + IGREYN
 881          DSOPP  = SOPPP (1) + SOPPP (2)
 882 *  |  Now it is not possible to reduce too heavily the available energy
 883          IF ( EKIN - TVGRE0 .LT. 0.5D+00 * EKE .OR.  PZCASC .GE.
 884      &        0.5D+00 * POO ) GO TO 206
 885 *  |  +----------------------------------------------------------------*
 886 *  |  |
 887          IF ( IGREYT .GE. NGREYT ) THEN
 888             AEXTRA = 2.D+00 * DSOPP / ( EKINAV (1) + EKINAV (2)
 889      &             + ESWELL (1) + ESWELL (2) )
 890      &             / ( IGREYT - NGREYT + 1 )
 891             CALL GRNDM(RNDM,1)
 892             RNDUMO = RNDM (1)
 893 *  |  |  +-------------------------------------------------------------*
 894 *  |  |  |
 895             IF ( RNDUMO .LT. AEXTRA .AND. LACCEP .AND. NINT (ANOW)
 896      &           .GT. 0 ) THEN
 897                REJE = ZNOW / ANOW
 898 *  |  |  |  +----------------------------------------------------------*
 899 *  |  |  |  |
 900                CALL GRNDM(RNDM,1)
 901                IF ( RNDM (1) .LT. REJE ) THEN
 902                   IEXTRP = MIN ( 1, NINT (ZNOW), NGREYP / 2 )
 903                   IF ( NGREYP + IEXTRP - IGREYP .LE. 0 ) GO TO 206
 904 *  |  |  |  |
 905 *  |  |  |  +----------------------------------------------------------*
 906 *  |  |  |  |
 907                ELSE
 908 *                 IEXTRN = MIN ( IEXTRN + 1, NINT (ANOW - ZNOW), 2 )
 909                   IEXTRN = MIN ( 1, NINT (ANOW - ZNOW), NGREYN / 2 )
 910                   IF ( NGREYN + IEXTRN - IGREYN .LE. 0 ) GO TO 206
 911                END IF
 912 *  |  |  |  |
 913 *  |  |  |  +----------------------------------------------------------*
 914 *  |  |  |
 915 *  |  |  +-------------------------------------------------------------*
 916 *  |  |  |
 917             ELSE
 918                GO TO 206
 919 *  |  |  |
 920 *  |  |  |-->-->-->-->--> We have finished particles!!
 921             END IF
 922 *  |  |  |
 923 *  |  |  +-------------------------------------------------------------*
 924 *  |  |
 925 *  |  +----------------------------------------------------------------*
 926 *  |  |
 927          ELSE IF ( DSOPP .LE. - EEXTRA ) THEN
 928             GO TO 206
 929 *  |  |
 930 *  |  |-->-->-->-->--> We have finished energy
 931          END IF
 932 *  |  |
 933 *  |  +----------------------------------------------------------------*
 934          LACCEP = .FALSE.
 935 *  |  +----------------------------------------------------------------*
 936 *  |  | Check how many nucleons are still available for cascade
 937          IF ( ANOW .LT. 1.5D+00 ) THEN
 938 *  |  |  +-------------------------------------------------------------*
 939 *  |  |  | Check if at least two high energy collisions are foreseen
 940 *  |  |  | if yes the proper energy / momentum balance for the last
 941 *  |  |  | two nucleons will be done inside Ferevv
 942             IF ( ANCOLL .GT. 1.5D+00 ) THEN
 943 *  |  |  |
 944 *  |  |  +-------------------------------------------------------------*
 945 *  |  |  |     Only the valence collision is foreseen: we are forced
 946 *  |  |  |     to perform here the balance for the 2 last nucleons!
 947 *  |  |  |     (One might be used here for cascading, the other
 948 *  |  |  |     will be used by Hadevv / Difevv)
 949             ELSE
 950                LLASTN = .TRUE.
 951                KTLAST = IJTARG (1)
 952                AMLAST = AM ( KTLAST )
 953 *  |  |  |  +----------------------------------------------------------*
 954 *  |  |  |  |
 955                IF ( NINT ( ZNOW ) .GT. 0 ) THEN
 956                   KTINC = 1
 957                   IGREYP = IGREYP + 1
 958 *  |  |  |  |
 959 *  |  |  |  +----------------------------------------------------------*
 960 *  |  |  |  |
 961                ELSE
 962                   KTINC = 8
 963                   IGREYN = IGREYN + 1
 964                END IF
 965 *  |  |  |  |
 966 *  |  |  |  +----------------------------------------------------------*
 967                AMINC = AM (KTINC)
 968                UMIN2 = ( AMLAST + AMINC )**2
 969                DSOPP = MAX ( DSOPP, AV0WEL - AEFRMA )
 970   216          CONTINUE
 971                EKIN   = EKIN  - DSOPP
 972 *  |  |  |  +----------------------------------------------------------*
 973 *  |  |  |  |  Check if we can get enough invariant mass
 974                IF ( EKIN .LE. 0.2D+00 * EKE ) THEN
 975                   WRITE ( LUNOUT, * )
 976      &            ' *** Eventv: impossible to get enough invariant',
 977      &            ' mass for the last two nucleons! ****'
 978                   WRITE ( LUNERR, * )
 979      &            ' *** Eventv: impossible to get enough invariant',
 980      &            ' mass for the last two nucleons! ****'
 981                   LRESMP = .TRUE.
 982                   GO TO 50
 983 *  |  |  |  |---<---<---<  Go to resampling
 984                END IF
 985 *  |  |  |  |
 986 *  |  |  |  +----------------------------------------------------------*
 987                ELEFT  = ETTOT - EINTR  - EKIN - AM (IJ)
 988                PLA    = SQRT ( EKIN * ( EKIN + 2.D+00 * AM (IJ) ) )
 989                PXLEFT = PXTTOT - PXINTR - PLA * TXI
 990                PYLEFT = PYTTOT - PYINTR - PLA * TYI
 991                PZLEFT = PZTTOT - PZINTR - PLA * TZI
 992 *  |  |  |  Now we will divide Eleft and Pleft between the two
 993 *  |  |  |  left nucleons!
 994                UMO2 = ELEFT**2 - PXLEFT**2 - PYLEFT**2 - PZLEFT**2
 995 *  |  |  |  +----------------------------------------------------------*
 996 *  |  |  |  |
 997                IF ( UMO2 .LE. UMIN2 ) THEN
 998 *  |  |  |  |  +-------------------------------------------------------*
 999 *  |  |  |  |  |
1000                   IF ( KTINC .EQ. 1 ) THEN
1001                      EINCP = EINCP + DSOPP
1002 *  |  |  |  |  |
1003 *  |  |  |  |  +-------------------------------------------------------*
1004 *  |  |  |  |  |
1005                   ELSE
1006                      EINCN = EINCN + DSOPP
1007                   END IF
1008 *  |  |  |  |  |
1009 *  |  |  |  |  +-------------------------------------------------------*
1010                   GO TO 216
1011                END IF
1012 *  |  |  |  |
1013 *  |  |  |  +----------------------------------------------------------*
1014                UMO = SQRT ( UMO2 )
1015                ELACMS = 0.5D+00 * ( UMO2 + AMLAST**2 - AMINC**2 )
1016      &                / UMO
1017                EINCMS = UMO - ELACMS
1018                PCMS   = SQRT ( ELACMS**2 - AMLAST**2 )
1019                GAMCM = ELEFT  / UMO
1020                ETAX  = PXLEFT / UMO
1021                ETAY  = PYLEFT / UMO
1022                ETAZ  = PZLEFT / UMO
1023                CALL RACO ( CXXINC, CYYINC, CZZINC )
1024                PCMSX = PCMS * CXXINC
1025                PCMSY = PCMS * CYYINC
1026                PCMSZ = PCMS * CZZINC
1027 *  |  |  |  Now go back from the CMS frame to the lab frame!!!
1028 *  |  |  |  First the inc nucleon:
1029                ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
1030                NP = NP + 1
1031                IF (NP .GT. MXP) GO TO 400
1032                KPART(NP)= KTINC
1033                WEI (NP) = WE
1034                TKI (NP) = GAMCM * EINCMS + ETAPCM - AMINC
1035                PHELP  = ETAPCM / (GAMCM + 1.D0) + EINCMS
1036                CXXINC = PCMSX + ETAX * PHELP
1037                CYYINC = PCMSY + ETAY * PHELP
1038                CZZINC = PCMSZ + ETAZ * PHELP
1039 *  |  |  |  Updating conservation counters
1040                EINTR  = EINTR  + TKI(NP) + AMINC
1041                PXINTR = PXINTR + CXXINC
1042                PYINTR = PYINTR + CYYINC
1043                PZINTR = PZINTR + CZZINC
1044                ICINTR = ICINTR + ICH  (KPART(NP))
1045                IBINTR = IBINTR + IBAR (KPART(NP))
1046                PLR (NP) = SQRT (CXXINC * CXXINC + CYYINC * CYYINC +
1047      &                          CZZINC * CZZINC)
1048                CXR (NP) = CXXINC / PLR (NP)
1049                CYR (NP) = CYYINC / PLR (NP)
1050                CZR (NP) = CZZINC / PLR (NP)
1051 *  |  |  |  Now the high energy collision nucleon
1052                EKLAST = GAMCM * ELACMS - ETAPCM - AMLAST
1053                PHELP  = - ETAPCM / (GAMCM + 1.D0) + ELACMS
1054                PXLAST = - PCMSX + ETAX * PHELP
1055                PYLAST = - PCMSY + ETAY * PHELP
1056                PZLAST = - PCMSZ + ETAZ * PHELP
1057 *  |  |  |  Now we perform a Lorentz transformation to the rest fra-
1058 *  |  |  |  me of the "last" nucleon, with the projectile momentum
1059 *  |  |  |  along the +z axis
1060                AMPROJ = AM (IJ)
1061                EPROJ  = EKIN + AMPROJ
1062                ECHCK  = EPROJ + AMLAST + EKLAST
1063                PXCHCK = PLA * TXI + PXLAST
1064                PYCHCK = PLA * TYI + PYLAST
1065                PZCHCK = PLA * TZI + PZLAST
1066                UMO  = SQRT ( ECHCK**2 - PXCHCK**2 - PYCHCK**2 -
1067      &                       PZCHCK**2 )
1068                EPROJX = 0.5D+00 * ( UMO**2 - AMPROJ**2 - AMLAST**2 )
1069      &                / AMLAST
1070                PPROJX = SQRT ( EPROJX**2 - AMPROJ**2 )
1071                ETOTX  = EPROJX + AMLAST
1072 *  |  |  |  Now set the parameters for the Lorentz transformation
1073                AAFACT = ECHCK  + ETOTX
1074                BBFACT = PPROJX - PZCHCK
1075                DDENOM = ETOTX * AAFACT - PPROJX * BBFACT
1076                GAMTRA = ( ECHCK * AAFACT + PPROJX * BBFACT ) / DDENOM
1077                ETAZ = - BBFACT * AAFACT / DDENOM
1078                ETAX = PXCHCK * ( GAMTRA + 1.D+00 ) / AAFACT
1079                ETAY = PYCHCK * ( GAMTRA + 1.D+00 ) / AAFACT
1080                PLABS = PPROJX
1081                ELABS = EPROJX
1082 *  |  |  |  +----------------------------------------------------------*
1083 *  |  |  |  |
1084                IF ( .NOT. LDIFFR (KPTOIP(IJ)) .OR. PLABS .LE. PTHDFF )
1085      &            THEN
1086                   RUUN = 1.D0
1087 *  |  |  |  |
1088 *  |  |  |  +----------------------------------------------------------*
1089 *  |  |  |  |
1090                ELSE
1091                   CALL GRNDM(RNDM,1)
1092                   RUUN = RNDM (1)
1093                END IF
1094 *  |  |  |  |
1095 *  |  |  |  +----------------------------------------------------------*
1096                NHAD = 0
1097 *  |  |  |  +----------------------------------------------------------*
1098 *  |  |  |  |   Diffractive event
1099                IF ( RUUN .LE. FRDIFF ) THEN
1100                   LEVDIF = .TRUE.
1101                   CALL DIFEVV ( NHAD, IJ, KTLAST, PLABS, ELABS, UMO )
1102 *  |  |  |  |
1103 *  |  |  |  +----------------------------------------------------------*
1104 *  |  |  |  |   Usual event
1105                ELSE
1106                   CALL HADEVV ( NHAD, IJ, KTLAST, PLABS, ELABS, UMO )
1107                END IF
1108 *  |  |  |  |
1109 *  |  |  |  +----------------------------------------------------------*
1110 *  |  |  |  Now we have the particles in the nucleon rest frame,
1111 *  |  |  |  transform back in the lab frame
1112 *  |  |  |  +----------------------------------------------------------*
1113 *  |  |  |  |           Looping over the produced particles
1114                DO 236 I = 1, NHAD
1115                   NP = NP + 1
1116                   IF (NP .GT. MXP) GO TO 400
1117                   CALL ALTRA ( GAMTRA, ETAX, ETAY, ETAZ, PXH (I),
1118      &                         PYH (I), PZH (I), HEPH (I), PLR (NP),
1119      &                         CXR (NP), CYR (NP), CZR (NP), ELR )
1120 *  |  |  |  |  Updating conservation counters
1121                   EUZ  = EUZ + ELR
1122                   PUX  = PUX + CXR (NP)
1123                   PUY  = PUY + CYR (NP)
1124                   PUZ  = PUZ + CZR (NP)
1125                   KPART (NP) = KPTOIP ( NREH (I) )
1126                   IF ( KPART (NP) .EQ. 0 ) THEN
1127                      WRITE (LUNERR,*)
1128      &               ' **** Charmed particle produced in Hadevv',
1129      &               NREH(I),ELR,AMH(I),' ****'
1130                      KPART (NP) = NREH (I)
1131                   END IF
1132                   ICU = ICU + ICH  (NREH(I))
1133                   IBU = IBU + IBAR (NREH(I))
1134                   CXR (NP) = CXR (NP) / PLR (NP)
1135                   CYR (NP) = CYR (NP) / PLR (NP)
1136                   CZR (NP) = CZR (NP) / PLR (NP)
1137                   TKI (NP) = ELR - AMH (I)
1138                   WEI (NP) = WE
1139   236          CONTINUE
1140 *  |  |  |  |
1141 *  |  |  |  +----------------------------------------------------------*
1142 *  |  |  | Now perform the standard tests for energy and momentum
1143 *  |  |  | conservation
1144                ERES  = ETTOT  - EUZ - ENUCR  - EINTR
1145                PXRES = PXTTOT - PUX - PXNUCR - PXINTR
1146                PYRES = PYTTOT - PUY - PYNUCR - PYINTR
1147                PZRES = PZTTOT - PUZ - PZNUCR - PZINTR
1148                ICRES = ICHTOT - ICU - ICNUCR - ICINTR
1149                IBRES = IBTOT  - IBU - IBNUCR - IBINTR
1150                PTRES = MAX ( ABS (PXRES), ABS (PYRES), ABS (PZRES) )
1151 *  |  |  |  +----------------------------------------------------------*
1152 *  |  |  |  |
1153                IF ( IBRES .NE. 0 .OR. ICRES .NE. 0 .OR. ABS (ERES)
1154      &              .GT. 1.D-10*EPROJ .OR. PTRES .GT. 1.D-10*PTTOT )
1155      &              THEN
1156                   WRITE ( LUNERR, * )
1157      &                   ' Eventq: last nucleon failure!!', ICRES,
1158      &                     IBRES, REAL (ERES), REAL (PXRES),
1159      &                     REAL (PYRES), REAL (PZRES)
1160                   LRESMP = .TRUE.
1161                   GO TO 50
1162 *  |  |  |--|--<--<--<--<--< go to resampling
1163                END IF
1164 *  |  |  |  |
1165 *  |  |  |  +----------------------------------------------------------*
1166                PTRES = 0.D+00
1167                RETURN
1168             END IF
1169 *  |  |  |
1170 *  |  |  +-------------------------------------------------------------*
1171          END IF
1172 *  |  |
1173 *  |  +----------------------------------------------------------------*
1174
1175 *  |  +----------------------------------------------------------------*
1176 *  |  |
1177          IF ( NINT ( ANOW - ZNOW ) .LE. 0 ) THEN
1178             REJE = 1.D+00
1179 *  |  |
1180 *  |  +----------------------------------------------------------------*
1181 *  |  |
1182          ELSE
1183             REJE = MAX ( ZNOW * ( NGREYP + IEXTRP - IGREYP ),
1184      &                   ZERZER )
1185             HELP = REJE + ( ANOW - ZNOW ) * ( NGREYN +
1186      &             IEXTRN - IGREYN )
1187             IF ( HELP .GT. 0.D+00 ) THEN
1188                REJE = REJE / HELP
1189             ELSE
1190                REJE = 1.D+00
1191                WRITE(LUNERR,*)' EVENTV: HELP=0, ANOW,ZNOW',ANOW,ZNOW
1192      &                       ,' NGREYN,IGREYN,IEXTRN',
1193      &                          NGREYN,IGREYN,IEXTRN,
1194      &                        ' NGREYP,IGREYP,IEXTRP',
1195      &                          NGREYP,IGREYP,IEXTRP
1196             END IF
1197          END IF
1198 *  |  |
1199 *  |  +----------------------------------------------------------------*
1200
1201 *  |  +----------------------------------------------------------------*
1202 *  |  |
1203          CALL GRNDM(RNDM,1)
1204          IF ( RNDM (1) .LT. REJE ) THEN
1205             ILO  = 1
1206             ILLO = 2
1207             KP   = 1
1208 *  |  |
1209 *  |  +----------------------------------------------------------------*
1210 *  |  |
1211          ELSE
1212             ILO  = 2
1213             ILLO = 1
1214             KP   = 8
1215          END IF
1216 *  |  |
1217 *  |  +----------------------------------------------------------------*
1218          CALL RAKEKV ( ILO, EXSOP, BBTAR, TKIN, TSTRCK, PLA, ARECL,
1219      &                 TKRECL, EFERMI, CDE, SDE )
1220          IF ( EKIN - TKIN .LT. 0.5D+00 * EKE ) GO TO 206
1221 *  |--|--<--<--<--<--< avoid to deplete too much the energy
1222 *  |  +----------------------------------------------------------------*
1223 *  |  |  Now check if the nucleon energy is acceptable:
1224          IF ( SOPPP (ILO) .LT. TKIN )  THEN
1225 *  |  |  +-------------------------------------------------------------*
1226 *  |  |  |
1227             IF ( TKIN - SOPPP (ILO) .LT. SOPPP (ILLO) + 1.5D+00 *
1228      &           EEXTRA ) THEN
1229                SOPPP (ILLO) = SOPPP (ILLO) + SOPPP (ILO) - TKIN
1230      &                      - EXSOP (ILO)
1231                SOPPP (ILO)  = 0.D+00
1232 *  |  |  |
1233 *  |  |  +-------------------------------------------------------------*
1234 *  |  |  |
1235             ELSE
1236                KREJE = KREJE + 1
1237                IF ( KREJE .GT. 10 ) GO TO 206
1238                SOPPP (ILLO) = SOPPP (ILLO) + SOPPP (ILO)
1239                SOPPP (ILO)  = 0.D+00
1240                GO TO 205
1241             END IF
1242 *  |  |  |
1243 *  |  |  +-------------------------------------------------------------*
1244 *  |  |
1245 *  |  +----------------------------------------------------------------*
1246 *  |  |
1247          ELSE
1248             SOPPP (ILO) = SOPPP (ILO) - TKIN - EXSOP (ILO)
1249          END IF
1250 *  |  |
1251 *  |  +----------------------------------------------------------------*
1252          IF ( NINT ( ANOW ) .LE. 0 ) GO TO 206
1253 *  |  +----------------------------------------------------------------*
1254 *  |  |   Update the current atomic and mass number
1255          IF ( ILO .EQ. 1 ) THEN
1256             IF ( NINT ( ZNOW ) .LE. 0 ) GO TO 206
1257             IGREYP = IGREYP + 1
1258             ZNOW   = ZNOW   - 1.D+00
1259 *  |  |
1260 *  |  +----------------------------------------------------------------*
1261 *  |  |
1262          ELSE
1263             IGREYN = IGREYN + 1
1264          END IF
1265 *  |  |
1266 *  |  +----------------------------------------------------------------*
1267          ANOW = ANOW - 1.D+00
1268 *  |  Now make the kinematical tests!!
1269          FRSURV = ARECL / BBTAR
1270          NP = NP + 1
1271          IF ( NP .GT. MXP ) GO TO 400
1272          EKIN   = EKIN - TKIN
1273 *  Note the change!!
1274          LINCCK = ARECL .LT. 30.D+00 .AND. FRSURV .LE. 0.33D+00
1275 *  |  +----------------------------------------------------------------*
1276 *  |  |
1277          IF ( .NOT. LINCCK ) THEN
1278             CALL SFECFE ( SFE, CFE )
1279 *  |  |
1280 *  |  +----------------------------------------------------------------*
1281 *  |  |
1282          ELSE
1283             IF ( CDE .GT. 0.D+00 ) THEN
1284                CDE = CDE * FRSURV / 0.33D+00
1285                SDE = SQRT ( 1.D+00 - CDE**2 )
1286             END IF
1287             CALL SFECFE ( SFE, CFE )
1288          END IF
1289 *  |  |
1290 *  |  +----------------------------------------------------------------*
1291          PTXINT = PTXINT + PLA * CFE * SDE
1292          PTYINT = PTYINT + PLA * SFE * SDE
1293          PPINTR = PPINTR + PLA * CDE
1294          CALL TTRANS ( TXI, TYI, TZI, CDE, SDE, SFE, CFE,
1295      &                 CXR (NP), CYR (NP), CZR (NP) )
1296 *  | Generated nucleon is acceptable - save it
1297          LACCEP = .TRUE.
1298          TVGRE0 = TVGRE0 + EXSOP (ILO)
1299          TVGREY = TVGREY + TKIN + TKRECL - TSTRCK - EBNDNG (ILO)
1300          EKRECL = EKRECL + TKRECL
1301          ARECL  = MAX ( ARECL - 1.D+00, 1.D+00 )
1302          KPART(NP)= KP
1303          WEI(NP)  = WE
1304          TKI(NP)  = TSTRCK
1305          PLR(NP)  = PLA
1306 *  |                        Updating the counters for spent momenta
1307 *  |                        of the cascade nucleons (the momentum
1308 *  |                        spent by the high energy particles)
1309          SINTH = SQRT ( 1.D+00 - COSTH * COSTH )
1310          PLA = SQRT ( ( EFERMI + TKIN )**2 - AMNUSQ (ILO) )
1311          PZLEFT = PLA * CDE - COSTH * PKFRMI
1312          IF ( PZLEFT + PZCASC .LT. EKE - EKIN + TVGRE0 ) THEN
1313             PZCASC = EKE - EKIN + TVGRE0
1314             PXCASC = PXCASC + CFE * ( PLA * SDE - SINTH * PKFRMI )
1315             PYCASC = PYCASC + SFE * ( PLA * SDE - SINTH * PKFRMI )
1316          ELSE
1317             PZCASC = PZCASC + PZLEFT
1318             PXCASC = PXCASC + CFE * ( PLA * SDE - SINTH * PKFRMI )
1319             PYCASC = PYCASC + SFE * ( PLA * SDE - SINTH * PKFRMI )
1320          END IF
1321 *  |                        Updating conservation counters
1322          EINTR  = EINTR  + TKI(NP) + AM (KPART(NP))
1323          PXINTR = PXINTR + PLR(NP) * CXR(NP)
1324          PYINTR = PYINTR + PLR(NP) * CYR(NP)
1325          PZINTR = PZINTR + PLR(NP) * CZR(NP)
1326          ICINTR = ICINTR + ICH  (KPART(NP))
1327          IBINTR = IBINTR + IBAR (KPART(NP))
1328       GO TO 205
1329 *  |
1330 *  +--<--<--<--<--< go to sample another nucleon
1331  206  CONTINUE
1332 *  +-------------------------------------------------------------------*
1333 *  |  First check charge and baryonic number conservation:
1334       IF ( IGREYP .NE. ICINTR .OR. ( IGREYP + IGREYN ) .NE.
1335      &     IBINTR ) THEN
1336          WRITE ( LUNERR,* )' Eventq: charge or baryon number',
1337      &                     ' conservation failure in the Inc',
1338      &                     ' section!!', IGREYP, ICINTR,
1339      &                       IGREYP + IGREYN, IBINTR
1340          LRESMP = .TRUE.
1341          GO TO 50
1342 *--|--<--<--<--<--< go to resampling
1343       END IF
1344 *  |
1345 *  +-------------------------------------------------------------------*
1346       TVTENT = TVGRE0
1347 * Tvgrey is already inside Eincp and Eincn since Rakekv updates the
1348 * estimated excitation energy (and Tvgre0, Eincp, Eincn and the
1349 * Soppp array)
1350       EINCP  = EINCP - SOPPP (1)
1351       EINCN  = EINCN - SOPPP (2)
1352       EKIN   = EKIN  - TVTENT
1353 *  +-------------------------------------------------------------------*
1354 *  |   Check if we have not spent too much energy!!!
1355       IF ( EKIN .LT. 0.5D+00 * EKE ) THEN
1356          LRESMP = .TRUE.
1357          WRITE ( LUNERR,* ) ' Eventq: Ekin after inc too low!! ',
1358      &                        EKIN, EKE, IJIJ
1359       END IF
1360 *  |
1361 *  +-------------------------------------------------------------------*
1362       IF ( LRESMP ) GO TO 50
1363 *
1364 *--<--<--<--<--<--< go to resampling if something was wrong
1365 *   Here the modification to take into account properly the
1366 *   cascade nucleon momenta
1367 *   New version: rotate the spent momentum in the lab frame
1368       IF ( IGREYT .GT. 0 ) THEN
1369          PZCASC = MAX ( PZCASC, ZERZER )
1370          PTH = PXCASC * PXCASC + PYCASC * PYCASC
1371          PLA = SQRT ( PTH + PZCASC * PZCASC )
1372          PTH = SQRT ( PTH )
1373          CDE = PZCASC / PLA
1374          SDE = PTH    / PLA
1375          IF (SDE .GE. ANGLGB) THEN
1376             CFE = PXCASC / PTH
1377             SFE = PYCASC / PTH
1378          ELSE
1379             CFE = 1.D+00
1380             SFE = 0.D+00
1381          END IF
1382          CALL TTRANS ( TXI, TYI, TZI, CDE, SDE, SFE, CFE,
1383      &                 CXXINC, CYYINC, CZZINC )
1384          PXCASC = PLA * CXXINC
1385          PYCASC = PLA * CYYINC
1386          PZCASC = PLA * CZZINC
1387       END IF
1388       PXLEFT = PXTTOT - PXCASC - TXI * TVGRE0
1389       PYLEFT = PYTTOT - PYCASC - TYI * TVGRE0
1390       PZLEFT = PZTTOT - PZCASC - TZI * TVGRE0
1391       PLA = SQRT ( PXLEFT * PXLEFT + PYLEFT * PYLEFT + PZLEFT
1392      &           * PZLEFT )
1393       DEKVSP = PLA - EKIN - AM (IJ)
1394 *  +-------------------------------------------------------------------*
1395 *  |  Check the momentum versus the total energy
1396       IF ( DEKVSP .GE. 0.D+00 ) THEN
1397 *  |  +----------------------------------------------------------------*
1398 *  |  |  There are good reasons to believe that the following attempt
1399 *  |  |  is dangerous rather than useful, leading to a Dp > DEk
1400          IF ( TVGRE0 .GT. 0.D+00 ) THEN
1401             TVTENT = TVTENT - TVGRE0
1402             EKIN   = EKIN + TVGRE0
1403             TVGRE0 = 0.D+00
1404             PXLEFT = PXTTOT - PXCASC
1405             PYLEFT = PYTTOT - PYCASC
1406             PZLEFT = PZTTOT - PZCASC
1407             PLA = SQRT ( PXLEFT * PXLEFT + PYLEFT * PYLEFT + PZLEFT
1408      &                 * PZLEFT )
1409             IF ( EKIN + AM (IJ) .GT. PLA ) GO TO 300
1410 *  |  |  This part is new (1-10-91)
1411             DEEXTR = EKE - EKIN - EINTR + IGREYP * AM (1)
1412      &             + IGREYN * AM (8) - 1.5D+00
1413      &             * (IGREYP+IGREYN) * EBNDAV
1414             IF ( DEEXTR .GT. 0.D+00 ) THEN
1415                EKIN  = EKIN  + 0.5D+00 * DEEXTR
1416                EINCP = EINCP - 0.5D+00 * DEEXTR
1417                EINCN = EINCN - 0.5D+00 * DEEXTR
1418                PXLEFT = PXTTOT - PXCASC + 0.5D+00 * TXI * DEEXTR
1419                PYLEFT = PYTTOT - PYCASC + 0.5D+00 * TYI * DEEXTR
1420                PZLEFT = PZTTOT - PZCASC + 0.5D+00 * TZI * DEEXTR
1421                PLA = SQRT ( PXLEFT * PXLEFT + PYLEFT * PYLEFT + PZLEFT
1422      &                    * PZLEFT )
1423                IF ( EKIN + AM (IJ) .GT. PLA ) GO TO 300
1424             END IF
1425          END IF
1426 *  |  |
1427 *  |  +----------------------------------------------------------------*
1428          IF ( PLA - EKIN - AM (IJ) .LE. 1.D-04 * PLA ) GO TO 300
1429 *  |  +----------------------------------------------------------------*
1430 *  |  |  Printing condition relaxed, A.F., 23-12-92
1431          IF ( PLA - EKIN - AM (IJ) .GT. 1.D-02 * PLA ) THEN
1432             WRITE ( LUNERR,* )' Eventv: ekin+am < pla,ij,igreyt',
1433      &                          EKIN+AM(IJ),PLA,IJ,IGREYT
1434          END IF
1435 *  |  |
1436 *  |  +----------------------------------------------------------------*
1437          PTH = PLA
1438          PLA = EKIN + AM(IJ)
1439          PXLEFT = PXLEFT * PLA / PTH
1440          PYLEFT = PYLEFT * PLA / PTH
1441          PZLEFT = PZLEFT * PLA / PTH
1442       END IF
1443 *  |
1444 *  +-------------------------------------------------------------------*
1445   300 CONTINUE
1446 *  end new version
1447 *  +-------------------------------------------------------------------*
1448 *  |  Check if we have enough energy to enter the high energy module
1449 *  |  if not reset the accumulators and go to nucrin
1450       IF ( PLA .LT. PTHRSH (IJIJ) ) THEN
1451          PREF = 0.45D+00 * POO
1452 *  |  +----------------------------------------------------------------*
1453 *  |  |  Check if the momentum is much smaller than the original one
1454          IF ( PLA .LE. PREF ) THEN
1455 *  |  |  +-------------------------------------------------------------*
1456 *  |  |  |
1457             IF ( PLA .LE. 0.4D+00 * POO ) THEN
1458                WRITE ( LUNERR,* )' Eventv: Pla < 0.4 Poo',PLA,POO,IJIJ,
1459      &                             TKIN,IMAT
1460                LRESMP = .TRUE.
1461                GO TO 50
1462             END IF
1463 *  |  |  |
1464 *  |  |  +-------------------------------------------------------------*
1465          END IF
1466 *  |  |
1467 *  |  +----------------------------------------------------------------*
1468          PREF = MAX ( PREF, TWOTWO )
1469          PREF = MIN ( PREF, FOUFOU )
1470 *  |  +----------------------------------------------------------------*
1471 *  |  |  Originally it was < 10, but with Fermi momentum sometimes
1472 *  |  |  we were calling Hadrin with p > 10 GeV/c
1473 *        IF ( POO .LE. 9.75D+00 ) THEN
1474          IF ( ( POO .LE. 9.75D+00 .AND. IJNUCR (IJIJ) .LE. 0 )
1475      &        .OR. PLA .LT. PREF ) THEN
1476             IGREYP = 0
1477             IGREYN = 0
1478             KTARP  = 0
1479             KTARN  = 0
1480             TVGRE0 = 0.D+00
1481             TVGREY = 0.D+00
1482             EINTR  = 0.D+00
1483             PXINTR = 0.D+00
1484             PYINTR = 0.D+00
1485             PZINTR = 0.D+00
1486             ICINTR = 0
1487             IBINTR = 0
1488             ANOW   = BBTAR
1489             ZNOW   = ZZTAR
1490             IF ( IJ .EQ. 26 ) IJ = 23
1491             GO TO 100
1492          END IF
1493 *  |  |
1494 *  |  +----------------------------------------------------------------*
1495       END IF
1496 *  |
1497 *  +-------------------------------------------------------------------*
1498       TXX = PXLEFT / PLA
1499       TYY = PYLEFT / PLA
1500       TZZ = PZLEFT / PLA
1501       ZTEMP = ZZTAR - IGREYP
1502       ATEMP = BBTAR - IGREYN
1503       ERES  = ETTOT - EKIN - AM (IJ) - EINTR
1504       AMNRES = AMUAMU * ATEMP + 1.D-03 * FKENER ( ATEMP, ZTEMP ) -
1505      &         ZTEMP * AMELEC + ELBNDE ( NINT (ZTEMP) )
1506 C
1507 C********************************************************************
1508 C     FOR MOMENTA ABOVE 5.0 GEV/C USE NUCEVT
1509 C********************************************************************
1510 C
1511 *
1512 *  From here the high energy model....
1513 *
1514       NNHAD = 0
1515       CALL NUCEVV ( NNHAD, IJ, PLA, EKIN, TXX, TYY, TZZ )
1516       IF ( LRESMP ) GO TO 50
1517 *--<--<--<--<--< go to resampling if something was wrong
1518 *  +-------------------------------------------------------------------*
1519 *  |
1520       IF ( LLASTN .AND. NINT ( ANOW ) .EQ. 1 ) THEN
1521 *  |  +----------------------------------------------------------------*
1522 *  |  |
1523          IF ( KTINC .EQ. 1 ) THEN
1524             IGREYP = IGREYP + 1
1525             EINCP  = EINCP  + EKINC
1526             ZNOW   = ZNOW - 1.D+00
1527 *  |  |
1528 *  |  +----------------------------------------------------------------*
1529 *  |  |
1530          ELSE
1531             IGREYN = IGREYN + 1
1532             EINCN  = EINCN  + EKINC
1533          END IF
1534 *  |  |
1535 *  |  +----------------------------------------------------------------*
1536          ANOW = ANOW - 1.D+00
1537          NP = NP + 1
1538          IF (NP .GT. MXP) GO TO 400
1539          KPART(NP)= KPTOIP (KTINC)
1540          WEI (NP) = WE
1541          TKI (NP) = EKINC
1542 *  |  Updating conservation counters
1543          EINTR  = EINTR  + TKI (NP) + AMINC
1544          PXINTR = PXINTR + PXXINC
1545          PYINTR = PYINTR + PYYINC
1546          PZINTR = PZINTR + PZZINC
1547          ICINTR = ICINTR + ICH  (KTINC)
1548          IBINTR = IBINTR + IBAR (KTINC)
1549          PLR (NP) = SQRT (PXXINC * PXXINC + PYYINC * PYYINC +
1550      &                    PZZINC * PZZINC)
1551          CXR (NP) = PXXINC / PLR (NP)
1552          CYR (NP) = PYYINC / PLR (NP)
1553          CZR (NP) = PZZINC / PLR (NP)
1554       END IF
1555 *  |
1556 *  +-------------------------------------------------------------------*
1557
1558 *  +-------------------------------------------------------------------*
1559 *  |                          Loop over particles produced in nucevt
1560       DO 301 I=1,NNHAD
1561          NP = NP+1
1562          IF (NP .GT. MXP) GO TO 400
1563 *  |  Get the "paprop" index for Nucevt produced particles
1564          KPART(NP) = KPTOIP (NRENU(I))
1565          IF ( KPART (NP) .EQ. 0 ) THEN
1566             WRITE (LUNERR,*)' **** Charmed particle produced in Nucevv',
1567      &      NRENU(I),HEPNU(I),AMNU(I),' ****'
1568             KPART (NP) = NRENU (I)
1569          END IF
1570          PLR  (NP) = SQRT ( PXNU (I)**2 + PYNU (I)**2 + PZNU (I)**2 )
1571          CXR  (NP) = PXNU (I) / PLR (NP)
1572          CYR  (NP) = PYNU (I) / PLR (NP)
1573          CZR  (NP) = PZNU (I) / PLR (NP)
1574          TKI(NP) = HEPNU(I) - AMNU(I)
1575          WEI(NP) = WE
1576 *  |  +----------------------------------------------------------------*
1577 *  |  |
1578          IF (TKI(NP) .LE. 0.D+00) THEN
1579 *  |  | Is this the only check on the generated particle energy??
1580             EUZ = EUZ - HEPNU(I)
1581             PUX = PUX - PXNU(I)
1582             PUY = PUY - PYNU(I)
1583             PUZ = PUZ - PZNU(I)
1584             ICU = ICU - ICHNU(I)
1585             IBU = IBU - IBARNU(I)
1586             WRITE (LUNERR,*) ' Eventq: Kin en. < 0 from nucevt',TKI(NP)
1587             NP  = NP - 1
1588          END IF
1589 *  |  |
1590 *  |  +----------------------------------------------------------------*
1591  301  CONTINUE
1592 *  |
1593 *  +-------------------------------------------------------------------*
1594  500  CONTINUE
1595 *
1596 *  Now we have to compute the excitation energy: we have used Eincp and
1597 *  Eincn for cascade protons and neutrons, Tvgre0 and Tvgrey have been
1598 *  generated by cascade nucleons under threshold, Tveuz has been gene-
1599 *  rated by the high energy collisions: however Tveuz and Tvgrey are
1600 *  only approximate since they have been computed starting from an
1601 *  average binding energy, furthermore Tvgrey is already accounted for
1602 *  inside Eincp and Eincn. The actual energy spent inside high energy
1603 *  collisions was Eke - Eincp - Eincn - Tvgre0 + Efrm, the one in ca-
1604 *  scade nucleons (approximately) Eincp + Eincn - Tvgrey: so first
1605 *  check the energy balance without excitation energy
1606 *  and then compute Tv!!!
1607 *
1608 *  Now the balance!!!!!!!!!!!!!!!
1609 *
1610       ERES  = ETTOT  - EUZ - ENUCR  - EINTR
1611       PXRES = PXTTOT - PUX - PXNUCR - PXINTR
1612       PYRES = PYTTOT - PUY - PYNUCR - PYINTR
1613       PZRES = PZTTOT - PUZ - PZNUCR - PZINTR
1614       ICRES = ICHTOT - ICU - ICNUCR - ICINTR
1615       IBRES = IBTOT  - IBU - IBNUCR - IBINTR
1616 *  +-------------------------------------------------------------------*
1617 *  |
1618       IF ( NINT (ZNOW) .NE. ICRES .OR. NINT (ANOW) .NE. IBRES ) THEN
1619 *  |  +----------------------------------------------------------------*
1620 *  |  |
1621          IF ( LNUCRI ) THEN
1622             WRITE (LUNERR,*)' Eventq: charge/baryon conservation',
1623      &                      ' failure with Nucrin',
1624      &                      NINT (ZNOW), ICRES, NINT (ANOW), IBRES
1625 *  |  |
1626 *  |  +----------------------------------------------------------------*
1627 *  |  |
1628          ELSE
1629             WRITE (LUNERR,*)' Eventq: charge/baryon conservation',
1630      &                      ' failure with Nucevt',
1631      &                      NINT (ZNOW), ICRES, NINT (ANOW), IBRES
1632          END IF
1633 *  |  |
1634 *  |  +----------------------------------------------------------------*
1635          LRESMP = .TRUE.
1636          GO TO 50
1637 *  |--<--<--<--<--< go to resampling
1638       END IF
1639 *  |
1640 *  +-------------------------------------------------------------------*
1641 *  +-------------------------------------------------------------------*
1642 *  |
1643       IF ( IBRES .GT. 0 ) THEN
1644          AMMRES = ANOW * AMUAMU + 0.001D+00 * FKENER ( ANOW, ZNOW )
1645          AMNRES = AMMRES - ZNOW * AMELEC + ELBNDE ( ICRES )
1646 *        AMNRES = ANOW * AMUC12
1647          EKR0   = ERES - AMNRES
1648 *  |  Now switch to atomic masses:
1649          ERES   = ERES + AMMTAR - AMNTAR - ( ZZTAR - ZNOW ) * AMELEC
1650          EKRES  = ERES - AMMRES
1651          TVTENT = TVGRE0 + TVGREY + TVEUZ
1652          TVCOMP = ERES - ANOW * AMUC12
1653          IF ( LNUCRI ) TVCOMP = TVCOMP + ( AMNTAR - BBTAR * AMUC12 )
1654 *  |
1655 *  +-------------------------------------------------------------------*
1656 *  |
1657       ELSE
1658          AMMRES = 0.D+00
1659          AMNRES = 0.D+00
1660          ERES   = 0.D+00
1661          EKR0   = 0.D+00
1662          EKRES  = 0.D+00
1663          TVTENT = 0.D+00
1664          GO TO 600
1665       END IF
1666 *  |
1667 *  +-------------------------------------------------------------------*
1668 *  Check that the kinetic energy of the residual nuclei is not much
1669 *  different from our prevision and kinematically consistent with
1670 *  its momentum
1671 *  +-------------------------------------------------------------------*
1672 *  |
1673       IF ( EKRES .LE. 0.D+00 ) THEN
1674          WRITE ( LUNERR,* )' Eventq: negative kinetic energy for',
1675      &                     ' the residual nucleus!!',ICRES,IBRES,
1676      &                       REAL (EKRES), LNUCRI
1677 *  |  +----------------------------------------------------------------*
1678 *  |  |
1679          IF ( EKRES .LT. -3.D-3 ) THEN
1680             LRESMP = .TRUE.
1681             GO TO 50
1682 *--|--|--<--<--<--<--< go to resampling
1683          END IF
1684 *  |  |
1685 *  |  +----------------------------------------------------------------*
1686          EKRES  = 0.D+00
1687          TVRECL = 0.D+00
1688          AMSTAR = AMMRES
1689          TVCMS  = 0.D+00
1690          PTRES2 = 0.D+00
1691          PXRES  = 0.D+00
1692          PYRES  = 0.D+00
1693          PZRES  = 0.D+00
1694 *  |
1695 *  +-------------------------------------------------------------------*
1696 *  |
1697       ELSE
1698          PTRES2 = PXRES**2 + PYRES**2 + PZRES**2
1699          AMSTAR = ERES**2  - PTRES2
1700 *  |  +----------------------------------------------------------------*
1701 *  |  |
1702          IF ( AMSTAR .GE. AMMRES**2 ) THEN
1703             AMSTAR = SQRT ( AMSTAR )
1704             TVCMS  = AMSTAR - AMMRES
1705 *  |  |
1706 *  |  +----------------------------------------------------------------*
1707 *  |  |              If the following condition is satisfied it is only
1708 *  |  |              a rounding problem, set the excitation energy to 0
1709 *  |  |              and continue
1710          ELSE IF ( AMMRES**2 - AMSTAR .LT. 2.D+00 * AMSTAR * TVEPSI
1711      &             ) THEN
1712             AMSTAR = AMMRES
1713             TVCMS  = 0.D+00
1714 *  |  |
1715 *  |  +----------------------------------------------------------------*
1716 *  |  |
1717          ELSE IF ( AMSTAR .LE. 0.D+00 ) THEN
1718             WRITE ( LUNERR,* )' Eventq: immaginary mass for',
1719      &                        ' the residual nucleus!!',ICRES,IBRES,
1720      &                          REAL (AMSTAR)
1721             LRESMP = .TRUE.
1722             GO TO 50
1723 *--|--|--<--<--<--<--< go to resampling
1724 *           AMSTAR = AMMRES
1725 *           TVCMS  = 0.D+00
1726 *  |  |
1727 *  |  +----------------------------------------------------------------*
1728 *  |  |
1729          ELSE
1730             AMSTAR = SQRT ( AMSTAR )
1731 *  |  |  +-------------------------------------------------------------*
1732 *  |  |  |           If the following condition is satisfied it is only
1733 *  |  |  |           a rounding problem, set the excitation energy to 0
1734 *  |  |  |           and continue
1735             IF ( AMMRES - AMSTAR .LT. TVEPSI ) THEN
1736                AMSTAR = AMMRES
1737                TVCMS  = 0.D+00
1738                TVRECL = ERES - AMSTAR
1739                GO TO 550
1740             END IF
1741 *  |  |  |
1742 *  |  |  +-------------------------------------------------------------*
1743             WRITE ( LUNERR,* )' Eventq: negative excitation energy for',
1744      &                        ' the residual nucleus!!',ICRES,IBRES,
1745      &                          REAL ( AMSTAR - AMMRES ), LNUCRI
1746 *  |  |  +-------------------------------------------------------------*
1747 *  |  |  |
1748             IF ( AMSTAR - AMMRES .LT. -3.D-3 ) THEN
1749                LRESMP = .TRUE.
1750                GO TO 50
1751 *--|--|--|--<--<--<--<--< go to resampling
1752             END IF
1753 *  |  |  |
1754 *  |  |  +-------------------------------------------------------------*
1755             AMSTAR = AMMRES
1756             TVCMS  = 0.D+00
1757             AMSTAR = AMMRES
1758             TVCMS  = 0.D+00
1759             HELP   = SQRT ( ( ERES - AMMRES ) * ( ERES + AMMRES )
1760      &             / PTRES2 )
1761             PXRES = PXRES * HELP
1762             PYRES = PYRES * HELP
1763             PZRES = PZRES * HELP
1764             PTRES2 = PTRES2 * HELP * HELP
1765          END IF
1766 *  |  |
1767 *  |  +----------------------------------------------------------------*
1768          TVRECL = ERES - AMSTAR
1769       END IF
1770 *  |
1771 *  +-------------------------------------------------------------------*
1772   550 CONTINUE
1773 *  The following two cards are equivalent providing the kinematical
1774 *  limits are ok and we use for both amnres or ammres! Now Tv is left
1775 *  = 0
1776       TV     = 0.D+00
1777 *     TV     = EKRES
1778 *     TV     = TVRECL + TVCMS
1779 *  +-------------------------------------------------------------------*
1780 *  |
1781       IF ( .NOT. LEVPRT .AND. ABS ( ( TVTENT - TVCOMP ) / TVCOMP )
1782      &     .GT. 30000000.D+00 ) THEN
1783 *  |  +----------------------------------------------------------------*
1784 *  |  |
1785          IF ( LINCTV ) THEN
1786             WRITE ( LUNERR,* )
1787      &                      ' Eventq: excitation energy very different',
1788      &                      ' from the approximate one!!', ICRES, IBRES,
1789      &                        REAL (TVTENT), REAL (TVCOMP), LNUCRI
1790          END IF
1791 *  |  |
1792 *  |  +----------------------------------------------------------------*
1793       END IF
1794 *  |
1795 *  +-------------------------------------------------------------------*
1796       EKRES = TVRECL
1797  600  CONTINUE
1798       EOTEST = EOTEST - EUZ - ENUCR - EINTR - EKR0 - AMNRES
1799 *  +-------------------------------------------------------------------*
1800 *  |
1801       IF ( ABS (EOTEST) .GT. ETEPS ) THEN
1802 *  |  +----------------------------------------------------------------*
1803 *  |  |
1804          IF ( LNUCRI ) THEN
1805             WRITE (LUNERR,*)' Eventq: eotest failure with Nucrin',
1806      &                      EOTEST
1807 *  |  |
1808 *  |  +----------------------------------------------------------------*
1809 *  |  |
1810          ELSE
1811             WRITE (LUNERR,*)' Eventq: eotest failure with Nucevt',
1812      &                      EOTEST
1813          END IF
1814 *  |  |
1815 *  |  +----------------------------------------------------------------*
1816          LRESMP = .TRUE.
1817          GO TO 50
1818 *  |--|--<--<--<--<--< go to resampling
1819       END IF
1820 *  |
1821 *  +-------------------------------------------------------------------*
1822       IF ( IBRES .EQ. 0 ) RETURN
1823 *-->-->-->-->--> Here for the evaporation step!!!
1824 *  +-------------------------------------------------------------------*
1825 *  |                      Check if the evaporation is requested
1826       IF ( LEVPRT ) THEN
1827          PTRES = SQRT ( PTRES2 )
1828          CALL EVEVAP ( WE )
1829          IF ( LRESMP ) GO TO 50
1830 *
1831 *--|--<--<--<--<--<--< go to resampling
1832 *  |
1833 *  +-------------------------------------------------------------------*
1834 *  |                      No evaporation
1835       ELSE
1836          TVHEAV = 0.D+00
1837       END IF
1838 *  |
1839 *  +-------------------------------------------------------------------*
1840       IF (IPRI.NE.1) RETURN
1841 C
1842 C********************************************************************
1843 C     TEST PRINTOUT
1844 C********************************************************************
1845 C
1846       WRITE(LUNOUT,590)NP-NP0,NNHAD,IJ,IMAT,POO,EKE,TXI,TYI,TZI,WE
1847  590  FORMAT (4I7,6F12.6)
1848       DO 501 I=NP0+1,NP
1849       WRITE(LUNOUT,591)I,KPART(I),CXR(I),CYR(I),CZR (I),TKI(I),PLR (I),
1850      *                         WEI(I)
1851  591  FORMAT (2I5,6F12.6)
1852  501  CONTINUE
1853       RETURN
1854 C
1855 C********************************************************************
1856 C     FINUC FLOWS OVER - THIS IS FATAL - INCREASE THE SIZE OF IT
1857 C********************************************************************
1858 C
1859  400  CONTINUE
1860       WRITE(LUNOUT,490)
1861  490  FORMAT(1X,'OVERFLOW IN EVENTQ - INCREASE THE SIZE OF THE',
1862      1' COMMON BLOCK FINUC.')
1863       STOP
1864       END