GEANT321/gheisha/nucrec.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:58  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.38  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE NUCREC(NOPT,IREC)
  13 C
  14 C *** NUCLEAR REACTION KINEMATICS AT LOW ENERGIES ***
  15 C *** NVE 18-MAY-1988 CERN GENEVA ***
  16 C
  17 C CALLED BY : GHEISH, GNSLWD
  18 C ORIGIN    : H.FESEFELDT (12-FEB-1987)
  19 C
  20 C NOPT=1   N M(A,Z) --> G (G) M(A+1,Z  )    NEUTRON CAPTURE
  21 C NOPT=2   N M(A,Z) --> N (G) M(A  ,Z  )    INELASTIC NEUTRON SCATT.
  22 C NOPT=3   N M(A,Z) --> P (G) M(A  ,Z-1)
  23 C NOPT=4   N M(A,Z) --> D (G) M(A-1,Z-1)
  24 C NOPT=5   N M(A,Z) --> T (G) M(A-2,Z-1)
  25 C NOPT=6   N M(A,Z) --> ALP.  M(A-3,Z-2)
  26 C NOPT=7   N M(A,Z) --> N N   M(A-1,Z  )
  27 C NOPT=8   N M(A,Z) --> N P   M(A-1,Z-1)
  28 C NOPT=9   N M(A,Z) --> P P   M(A-1,Z-2)
  29 C NOPT=11  P M(A,Z) --> G (G) M(A+1,Z+1)    PROTON CAPTURE
  30 C NOPT=12  P M(A,Z) --> N (G) M(A  ,Z  )    INELASTIC PROTON SCATT.
  31 C NOPT=13  P M(A,Z) --> P (G) M(A  ,Z+1)
  32 C NOPT=14  P M(A,Z) --> D (G) M(A-1,Z  )
  33 C NOPT=15  P M(A,Z) --> T (G) M(A-2,Z  )
  34 C NOPT=16  P M(A,Z) --> ALP.  M(A-3,Z-1)
  35 C NOPT=17  P M(A,Z) --> N N   M(A-1,Z+1)
  36 C NOPT=18  P M(A,Z) --> N P   M(A-1,Z  )
  37 C NOPT=19  P M(A,Z) --> P P   M(A-1,Z-1)
  38 C SIMILAR FOR D,T,ALPHA SCATTERING ON NUCLEI
  39 C
  40 C NOTE : DOUBLE PRECISION CALCULATIONS ARE VITAL FOR THESE LOW
  41 C        ENERGY PROCESSES
  42 C        THEREFORE THE VARS OF /GENIO/ ARE DECLARED DOUBLE PRECISION
  43 C        ALSO A DOUBLE PRECISION VERSION OF THE PHASE SPACE PACKAGE
  44 C        "PHPNUC" HAS BEEN INTRODUCED
  45 C *** HMF 29-AUG-1989 RWTH AACHEN ***
  46 C
  47 #include "geant321/s_defcom.inc"
  48 #include "geant321/s_nucio.inc"
  49 C
  50       DIMENSION QVAL(10),TCH(10)
  51       DIMENSION RNDM(2)
  52 C
  53 C** PROGRAM RETURNS WITH NOPT=0, IF INELASTIC SCATTERING ENERGETICALLY
  54 C** NOT POSSIBLE, OR IF WRONG PARTICLES ENTER THIS ROUTINE: ONLY FOR
  55 C** PROTONS,NEUTRONS, DEUTERIUM, TRITIUM AND ALPHAS.
  56 C** IF EK > 100 MEV, THIS ROUTINE IS CERTAINLY NOT ADEQUATE.
  57 C
  58       NOPT=0
  59       IF (IREC .EQ. 0) GO TO 9999
  60 C
  61       IF (NPRT(9) .AND. (EK .GT. 0.1)) PRINT 9000,EK,IPART
  62  9000 FORMAT(' *NUCREC* ENERGY TOO HIGH EK = ',G12.5,' GEV ',
  63      $ ' KPART = ',I3)
  64       IF (EK .GT. 0.1) GO TO 9999
  65 C
  66 C%%%      IF(IPART.EQ.16) GOTO 2
  67 C%%%      IF(IPART.EQ.14) GOTO 3
  68 C%%%      IF(IPART.EQ.30) GOTO 4
  69 C%%%      IF(IPART.EQ.31) GOTO 5
  70 C%%%      IF(IPART.EQ.32) GOTO 6
  71 C%%%      GO TO 9999
  72 C%%%    2 AMAS = ATOMAS(1.,0.)
  73 C%%%      GOTO 8
  74 C%%%    3 AMAS = ATOMAS(1.,1.)
  75 C%%%      GOTO 8
  76 C%%%    4 AMAS = ATOMAS(2.,1.)
  77 C%%%      GOTO 8
  78 C%%%    5 AMAS = ATOMAS(3.,1.)
  79 C%%%      GOTO 8
  80 C%%%    6 AMAS = ATOMAS(4.,2.)
  81 C
  82       IF (IPART .EQ. 16) GO TO 8
  83       IF (IPART .EQ. 14) GO TO 8
  84       IF (IPART .EQ. 30) GO TO 8
  85       IF (IPART .EQ. 31) GO TO 8
  86       IF( IPART .EQ. 32) GO TO 8
  87       GO TO 9999
  88 C** SET BEAM PARTICLE, TAKE EK AS FUNDAMENTAL QUANTITY
  89 C** DUE TO THE DIFFICULT KINEMATIC, ALL MASSES HAVE TO BE ASSIGNED
  90 C** THE BEST MEASURED VALUES.
  91  8    CONTINUE
  92       CALL VZERO(QVAL,10)
  93       CALL VZERO(TCH ,10)
  94 C --- GET MASS WHICH MATCHES GEANT ---
  95       AMAS=RMASS(IPART)
  96       EN=EK+ABS(AMAS)
  97       P =SQRT(ABS(EN*EN-AMAS*AMAS))
  98       PP=SQRT(PX*PX+PY*PY+PZ*PZ)
  99       IF (PP .GT. 1.0E-6) GO TO 8000
 100       CALL GRNDM(RNDM,2)
 101       PHINVE=TWPI*RNDM(1)
 102       COST=-1.+2.*RNDM(2)
 103       IF (COST .LE. -1.) COST=-1.
 104       IF (COST .GE.  1.) COST= 1.
 105       RTHNVE=ACOS(COST)
 106       PX=SIN(RTHNVE)*COS(PHINVE)
 107       PY=SIN(RTHNVE)*SIN(PHINVE)
 108       PZ=COS(RTHNVE)
 109       PP=1.
 110  8000 CONTINUE
 111       PX=PX/PP
 112       PY=PY/PP
 113       PZ=PZ/PP
 114       CALL VZERO(PV,10*MXGKPV)
 115       PV(1,1) =PX*P
 116       PV(2,1) =PY*P
 117       PV(3,1) =PZ*P
 118       PV(4,1) =EN
 119       PV(5,1) =AMAS
 120       PV(6,1) =NCH
 121       PV(7,1) =TOF
 122       PV(8,1) =IPART
 123       PV(9,1) =0.
 124       PV(10,1)=USERW
 125       PV(1,2) =0.
 126       PV(2,2) =0.
 127       PV(3,2) =0.
 128       PV(4,2) =0.
 129       PV(5,2) =ATOMAS(ATNO2,ZNO2)
 130       PV(6,2) =ZNO2
 131       PV(7,2) =TOF
 132       PV(8,2) =0.
 133       PV(9,2) =0.
 134       PV(10,2)=0.
 135 C** CALCULATE Q-VALUE OF REACTIONS
 136       IF(IPART.EQ.16) GOTO 20
 137       IF(IPART.EQ.14) GOTO 30
 138       IF(IPART.EQ.30) GOTO 40
 139       IF(IPART.EQ.31) GOTO 50
 140       IF(IPART.EQ.32) GOTO 60
 141    20 PV(5,11)=ATOMAS(ATNO2+1.,ZNO2   )
 142       PV(6,11)=ZNO2
 143       PV(5,21)=0.
 144       PV(6,21)=0.
 145       PV(8,21)=1.
 146       PV(5,31)=0.
 147       PV(6,31)=0.
 148       PV(8,31)=1.
 149 C
 150       PV(5,12)=PV(5,2)
 151       PV(6,12)=PV(6,2)
 152       PV(5,22)=RMASS(16)
 153       PV(6,22)=0.
 154       PV(8,22)=16.
 155       PV(5,32)=0.
 156       PV(6,32)=0.
 157       PV(8,32)=1.
 158 C
 159       PV(5,13)=ATOMAS(ATNO2   ,ZNO2-1.)
 160       PV(6,13)=ZNO2-1.
 161       PV(5,23)=RMASS(14)
 162       PV(6,23)=1.
 163       PV(8,23)=14.
 164       PV(5,33)=0.
 165       PV(6,33)=0.
 166       PV(8,33)=1.
 167 C
 168       PV(5,14)=ATOMAS(ATNO2-1.,ZNO2-1.)
 169       PV(6,14)=ZNO2-1.
 170       PV(5,24)=RMASS(30)
 171       PV(6,24)=1.
 172       PV(8,24)=30.
 173       PV(5,34)=0.
 174       PV(6,34)=0.
 175       PV(8,34)=1.
 176 C
 177       PV(5,15)=ATOMAS(ATNO2-2.,ZNO2-1.)
 178       PV(6,15)=ZNO2-1.
 179       PV(5,25)=RMASS(31)
 180       PV(6,25)=1.
 181       PV(8,25)=31.
 182       PV(5,35)=0.
 183       PV(6,35)=0.
 184       PV(8,35)=1.
 185 C
 186       PV(5,16)=ATOMAS(ATNO2-3.,ZNO2-2.)
 187       PV(6,16)=ZNO2-2.
 188       PV(5,26)=RMASS(32)
 189       PV(6,26)=2.
 190       PV(8,26)=32.
 191       PV(5,36)=0.
 192       PV(6,36)=0.
 193       PV(8,36)=1.
 194 C
 195       PV(5,17)=ATOMAS(ATNO2-1.,ZNO2   )
 196       PV(6,17)=ZNO2
 197       PV(5,27)=PV(5,22)
 198       PV(6,27)=0.
 199       PV(8,27)=16.
 200       PV(5,37)=PV(5,22)
 201       PV(6,37)=0.
 202       PV(8,37)=16.
 203 C
 204       PV(5,18)=PV(5,14)
 205       PV(6,18)=PV(6,14)
 206       PV(5,28)=PV(5,22)
 207       PV(6,28)=0.
 208       PV(8,28)=16.
 209       PV(5,38)=PV(5,23)
 210       PV(6,38)=1.
 211       PV(8,38)=14.
 212 C
 213       PV(5,19)=ATOMAS(ATNO2-1.,ZNO2-2.)
 214       PV(6,19)=ZNO2-2.
 215       PV(5,29)=PV(5,23)
 216       PV(6,29)=1.
 217       PV(8,29)=14.
 218       PV(5,39)=PV(5,23)
 219       PV(6,39)=1.
 220       PV(8,39)=14.
 221 C
 222       GOTO 70
 223    30 PV(5,11)=ATOMAS(ATNO2+1.,ZNO2+1.)
 224       PV(6,11)=ZNO2+1.
 225       PV(5,21)=0.
 226       PV(6,21)=0.
 227       PV(8,21)=1.
 228       PV(5,31)=0.
 229       PV(6,31)=0.
 230       PV(8,31)=1.
 231 C
 232       PV(5,12)=ATOMAS(ATNO2   ,ZNO2+1.)
 233       PV(6,12)=ZNO2+1.
 234       PV(5,22)=RMASS(16)
 235       PV(6,22)=0.
 236       PV(8,22)=16.
 237       PV(5,32)=0.
 238       PV(6,32)=0.
 239       PV(8,32)=1.
 240 C
 241       PV(5,13)=PV(5,2)
 242       PV(6,13)=PV(6,2)
 243       PV(5,23)=RMASS(14)
 244       PV(6,23)=1.
 245       PV(8,23)=14.
 246       PV(5,33)=0.
 247       PV(6,33)=0.
 248       PV(8,33)=1.
 249 C
 250       PV(5,14)=ATOMAS(ATNO2-1.,ZNO2   )
 251       PV(6,14)=ZNO2
 252       PV(5,24)=RMASS(30)
 253       PV(6,24)=1.
 254       PV(8,24)=30.
 255       PV(5,34)=0.
 256       PV(6,34)=0.
 257       PV(8,34)=1.
 258 C
 259       PV(5,15)=ATOMAS(ATNO2-2.,ZNO2   )
 260       PV(6,15)=ZNO2
 261       PV(5,25)=RMASS(31)
 262       PV(6,25)=1.
 263       PV(8,25)=31.
 264       PV(5,35)=0.
 265       PV(6,35)=0.
 266       PV(8,35)=1.
 267 C
 268       PV(5,16)=ATOMAS(ATNO2-3.,ZNO2-1.)
 269       PV(6,16)=ZNO2-1.
 270       PV(5,26)=RMASS(32)
 271       PV(6,26)=2.
 272       PV(8,26)=32.
 273       PV(5,36)=0.
 274       PV(6,36)=0.
 275       PV(8,36)=1.
 276 C
 277       PV(5,17)=ATOMAS(ATNO2-1.,ZNO2+1.)
 278       PV(6,17)=ZNO2+1.
 279       PV(5,27)=PV(5,22)
 280       PV(6,27)=0.
 281       PV(8,27)=16.
 282       PV(5,37)=PV(5,22)
 283       PV(6,37)=0.
 284       PV(8,37)=16.
 285 C
 286       PV(5,18)=PV(5,14)
 287       PV(6,18)=PV(6,14)
 288       PV(5,28)=PV(5,22)
 289       PV(6,28)=0.
 290       PV(8,28)=16.
 291       PV(5,38)=PV(5,23)
 292       PV(6,38)=1.
 293       PV(8,38)=14.
 294 C
 295       PV(5,19)=ATOMAS(ATNO2-1.,ZNO2-1.)
 296       PV(6,19)=ZNO2-1.
 297       PV(5,29)=PV(5,23)
 298       PV(6,29)=1.
 299       PV(8,29)=14.
 300       PV(5,39)=PV(5,23)
 301       PV(6,39)=1.
 302       PV(8,39)=14.
 303 C
 304       NOPT=10
 305       GOTO 70
 306    40 PV(5,11)=ATOMAS(ATNO2+2.,ZNO2+1.)
 307       PV(6,11)=ZNO2+1.
 308       PV(5,21)=0.
 309       PV(6,21)=0.
 310       PV(8,21)=1.
 311       PV(5,31)=0.
 312       PV(6,31)=0.
 313       PV(8,31)=1.
 314 C
 315       PV(5,12)=ATOMAS(ATNO2+1.,ZNO2+1.)
 316       PV(6,12)=ZNO2+1.
 317       PV(5,22)=RMASS(16)
 318       PV(6,22)=0.
 319       PV(8,22)=16.
 320       PV(5,32)=0.
 321       PV(6,32)=0.
 322       PV(8,32)=1.
 323 C
 324       PV(5,13)=ATOMAS(ATNO2+1.,ZNO2   )
 325       PV(6,13)=ZNO2
 326       PV(5,23)=RMASS(14)
 327       PV(6,23)=1.
 328       PV(8,23)=14.
 329       PV(5,33)=0.
 330       PV(6,33)=0.
 331       PV(8,33)=1.
 332 C
 333       PV(5,14)=PV(5,2)
 334       PV(6,14)=PV(6,2)
 335       PV(5,24)=RMASS(30)
 336       PV(6,24)=1.
 337       PV(8,24)=30.
 338       PV(5,34)=0.
 339       PV(6,34)=0.
 340       PV(8,34)=1.
 341 C
 342       PV(5,15)=ATOMAS(ATNO2-1.,ZNO2   )
 343       PV(6,15)=ZNO2
 344       PV(5,25)=RMASS(31)
 345       PV(6,25)=1.
 346       PV(8,25)=31.
 347       PV(5,35)=0.
 348       PV(6,35)=0.
 349       PV(8,35)=1.
 350 C
 351       PV(5,16)=ATOMAS(ATNO2-2.,ZNO2-1.)
 352       PV(6,16)=ZNO2-1.
 353       PV(5,26)=RMASS(32)
 354       PV(6,26)=2.
 355       PV(8,26)=32.
 356       PV(5,36)=0.
 357       PV(6,36)=0.
 358       PV(8,36)=1.
 359 C
 360       PV(5,17)=ATOMAS(ATNO2   ,ZNO2+1.)
 361       PV(6,17)=ZNO2+1.
 362       PV(5,27)=PV(5,22)
 363       PV(6,27)=0.
 364       PV(8,27)=16.
 365       PV(5,37)=PV(5,22)
 366       PV(6,37)=0.
 367       PV(8,37)=16.
 368 C
 369       PV(5,18)=PV(5,14)
 370       PV(6,18)=PV(6,14)
 371       PV(5,28)=PV(5,22)
 372       PV(6,28)=0.
 373       PV(8,28)=16.
 374       PV(5,38)=PV(5,23)
 375       PV(6,38)=1.
 376       PV(8,38)=14.
 377 C
 378       PV(5,19)=ATOMAS(ATNO2   ,ZNO2-1.)
 379       PV(6,19)=ZNO2-1.
 380       PV(5,29)=PV(5,23)
 381       PV(6,29)=1.
 382       PV(8,29)=14.
 383       PV(5,39)=PV(5,23)
 384       PV(6,39)=1.
 385       PV(8,39)=14.
 386 C
 387       NOPT=20
 388       GOTO 70
 389    50 PV(5,11)=ATOMAS(ATNO2+3.,ZNO2+1.)
 390       PV(6,11)=ZNO2+1.
 391       PV(5,21)=0.
 392       PV(6,21)=0.
 393       PV(8,21)=1.
 394       PV(5,31)=0.
 395       PV(6,31)=0.
 396       PV(8,31)=1.
 397 C
 398       PV(5,12)=ATOMAS(ATNO2+2.,ZNO2+1.)
 399       PV(6,12)=ZNO2+1.
 400       PV(5,22)=RMASS(16)
 401       PV(6,22)=0.
 402       PV(8,22)=16.
 403       PV(5,32)=0.
 404       PV(6,32)=0.
 405       PV(8,32)=1.
 406 C
 407       PV(5,13)=ATOMAS(ATNO2+2.,ZNO2   )
 408       PV(6,13)=ZNO2
 409       PV(5,23)=RMASS(14)
 410       PV(6,23)=1.
 411       PV(8,23)=14.
 412       PV(5,33)=0.
 413       PV(6,33)=0.
 414       PV(8,33)=1.
 415 C
 416       PV(5,14)=ATOMAS(ATNO2+1.,ZNO2   )
 417       PV(6,14)=ZNO2
 418       PV(5,24)=RMASS(30)
 419       PV(6,24)=1.
 420       PV(8,24)=30.
 421       PV(5,34)=0.
 422       PV(6,34)=0.
 423       PV(8,34)=1.
 424 C
 425       PV(5,15)=PV(5,2)
 426       PV(6,15)=PV(6,2)
 427       PV(5,25)=RMASS(31)
 428       PV(6,25)=1.
 429       PV(8,25)=31.
 430       PV(5,35)=0.
 431       PV(6,35)=0.
 432       PV(8,35)=1.
 433 C
 434       PV(5,16)=ATOMAS(ATNO2-1.,ZNO2-1.)
 435       PV(6,16)=ZNO2-1.
 436       PV(5,26)=RMASS(32)
 437       PV(6,26)=2.
 438       PV(8,26)=32.
 439       PV(5,36)=0.
 440       PV(6,36)=0.
 441       PV(8,36)=1.
 442 C
 443       PV(5,17)=ATOMAS(ATNO2+1.,ZNO2+1.)
 444       PV(6,17)=ZNO2+1.
 445       PV(5,27)=PV(5,22)
 446       PV(6,27)=0.
 447       PV(8,27)=16.
 448       PV(5,37)=PV(5,22)
 449       PV(6,37)=0.
 450       PV(8,37)=16.
 451 C
 452       PV(5,18)=PV(5,14)
 453       PV(6,18)=PV(6,14)
 454       PV(5,28)=PV(5,22)
 455       PV(6,28)=0.
 456       PV(8,28)=16.
 457       PV(5,38)=PV(5,23)
 458       PV(6,38)=1.
 459       PV(8,38)=14.
 460 C
 461       PV(5,19)=ATOMAS(ATNO2+1.,ZNO2-1.)
 462       PV(6,19)=ZNO2-1.
 463       PV(5,29)=PV(5,23)
 464       PV(6,29)=1.
 465       PV(8,29)=14.
 466       PV(5,39)=PV(5,23)
 467       PV(6,39)=1.
 468       PV(8,39)=14.
 469 C
 470       NOPT=30
 471       GOTO 70
 472    60 PV(5,11)=ATOMAS(ATNO2+4.,ZNO2+2.)
 473       PV(6,11)=ZNO2+2.
 474       PV(5,21)=0.
 475       PV(6,21)=0.
 476       PV(8,21)=1.
 477       PV(5,31)=0.
 478       PV(6,31)=0.
 479       PV(8,31)=1.
 480 C
 481       PV(5,12)=ATOMAS(ATNO2+3.,ZNO2+2.)
 482       PV(6,12)=ZNO2+2.
 483       PV(5,22)=RMASS(16)
 484       PV(6,22)=0.
 485       PV(8,22)=16.
 486       PV(5,32)=0.
 487       PV(6,32)=0.
 488       PV(8,32)=1.
 489 C
 490       PV(5,13)=ATOMAS(ATNO2+3.,ZNO2+1.)
 491       PV(6,13)=ZNO2+1.
 492       PV(5,23)=RMASS(14)
 493       PV(6,23)=1.
 494       PV(8,23)=14.
 495       PV(5,33)=0.
 496       PV(6,33)=0.
 497       PV(8,33)=1.
 498 C
 499       PV(5,14)=ATOMAS(ATNO2+2.,ZNO2+1.)
 500       PV(6,14)=ZNO2+1.
 501       PV(5,24)=RMASS(30)
 502       PV(6,24)=1.
 503       PV(8,24)=30.
 504       PV(5,34)=0.
 505       PV(6,34)=0.
 506       PV(8,34)=1.
 507 C
 508       PV(5,15)=ATOMAS(ATNO2+1.,ZNO2+1.)
 509       PV(6,15)=ZNO2+1.
 510       PV(5,25)=RMASS(31)
 511       PV(6,25)=1.
 512       PV(8,25)=31.
 513       PV(5,35)=0.
 514       PV(6,35)=0.
 515       PV(8,35)=1.
 516 C
 517       PV(5,16)=PV(5,2)
 518       PV(6,16)=PV(6,2)
 519       PV(5,26)=RMASS(32)
 520       PV(6,26)=2.
 521       PV(8,26)=32.
 522       PV(5,36)=0.
 523       PV(6,36)=0.
 524       PV(8,36)=1.
 525 C
 526       PV(5,17)=ATOMAS(ATNO2+2.,ZNO2+2.)
 527       PV(6,17)=ZNO2+2.
 528       PV(5,27)=PV(5,22)
 529       PV(6,27)=0.
 530       PV(8,27)=16.
 531       PV(5,37)=PV(5,22)
 532       PV(6,37)=0.
 533       PV(8,37)=16.
 534 C
 535       PV(5,18)=PV(5,14)
 536       PV(6,18)=PV(6,14)
 537       PV(5,28)=PV(5,22)
 538       PV(6,28)=0.
 539       PV(8,28)=16.
 540       PV(5,38)=PV(5,23)
 541       PV(6,38)=1.
 542       PV(8,38)=14.
 543 C
 544       PV(5,19)=ATOMAS(ATNO2+2.,ZNO2   )
 545       PV(6,19)=ZNO2
 546       PV(5,29)=PV(5,23)
 547       PV(6,29)=1.
 548       PV(8,29)=14.
 549       PV(5,39)=PV(5,23)
 550       PV(6,39)=1.
 551       PV(8,39)=14.
 552 C
 553       NOPT=40
 554    70 QV     =EK+PV(5,2)+PV(5,1)
 555       TC     =   PV(6,2)+PV(6,1)
 556       QVAL(1)=QV - PV(5,11)
 557       TCH (1)=TC - PV(6,11)
 558       QVAL(2)=QV - PV(5,12) - PV(5,22)
 559       TCH (2)=TC - PV(6,12) - PV(6,22)
 560       QVAL(3)=QV - PV(5,13) - PV(5,23)
 561       TCH (3)=TC - PV(6,13) - PV(6,23)
 562       QVAL(4)=QV - PV(5,14) - PV(5,24)
 563       TCH (4)=TC - PV(6,14) - PV(6,24)
 564       QVAL(5)=QV - PV(5,15) - PV(5,25)
 565       TCH (5)=TC - PV(6,15) - PV(6,25)
 566       QVAL(6)=QV - PV(5,16) - PV(5,26)
 567       TCH (6)=TC - PV(6,16) - PV(6,26)
 568       QVAL(7)=QV - PV(5,17) - PV(5,27) - PV(5,37)
 569       TCH (7)=TC - PV(6,17) - PV(6,27) - PV(6,37)
 570       QVAL(8)=QV - PV(5,18) - PV(5,28) - PV(5,38)
 571       TCH (8)=TC - PV(6,18) - PV(6,28) - PV(6,38)
 572       QVAL(9)=QV - PV(5,19) - PV(5,29) - PV(5,39)
 573       TCH (9)=TC - PV(6,19) - PV(6,29) - PV(6,39)
 574    74 QV = 0
 575       IF(IREC.EQ.2) QVAL(1)=0.
 576       IF(IPART.NE.16) GOTO 75
 577       CALL GRNDM(RNDM,2)
 578       IF(RNDM(1).GT.((ATNO2-1.)/230.)**2) QVAL(1)=0.
 579       EKA=7.9254/ATNO2
 580       IF(RNDM(2).LT.EK/EKA) GOTO 75
 581       QVAL(3)=0.
 582       QVAL(4)=0.
 583       QVAL(5)=0.
 584       QVAL(6)=0.
 585       QVAL(9)=0.
 586    75 DO 71 I=1,9
 587       IF(PV(5,10+I).LT.0.5) QVAL(I)=0.
 588       IF(QVAL(I).LT.0.    ) QVAL(I)=0.
 589       IF(ABS(TCH(I)-0.1).GT.0.5 ) QVAL(I)=0.
 590       QV=QV+QVAL(I)
 591    71 CONTINUE
 592       CALL GRNDM(RNDM,1)
 593       RAN=RNDM(1)
 594       QV1=0.
 595       DO 72 I=1,9
 596       IF(QVAL(I).EQ.0.) GOTO 72
 597       QV1=QV1+QVAL(I)/QV
 598       IF(RAN.LE.QV1) GOTO 73
 599    72 CONTINUE
 600 C** REACTION KINEMATICALLY NOT POSSIBLE
 601       NOPT=0
 602       GO TO 9999
 603    73 NOPT=NOPT+I
 604       PV(5,3)=PV(5,10+I)
 605       PV(6,3)=PV(6,10+I)
 606       PV(8,3)=0.
 607       PV(5,4)=PV(5,20+I)
 608       PV(6,4)=PV(6,20+I)
 609       PV(8,4)=PV(8,20+I)
 610       PV(5,5)=PV(5,30+I)
 611       PV(6,5)=PV(6,30+I)
 612       PV(8,5)=PV(8,30+I)
 613       NT=2
 614       RAN=EK*10.
 615       IF(RAN.GT.0.5) RAN=0.5
 616       CALL GRNDM(RNDM,1)
 617       IF(RNDM(1).LT.RAN) NT=3
 618       IF(MOD(NOPT,10).GE.7) NT=3
 619 C** CALCULATE CMS ENERGY
 620    80 PV(4,2)=PV(5,2)
 621       CALL ADD(1,2,MXGKPV)
 622       PV(1,MXGKPV)=-PV(1,MXGKPV)
 623       PV(2,MXGKPV)=-PV(2,MXGKPV)
 624       PV(3,MXGKPV)=-PV(3,MXGKPV)
 625 C** SET QUANTITIES FOR PHASE SPACE ROUTINE IN CMS
 626       TECM=PV(5,MXGKPV)
 627       NPG=NT
 628       KGENEV=1
 629       DO 81 I=1,NPG
 630    81 AMASS(I)=PV(5,2+I)
 631 C --- Invoke double precision version of the phase space package ---
 632       CALL PHPNUC
 633       DO 83 I=1,NPG
 634       DO 82 J=1,4
 635    82 PV(J,2+I)=PCM(J,I)
 636 C** TRANSFORM INTO LAB.SYSTEM
 637       CALL LOR(2+I,MXGKPV,2+I)
 638       PV(7,2+I)=TOF
 639    83 CONTINUE
 640 C** SET CHARGES AND PARTICLE INDEX FOR LOW MASS FRAGMENTS
 641       IF (ABS(PV(5,3)-RMASS(14)) .LT. 0.0001) GO TO 84
 642       IF (ABS(PV(5,3)-RMASS(16)) .LT. 0.0001) GO TO 85
 643       IF (ABS(PV(5,3)-RMASS(30)) .LT. 0.0001) GO TO 86
 644       IF (ABS(PV(5,3)-RMASS(31)) .LT. 0.0001) GO TO 87
 645       IF (ABS(PV(5,3)-RMASS(32)) .LT. 0.0001) GO TO 88
 646       GOTO 89
 647    84 PV(6,3)=1.
 648       PV(8,3)=14.
 649       GOTO 89
 650    85 PV(6,3)=0.
 651       PV(8,3)=16.
 652       GOTO 89
 653    86 PV(6,3)=1.
 654       PV(8,3)=30.
 655       GOTO 89
 656    87 PV(6,3)=1.
 657       PV(8,3)=31.
 658       GOTO 89
 659    88 PV(6,3)=2.
 660       PV(8,3)=32.
 661    89 NTT=2+NT
 662       DO 90 I=1,NTT
 663       IPP=IFIX(PV(8,I)+0.01)
 664       IF(IPP.EQ.0) GOTO 90
 665       EK=PV(4,I)-PV(5,I)
 666       IF(I.LT.3) GOTO 92
 667       IF(IPP.LT.30) GOTO 92
 668       CALL GRNDM(RNDM,1)
 669       EK=EK*0.5*RNDM(1)
 670    92 IF(EK.LT.1.E-6) EK=1.E-6
 671       PV(5,I)=RMASS(IPP)
 672       PV(4,I)=EK+PV(5,I)
 673       P=SQRT(ABS(PV(4,I)**2-PV(5,I)**2))
 674       PP=SQRT(PV(1,I)**2+PV(2,I)**2+PV(3,I)**2)
 675       IF(PP.GT.1.E-6) GOTO 91
 676       CALL GRNDM(RNDM,2)
 677       PHINVE=TWPI*RNDM(1)
 678       COST=-1.+2.*RNDM(2)
 679       IF (COST .LE. -1.) COST=-1.
 680       IF (COST .GE.  1.) COST= 1.
 681       RTHNVE=ACOS(COST)
 682       PV(1,I)=SIN(RTHNVE)*COS(PHINVE)
 683       PV(2,I)=SIN(RTHNVE)*SIN(PHINVE)
 684       PV(3,I)=COS(RTHNVE)
 685       PP=1.
 686    91 PV(1,I)=PV(1,I)*P/PP
 687       PV(2,I)=PV(2,I)*P/PP
 688       PV(3,I)=PV(3,I)*P/PP
 689    90 CONTINUE
 690       IF(.NOT.NPRT(4)) GOTO 100
 691       WRITE(NEWBCD,1000) XEND,YEND,ZEND,IND,NOPT
 692  1000 FORMAT(1H ,'Nuclear reaction at (X,Y,Z) ',3(G12.5,1X)
 693      +,' Material ',I5,' NOPT ',I5)
 694       DO 95 I=1,NTT
 695          WRITE(NEWBCD,1001) I,(PV(J,I),J=1,10)
 696    95 CONTINUE
 697  1001 FORMAT(1H ,I3,1X,10(G10.3,1X))
 698   100 INTCT=INTCT+1.
 699 C** SET INTERACTION MODE ACCORDING TO GHEISHA-CONVENTION
 700 C** N-CAPTURE
 701       IF(PV(8,3).GT.0.) GOTO 110
 702       CALL SETCUR(4)
 703       NTK=NTK+1
 704       IF(NT.EQ.3) CALL SETTRK(5)
 705       GO TO 9999
 706  110  CONTINUE
 707       CALL SETCUR(4)
 708       NTK=NTK+1
 709       CALL SETTRK(3)
 710       IF(NT.EQ.3) CALL SETTRK(5)
 711       CALL SETTRK(3)
 712 C
 713  9999 CONTINUE
 714       END