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