Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / 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