This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gheisha / nucrec.F
CommitLineData
fe4da5cc 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)
13C
14C *** NUCLEAR REACTION KINEMATICS AT LOW ENERGIES ***
15C *** NVE 18-MAY-1988 CERN GENEVA ***
16C
17C CALLED BY : GHEISH, GNSLWD
18C ORIGIN : H.FESEFELDT (12-FEB-1987)
19C
20C NOPT=1 N M(A,Z) --> G (G) M(A+1,Z ) NEUTRON CAPTURE
21C NOPT=2 N M(A,Z) --> N (G) M(A ,Z ) INELASTIC NEUTRON SCATT.
22C NOPT=3 N M(A,Z) --> P (G) M(A ,Z-1)
23C NOPT=4 N M(A,Z) --> D (G) M(A-1,Z-1)
24C NOPT=5 N M(A,Z) --> T (G) M(A-2,Z-1)
25C NOPT=6 N M(A,Z) --> ALP. M(A-3,Z-2)
26C NOPT=7 N M(A,Z) --> N N M(A-1,Z )
27C NOPT=8 N M(A,Z) --> N P M(A-1,Z-1)
28C NOPT=9 N M(A,Z) --> P P M(A-1,Z-2)
29C NOPT=11 P M(A,Z) --> G (G) M(A+1,Z+1) PROTON CAPTURE
30C NOPT=12 P M(A,Z) --> N (G) M(A ,Z ) INELASTIC PROTON SCATT.
31C NOPT=13 P M(A,Z) --> P (G) M(A ,Z+1)
32C NOPT=14 P M(A,Z) --> D (G) M(A-1,Z )
33C NOPT=15 P M(A,Z) --> T (G) M(A-2,Z )
34C NOPT=16 P M(A,Z) --> ALP. M(A-3,Z-1)
35C NOPT=17 P M(A,Z) --> N N M(A-1,Z+1)
36C NOPT=18 P M(A,Z) --> N P M(A-1,Z )
37C NOPT=19 P M(A,Z) --> P P M(A-1,Z-1)
38C SIMILAR FOR D,T,ALPHA SCATTERING ON NUCLEI
39C
40C NOTE : DOUBLE PRECISION CALCULATIONS ARE VITAL FOR THESE LOW
41C ENERGY PROCESSES
42C THEREFORE THE VARS OF /GENIO/ ARE DECLARED DOUBLE PRECISION
43C ALSO A DOUBLE PRECISION VERSION OF THE PHASE SPACE PACKAGE
44C "PHPNUC" HAS BEEN INTRODUCED
45C *** HMF 29-AUG-1989 RWTH AACHEN ***
46C
47#include "geant321/s_defcom.inc"
48#include "geant321/s_nucio.inc"
49C
50 DIMENSION QVAL(10),TCH(10)
51 DIMENSION RNDM(2)
52C
53C** PROGRAM RETURNS WITH NOPT=0, IF INELASTIC SCATTERING ENERGETICALLY
54C** NOT POSSIBLE, OR IF WRONG PARTICLES ENTER THIS ROUTINE: ONLY FOR
55C** PROTONS,NEUTRONS, DEUTERIUM, TRITIUM AND ALPHAS.
56C** IF EK > 100 MEV, THIS ROUTINE IS CERTAINLY NOT ADEQUATE.
57C
58 NOPT=0
59 IF (IREC .EQ. 0) GO TO 9999
60C
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
65C
66C%%% IF(IPART.EQ.16) GOTO 2
67C%%% IF(IPART.EQ.14) GOTO 3
68C%%% IF(IPART.EQ.30) GOTO 4
69C%%% IF(IPART.EQ.31) GOTO 5
70C%%% IF(IPART.EQ.32) GOTO 6
71C%%% GO TO 9999
72C%%% 2 AMAS = ATOMAS(1.,0.)
73C%%% GOTO 8
74C%%% 3 AMAS = ATOMAS(1.,1.)
75C%%% GOTO 8
76C%%% 4 AMAS = ATOMAS(2.,1.)
77C%%% GOTO 8
78C%%% 5 AMAS = ATOMAS(3.,1.)
79C%%% GOTO 8
80C%%% 6 AMAS = ATOMAS(4.,2.)
81C
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
88C** SET BEAM PARTICLE, TAKE EK AS FUNDAMENTAL QUANTITY
89C** DUE TO THE DIFFICULT KINEMATIC, ALL MASSES HAVE TO BE ASSIGNED
90C** THE BEST MEASURED VALUES.
91 8 CONTINUE
92 CALL VZERO(QVAL,10)
93 CALL VZERO(TCH ,10)
94C --- 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.
135C** 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.
149C
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.
158C
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.
167C
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.
176C
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.
185C
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.
194C
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.
203C
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.
212C
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.
221C
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.
231C
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.
240C
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.
249C
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.
258C
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.
267C
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.
276C
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.
285C
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.
294C
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.
303C
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.
314C
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.
323C
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.
332C
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.
341C
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.
350C
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.
359C
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.
368C
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.
377C
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.
386C
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.
397C
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.
406C
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.
415C
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.
424C
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.
433C
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.
442C
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.
451C
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.
460C
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.
469C
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.
480C
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.
489C
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.
498C
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.
507C
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.
516C
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.
525C
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.
534C
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.
543C
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.
552C
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
600C** 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
619C** 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)
625C** 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)
631C --- 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)
636C** TRANSFORM INTO LAB.SYSTEM
637 CALL LOR(2+I,MXGKPV,2+I)
638 PV(7,2+I)=TOF
639 83 CONTINUE
640C** 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.
699C** SET INTERACTION MODE ACCORDING TO GHEISHA-CONVENTION
700C** 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)
712C
713 9999 CONTINUE
714 END