]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/fluka/dres.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Default compile option changed to -g (Alpha)
[u/mrichter/AliRoot.git] / GEANT321 / fluka / dres.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1995/10/24 10:19:54 cernlib
6* Geant
7*
8*
9#include "geant321/pilot.h"
10#if defined(CERNLIB_OLDNAME)
11*CMZ : 3.21/02 29/03/94 15.41.42 by S.Giani
12*-- Author :
13*=== dres =============================================================*
14* *
15 SUBROUTINE DRES ( M2, M3, T1, U, EREC, LOPPAR, JFISS )
16
17#include "geant321/dblprc.inc"
18#include "geant321/dimpar.inc"
19#include "geant321/iounit.inc"
20*
21*----------------------------------------------------------------------*
22* *
23* New version of DRES created by A.Ferrari & P.Sala, INFN - Milan *
24* *
25* Last change on 10-apr-93 by Alfredo Ferrari, INFN - Milan *
26* *
27* Dres93: Dres91 plus the RAL fission model taken from LAHET thanks *
28* to R.E.Prael *
29* Dres91: new version from A. Ferrari and P. Sala, INFN - Milan *
30* This routine has been adapted from the original one of the *
31* Evap-5 module (KFA - Julich). Main modifications concern *
32* with kinematics which is now fully relativistic and with *
33* the treatment of few nucleons nuclei, which are now frag- *
34* mented, even though in a very rough manner. Changes have *
35* been made also to other routines of the Evap-5 package *
36* *
37*----------------------------------------------------------------------*
38*
39*----------------------------------------------------------------------*
40* *
41* Input variables: *
42* M2 = Mass number of the residual nucleus *
43* M3 = Atomic number of the residual nucleus *
44* T1 = Excitation energy of the residual nucleus before evaporation*
45* U = Excitation energy of the residual nucleus after evaporation *
46* Erec = Recoil kinetic energy of the residual nucleus *
47* The recoil direction is given by Coslbr (i) *
48* *
49* Significant variables: *
50* JA = Present mass number of the residual nucleus *
51* JZ = Present atomic number of the residual nucleus *
52* Smom1 = Energy accumulators for the six types of evaporated *
53* particles *
54* *
55* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *
56* !!!! Please note that the following variables concerning !!!! *
57* !!!! with the present residual nucleus must be set before!!!! *
58* !!!! entering DRES91: Ammres, Ptres !!!! *
59* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *
60* *
61*----------------------------------------------------------------------*
62*
63C---------------------------------------------------------------------
64C SUBNAME = DRES --- EVAPORATION
65C EVAPORATION DATA SHOUD BE READ ON INPUT STAGE
66C---------------------------------------------------------------------
67C*****A LA EVAP III(TWA,8-68)
68#include "geant321/eva0.inc"
69#include "geant321/eva1.inc"
70#include "geant321/forcn.inc"
71#include "geant321/higfis.inc"
72#include "geant321/inpflg.inc"
73#include "geant321/labcos.inc"
74#include "geant321/nucdat.inc"
75#include "geant321/resnuc.inc"
76 DIMENSION ZMASS (6), Z2MASS(6), C(3), Q(0:6), FLKCOU(6), CCOUL(6),
77 & THRESH(6), SMALLA(6), R(6), S (6), SOS (6), STRUN(6),
78 & EYE1 (6), EYE0 (6), SMOM1 (6), BNMASS(6)
79 DIMENSION CORRRR(6)
80 REAL RNDM(2)
81 LOGICAL LOPPAR, PENBAR, LFIRST
82 SAVE ZMASS, Z2MASS, EMHN, EMNUM, UM, AMUMEV, AMEMEV, QBRBE8,
83 & BNMASS, IEVEVP, NBE8, NRNEEP, LFIRST
84 DATA IEVEVP / 0 /
85 DATA LFIRST / .TRUE. /
86*
87 IEVEVP = IEVEVP + 1
88C-------------------------------------- 1.ST CALL INIT
89 IF ( LFIRST ) THEN
90 LFIRST = .FALSE.
91 FKEY = ZERZER
92 NBE8 = 0
93 NRNEEP = 0
94 EXMASS(1) = 1.D+03 * ( AMNEUT - AMUAMU )
95 EXMASS(2) = ENERGY ( ONEONE, ONEONE )
96 EXMASS(3) = ENERGY ( TWOTWO, ONEONE )
97 EXMASS(4) = ENERGY ( THRTHR, ONEONE )
98 EXMASS(5) = ENERGY ( THRTHR, TWOTWO )
99 EXMASS(6) = ENERGY ( FOUFOU, TWOTWO )
100 ZMASS(1) = 1.D+03 * AMUAMU + EXMASS (1)
101 ZMASS(2) = 1.D+03 * AMUAMU + EXMASS (2)
102 ZMASS(3) = 2.D+03 * AMUAMU + EXMASS (3)
103 ZMASS(4) = 3.D+03 * AMUAMU + EXMASS (4)
104 ZMASS(5) = 3.D+03 * AMUAMU + EXMASS (5)
105 ZMASS(6) = 4.D+03 * AMUAMU + EXMASS (6)
106 BNMASS (1) = 0.D+00
107 BNMASS (2) = 0.D+00
108 BNMASS (3) = ZMASS (1) + ZMASS (2) - ZMASS (3)
109 BNMASS (4) = TWOTWO * ZMASS (1) + ZMASS (2) - ZMASS (4)
110 BNMASS (5) = ZMASS (1) + TWOTWO * ZMASS (2) - ZMASS (5)
111 BNMASS (6) = TWOTWO * ( ZMASS (1) + ZMASS (2) ) - ZMASS (6)
112 DO 1234 KK = 1,6
113 Z2MASS (KK) = ZMASS (KK) * ZMASS (KK)
1141234 CONTINUE
115 AMUMEV = 1.D+03 * AMUAMU
116 AMEMEV = 1.D+03 * AMELEC
117 QBRBE8 = ENERGY ( EIGEIG, FOUFOU ) - TWOTWO * EXMASS (6)
118 EMN = 1.D+03 * AMNEUT
119 EMH = ZMASS (2)
120 TMP16 = 16.D+00
121 UM = AMUMEV + ENERGY ( TMP16, EIGEIG ) / 16.D+00
122 EMHN = EMH - EMN
123 EMNUM = EMN - UM
124 END IF
125* | End of initialization:
126* +-------------------------------------------------------------------*
127C --------------------------------- START OF PROCESS
128* +-------------------------------------------------------------------*
129* | Initialize Npart and Smom if nothing has been already evaporated
130* | for this event
131 IF ( JFISS .LE. 0 ) THEN
132 DO 775 I=1,6
133 NPART(I) = 0
134 SMOM1(I) = ZERZER
135 775 CONTINUE
136 END IF
137* |
138* +-------------------------------------------------------------------*
139 JA = M2
140 JZ = M3
141 U = T1
142 RNMASS = 1.D+03 * AMMRES + U
143* P2res and Ptres are the squared momentum and the momentum of the
144* residual nucleus (now in relativistic kinematics), Umo the
145* invariant mass of the system!
146 UMO = RNMASS
147 UMO2 = UMO * UMO
148 ELBTOT = RNMASS + EREC
149 GAMCM = ELBTOT / RNMASS
150 ETACM = 1.D+03 * PTRES / RNMASS
151 HEVSUM = ZERZER
152 1000 CONTINUE
153 LOPPAR = .FALSE.
154* +-------------------------------------------------------------------*
155* | Check for starting data inconsistencies
156 IF (JA-JZ .LT. 0) THEN
157 WRITE(LUNOUT,6401)
158 WRITE(LUNERR,6401)
159 6401 FORMAT('1 Dres: cascade residual nucleus has mass no. less',
160 & ' than Z!!')
161 RETURN
162* |
163* +-------------------------------------------------------------------*
164* | Rough treatment for very few nucleon residual
165* | nuclei. The basic ideas are:
166* | a) as many as possible alpha particles are emitted
167* | b) particles are emitted one per time leaving a residual
168* | excitation energy proportional to number of nucleons
169* | left in the residual nucleus (so we deal only with
170* | two body kinematics)
171* | T A K E I N T O A C C O U N T T H A T T H I S
172* | T R E A T M E N T I S E X T R E M E L Y R O U G H
173* | T H E T A S K B E I N G O N L Y T O S U P P L Y
174* | S O M E T H I N G T O S H A R E E N E R G Y A N D
175* | M O M E N T U M A M O N G A F E W F R A G M E N T S
176 ELSE IF ( JA .LE. 6 .OR. JZ .LE. 2 ) THEN
177* | 1000 continue moved above according to FCA suggestion
178*1000 CONTINUE
179 JRESID = 0
180 IF ( JA .GT. 4 ) GO TO 2000
181* | +----------------------------------------------------------------*
182* | | First check we are not concerning with a couple of neutrons or
183* | | protons
184 IF ( JA .EQ. 2 .AND. JZ .NE. 1 ) THEN
185 JEMISS = 1 + JZ / 2
186 JRESID = JEMISS
187 RNMASS = ZMASS (JRESID)
188 U = 0.D+00
189 DELTU = UMO - 2.D+00 * ZMASS (JEMISS)
190* | | +-------------------------------------------------------------*
191* | | |
192 IF ( DELTU .LE. 0.D+00 ) THEN
193 IF ( DELTU .LT. - 2.D+00 * ANGLGB * UMO ) THEN
194 WRITE ( LUNERR, * )' *** Dres: insufficient Umo for',
195 & ' a nucleon couple', UMO,
196 & 2.D+00 * ZMASS (JEMISS)
197 END IF
198 UMO = ( UMO + DELTU ) * ( 1.D+00 + ANGLGB )
199 END IF
200* | | |
201* | | +-------------------------------------------------------------*
202 GO TO 2500
203 END IF
204* | |
205* | +----------------------------------------------------------------*
206* | +----------------------------------------------------------------*
207* | | Then check we are not concerning with one of the six
208* | | standard particles
209 DO 1700 J = 6, 1, -1
210* | | +-------------------------------------------------------------*
211* | | |
212 IF ( JZ .EQ. IZ (J) .AND. JA .EQ. IA (J) ) THEN
213 HEVSUM = SMOM1(3) + SMOM1(5) + SMOM1(6) + SMOM1(4)
214 GO TO ( 1100, 1100, 1600, 1500, 1400, 1300 ), J
215* | | | +----------------------------------------------------------*
216* | | | | Proton or neutron, nothing can be done
217 1100 CONTINUE
218 RETURN
219* | | | +----------------------------------------------------------*
220* | | | | Alpha:
221 1300 CONTINUE
222 DEUDEU = MAX ( ZERZER, U + TWOTWO * BNMASS (3)
223 & - BNMASS (6) )
224 PROTRI = MAX ( ZERZER, U + BNMASS (4) - BNMASS (6) )
225 UEU3HE = MAX ( ZERZER, U + BNMASS (5) - BNMASS (6) )
226 QNORM = DEUDEU + PROTRI + UEU3HE
227* | | | | If we cannot split then return
228 IF ( QNORM .LE. ZERZER ) RETURN
229 CALL GRNDM(RNDM,1)
230 V = RNDM (1)
231* | | | | +-------------------------------------------------------*
232* | | | | | Split or into two deuterons or a triton and a proton
233* | | | | | or a 3-He and a neutron: no account is made for
234* | | | | | Coulomb effects, probability is simply assumed
235* | | | | | proportional to reaction Qs
236 IF ( V .LT. DEUDEU / QNORM ) THEN
237* | | | | | Two deuterons selected
238 JEMISS = 3
239 JRESID = 3
240 RNMASS = ZMASS (3)
241 U = ZERZER
242 GO TO 2500
243* | | | | |
244* | | | | +-------------------------------------------------------*
245* | | | | | Split into a triton and a proton
246 ELSE IF ( V .LT. ( DEUDEU + PROTRI ) / QNORM ) THEN
247 JEMISS = 2
248 JRESID = 4
249 RNMASS = ZMASS (4)
250 U = ZERZER
251 GO TO 2500
252* | | | | |
253* | | | | +-------------------------------------------------------*
254* | | | | | Split into a 3-He and a neutron
255 ELSE
256 JEMISS = 1
257 JRESID = 5
258 RNMASS = ZMASS (5)
259 U = ZERZER
260 GO TO 2500
261 END IF
262* | | | | |
263* | | | | +-------------------------------------------------------*
264* | | | |
265* | | | +----------------------------------------------------------*
266* | | | | 3-He:
267 1400 CONTINUE
268 DEUPRO = MAX ( ZERZER, U + BNMASS (3) - BNMASS (5) )
269 PRPRNE = MAX ( ZERZER, U - BNMASS (5) )
270 QNORM = DEUPRO + PRPRNE
271* | | | | If we cannot split then return
272 IF ( QNORM .LE. ZERZER ) RETURN
273 CALL GRNDM(RNDM,1)
274 V = RNDM (1)
275* | | | | +-------------------------------------------------------*
276* | | | | | Split or into a deuteron and a proton
277* | | | | | or into two protons and one neutron: no account is
278* | | | | | made for Coulomb effects, probability is simply assumed
279* | | | | | prportional to reaction Qs
280 IF ( V .LT. DEUPRO / QNORM ) THEN
281* | | | | | A deuteron and a proton selected
282 JEMISS = 2
283 JRESID = 3
284 RNMASS = ZMASS (3)
285 U = ZERZER
286 GO TO 2500
287* | | | | |
288* | | | | +-------------------------------------------------------*
289* | | | | | Split into 2 protons and 1 neutron: part of the exci-
290* | | | | | tation energy is conserved to allow the further
291* | | | | | splitting of the deuteron
292 ELSE
293 JEMISS = 2
294 JRESID = 0
295 FACT = ONEONE
296* | | | | | +----------------------------------------------------*
297* | | | | | | Loop to compute the residual excitation energy
298 1450 CONTINUE
299 FACT = FACT * 0.6666666666666667D+00
300* | | | | | | Erncm, Eepcm are the total energies of the residual
301* | | | | | | nucleus and of the emitted particle in the CMS frame
302 U = FACT * PRPRNE + BNMASS (3)
303 RNMASS = ZMASS (3) + U
304 ERNCM = HLFHLF * ( UMO2 + RNMASS**2
305 & - Z2MASS (JEMISS) ) / UMO
306 EEPCM = UMO - ERNCM
307 IF ( EEPCM .LE. ZMASS (JEMISS) ) GO TO 1450
308* | | | | | |
309* | | | | | +----------------------------------------------------*
310 GO TO 2600
311 END IF
312* | | | | |
313* | | | | +-------------------------------------------------------*
314* | | | |
315* | | | +----------------------------------------------------------*
316* | | | | Triton:
317 1500 CONTINUE
318 DEUNEU = MAX ( ZERZER, U + BNMASS (3) - BNMASS (4) )
319 PRNENE = MAX ( ZERZER, U - BNMASS (4) )
320 QNORM = DEUNEU + PRNENE
321* | | | | If we cannot split then return
322 IF ( QNORM .LE. ZERZER ) RETURN
323 CALL GRNDM(RNDM,1)
324 V = RNDM (1)
325* | | | | +-------------------------------------------------------*
326* | | | | | Split or into a deuteron and a neutron
327* | | | | | or into two protons and one neutron: no account is
328* | | | | | made for Coulomb effects, probability is simply assumed
329* | | | | | proportional to reaction Qs
330 IF ( V .LT. DEUNEU / QNORM ) THEN
331* | | | | | A deuteron and a proton selected
332 JEMISS = 1
333 JRESID = 3
334 RNMASS = ZMASS (3)
335 U = ZERZER
336 GO TO 2500
337* | | | | |
338* | | | | +-------------------------------------------------------*
339* | | | | | Split into 1 proton and 2 neutrons: part of the exci-
340* | | | | | tation energy is conserved to allow the further
341* | | | | | splitting of the deuteron
342 ELSE
343 JEMISS = 1
344 JRESID = 0
345 FACT = ONEONE
346* | | | | | +----------------------------------------------------*
347* | | | | | | Loop to compute the residual excitation energy
348 1550 CONTINUE
349 FACT = FACT * 0.6666666666666667D+00
350* | | | | | | Erncm, Eepcm are the total energies of the residual
351* | | | | | | nucleus and of the emitted particle in the CMS frame
352 U = FACT * PRNENE + BNMASS (3)
353 RNMASS = ZMASS (3) + U
354 ERNCM = HLFHLF * ( UMO2 + RNMASS**2
355 & - Z2MASS (JEMISS) ) / UMO
356 EEPCM = UMO - ERNCM
357 IF ( EEPCM .LE. ZMASS (JEMISS) ) GO TO 1550
358* | | | | | |
359* | | | | | +----------------------------------------------------*
360 GO TO 2600
361 END IF
362* | | | | |
363* | | | | +-------------------------------------------------------*
364* | | | |
365* | | | +----------------------------------------------------------*
366* | | | | Deuteron:
367 1600 CONTINUE
368* | | | | +-------------------------------------------------------*
369* | | | | | Split into a proton and a neutron if it is possible
370 IF ( U .GT. BNMASS (3) ) THEN
371 JEMISS = 1
372 JRESID = 2
373 RNMASS = ZMASS (2)
374 U = ZERZER
375 GO TO 2500
376* | | | | |
377* | | | | +-------------------------------------------------------*
378* | | | | | Energy too low to split the deuteron, return
379 ELSE
380 WRITE (LUNERR,*)' **Dres: energy too low to split',
381 & ' a deuteron! M2,M3',M2,M3
382 RETURN
383 END IF
384* | | | | |
385* | | | | +-------------------------------------------------------*
386 END IF
387* | | |
388* | | +-------------------------------------------------------------*
389 1700 CONTINUE
390* | |
391* | +----------------------------------------------------------------*
392 2000 CONTINUE
393 A = JA
394 Z = JZ
395 Q (0) = ZERZER
396 ENERG0 = ENERGY (A,Z)
397* | +----------------------------------------------------------------*
398* | | Note that Q(i) are not the reaction Qs but the remaining
399* | | energy after the reaction
400 DO 2100 K = 1, 6
401 JJA = JA - IA (K)
402 JJZ = JZ - IZ (K)
403 JJN = JJA - JJZ
404 IF ( JJN .LT. 0 .OR. JJZ .LT. 0 ) THEN
405 Q (K) = Q (K-1)
406 GO TO 2100
407 END IF
408 DDJJA = JJA
409 DDJJZ = JJZ
410 Q (K) = MAX ( U + ENERG0 - ENERGY ( DDJJA, DDJJZ )
411 & - EXMASS (K), ZERZER ) + Q (K-1)
412 2100 CONTINUE
413* | |
414* | +----------------------------------------------------------------*
415* | +----------------------------------------------------------------*
416* | | If no emission channel is open then return
417 IF ( Q (6) .LE. ZERZER ) THEN
418 HEVSUM = SMOM1(3) + SMOM1(5) + SMOM1(6) + SMOM1(4)
419 RETURN
420 END IF
421* | |
422* | +----------------------------------------------------------------*
423 CALL GRNDM(RNDM,1)
424 V = RNDM (1)
425 FACT = ONEONE
426* | +----------------------------------------------------------------*
427* | |
428 DO 2200 J = 1, 6
429* | | +-------------------------------------------------------------*
430* | | |
431 IF ( V .LT. Q (J) / Q (6) ) THEN
432 JEMISS = J
433 JJA = JA - IA (JEMISS)
434 JJZ = JZ - IZ (JEMISS)
435* | | | +----------------------------------------------------------*
436* | | | |
437 DO 2150 JJ = 1, 6
438* | | | | +-------------------------------------------------------*
439* | | | | |
440 IF ( JJA .EQ. IA (JJ) .AND. JJZ .EQ. IZ (JJ) ) THEN
441 JRESID = JJ
442 RNMASS = ZMASS (JRESID)
443 ERNCM = HLFHLF * ( UMO2 + Z2MASS (JRESID)
444 & - Z2MASS (JEMISS) ) / UMO
445 EEPCM = UMO - ERNCM
446 U = ZERZER
447 GO TO 2600
448 END IF
449* | | | | |
450* | | | | +-------------------------------------------------------*
451 2150 CONTINUE
452* | | | |
453* | | | +----------------------------------------------------------*
454 AJJA = JJA
455 ZJJZ = JJZ
456 RNMAS0 = AJJA * AMUMEV + ENERGY ( AJJA, ZJJZ )
457 GO TO 2300
458 END IF
459* | | |
460* | | +-------------------------------------------------------------*
461 2200 CONTINUE
462* | |
463* | +----------------------------------------------------------------*
464 WRITE ( LUNOUT,* )' **** error in Dres, few nucleon treatment',
465 & ' ****'
466 WRITE ( LUNERR,* )' **** error in Dres, few nucleon treatment',
467 & ' ****'
468 RETURN
469* | +----------------------------------------------------------------*
470* | | Loop to compute the residual excitation energy
471 2300 CONTINUE
472 FACT = FACT * AJJA / A
473 U = FACT * ( Q (JEMISS) - Q (JEMISS-1) )
474* | | Erncm, Eepcm are the total energies of the residual
475* | | nucleus and of the emitted particle in the CMS frame
476 RNMASS = RNMAS0 + U
477 ERNCM = HLFHLF * ( UMO2 + RNMASS**2
478 & - Z2MASS (JEMISS) ) / UMO
479 EEPCM = UMO - ERNCM
480 IF ( EEPCM .LE. ZMASS (JEMISS) ) THEN
481 IF ( Q (JEMISS) - Q (JEMISS-1) .GE. 1.D-06 ) GO TO 2300
482* | +--<--<--<--<--<--< Loop back
483* | | Actually there is no excitation energy available!
484 U = ANGLGB
485 RNMASS = ONEPLS * RNMAS0
486 ERNCM = ONEPLS * RNMASS
487 EEPCM = ONEPLS * ZMASS (JEMISS)
488 END IF
489* | |
490* | +----------------------------------------------------------------*
491 GO TO 2600
492* | From here standard two bodies kinematics with Jemiss, Rnmass
493* | (new excitation energy is U)
494 2500 CONTINUE
495* | Erncm, Eepcm are the total energies of the residual
496* | nucleus and of the emitted particle in the CMS frame
497 ERNCM = HLFHLF * ( UMO2 + RNMASS**2 - Z2MASS (JEMISS) ) / UMO
498 EEPCM = UMO - ERNCM
499 2600 CONTINUE
500* | C(i) are the direction cosines of the emitted particle
501* | (Jemiss) in the CMS frame, of course - C(i)
502* | are the ones of the residual nucleus (Jresid if one of the
503* | standard partcles, say the proton)
504 CALL RACO (C(1),C(2),C(3))
505 PCMS = SQRT ( EEPCM**2 - Z2MASS (JEMISS) )
506* | Now we perform the Lorentz transformation back to the original
507* | frame (lab frame)
508* | First the emitted particle:
509 ETAX = ETACM * COSLBR (1)
510 ETAY = ETACM * COSLBR (2)
511 ETAZ = ETACM * COSLBR (3)
512 PCMSX = PCMS * C (1)
513 PCMSY = PCMS * C (2)
514 PCMSZ = PCMS * C (3)
515 ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
516 EPS = GAMCM * EEPCM + ETAPCM - ZMASS (JEMISS)
517 PHELP = ETAPCM / ( GAMCM + ONEONE ) + EEPCM
518 PLBPX = PCMSX + ETAX * PHELP
519 PLBPY = PCMSY + ETAY * PHELP
520 PLBPZ = PCMSZ + ETAZ * PHELP
521 PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
522 COSLBP (1) = PLBPX / PHELP
523 COSLBP (2) = PLBPY / PHELP
524 COSLBP (3) = PLBPZ / PHELP
525* | Then the residual nucleus ( for it c (i) --> - c (i) ):
526 EREC = GAMCM * ERNCM - ETAPCM - RNMASS
527 EKRES = 1.D-03 * EREC
528 PHELP = - ETAPCM / ( GAMCM + ONEONE ) + ERNCM
529 PXRES = 1.D-03 * ( - PCMSX + ETAX * PHELP )
530 PYRES = 1.D-03 * ( - PCMSY + ETAY * PHELP )
531 PZRES = 1.D-03 * ( - PCMSZ + ETAZ * PHELP )
532 P2RES = PXRES * PXRES + PYRES * PYRES + PZRES * PZRES
533 PTRES = SQRT (P2RES)
534 COSLBR (1) = PXRES / PTRES
535 COSLBR (2) = PYRES / PTRES
536 COSLBR (3) = PZRES / PTRES
537* | Score the emitted particle
538 NPART (JEMISS) = NPART (JEMISS) + 1
539 SMOM1 (JEMISS) = SMOM1 (JEMISS) + EPS
540 ITEMP=NPART(JEMISS)
541 EPART(ITEMP,JEMISS)=EPS
542 COSEVP(1,ITEMP,JEMISS)=COSLBP(1)
543 COSEVP(2,ITEMP,JEMISS)=COSLBP(2)
544 COSEVP(3,ITEMP,JEMISS)=COSLBP(3)
545* | +----------------------------------------------------------------*
546* | | Check if the residual nucleus is one of the emitted particles
547 IF ( JRESID .GT. 0 ) THEN
548 J = JRESID
549 GO TO 6102
550 END IF
551* | |
552* | +----------------------------------------------------------------*
553 JA = JA - IA (JEMISS)
554 JZ = JZ - IZ (JEMISS)
555* Umo is the invariant mass of the system!!
556 UMO = RNMASS
557 UMO2 = UMO * UMO
558 ELBTOT = RNMASS + EREC
559 GAMCM = ELBTOT / RNMASS
560 ETACM = 1.D+03 * PTRES / RNMASS
561 GO TO 1000
562 END IF
563* |
564* +-------------------------------------------------------------------*
565* Come to 23 at the beginning and after the end of a "normal"
566* evaporation cycle
567 23 CONTINUE
568 A = JA
569 Z = JZ
570 IF(JA-8)8486,8488,8486
571 8488 CONTINUE
572 IF(JZ-4)8486,1224,8486
573 8486 CONTINUE
574 DO 2 K=1,6
575 IF((A-FLA(K)).LE.(Z-FLZ(K))) THEN
576 Q(K)=99999.D0
577 GO TO 2
578 END IF
579 IF(A-TWOTWO*FLA(K))2,727,727
580 727 CONTINUE
581 IF(Z-TWOTWO*FLZ(K))2,728,728
582 728 CONTINUE
583 Q(K) = QNRG(A-FLA(K),Z-FLZ(K),A,Z) + EXMASS(K)
584C 728 Q(K)=ENERGY(A-FLA(K),Z-FLZ(K))-ENERGY(A,Z)+EXMASS(K)
585 2 CONTINUE
586 FLKCOU(1)=ZERZER
587 FLKCOU(2)=DOST(1,Z-FLZ(2))
588 FLKCOU(3)=FLKCOU(2)+.06D+00
589 FLKCOU(4)=FLKCOU(2)+.12D+00
590 FLKCOU(6)=DOST(2,Z-FLZ(6))
591 FLKCOU(5)=FLKCOU(6)-.06D+00
592 CCOUL(1)=ONEONE
593 CCOU2=DOST(3,Z-FLZ(2))
594 CCOUL(2)=CCOU2+ONEONE
595 CCOUL(3)=CCOU2*1.5D0+THRTHR
596 CCOUL(4)=CCOU2+THRTHR
597 CCOUL(6)=DOST(4,Z-FLZ(6))*TWOTWO+TWOTWO
598 CCOUL(5)=TWOTWO*CCOUL(6)-ONEONE
599 SIGMA=ZERZER
600* Initialize the flag which checks for open particle decay with
601* zero excitation and pairing --> for particle unstable residual
602* nuclei
603 LOPPAR = .FALSE.
604 DO 33 J=1,6
605 IF(A-TWOTWO*FLA(J))5,725,725
606 725 CONTINUE
607 IF(Z-TWOTWO*FLZ(J))5,726,726
608 726 CONTINUE
609 MM=JA-IA(J)
610 ZZ=Z-FLZ(J)
611 AA=A-FLA(J)
612 IF(AA.LE.ZZ)GO TO 5
613* Energy threshold for the emission of the jth-particle
614 THRESH(J)=Q(J)+.88235D+00*FLKCOU(J)*FLZ(J)*
615 & ZZ/(RMASS(MM)+RHO(J))
616 LOPPAR = LOPPAR .OR. THRESH (J) .GT. ZERZER
617 IAA = NINT(AA)
618 IZZ = NINT(ZZ)
619 NN = IAA - IZZ
620* The residual nucleus excitation energy ranges from 0 up
621* to U - Q (J)
622 ILVMOD = IB0
623 UMXRES = U - THRESH (J)
624* This is the a lower case of the level density
625 SMALLA (J) = GETA ( UMXRES, IZZ, NN, ILVMOD, ISDUM, ASMMAX,
626 & ASMMIN )
627 CALL EEXLVL (IAA,IZZ,EEX1ST,EEX2ND,CORR)
628 EEX1ST = 1.D+03 * EEX1ST
629 EEX2ND = 1.D+03 * EEX2ND
630 CORR = 1.D+03 * MAX ( CORR, ZERZER )
631 IF ( NN .EQ. 4 .AND. IZZ .EQ. 4 ) THEN
632 IF ( U - THRESH (J) - 6.1D+00 .GT. ZERZER ) THEN
633 CORR = SIXSIX
634 ELSE
635 TMPVAR = U-THRESH(J)-0.1D+00
636 CORR = MAX ( ZERZER, TMPVAR )
637 END IF
638 END IF
639 IF (NINT(FKEY).EQ.1) CORR=ZERZER
640 CORRRR(J)=CORR
641* Standard calculation:
642 ARG=U-THRESH(J)-CORR
643 IF(ARG)5,4,4
644 5 CONTINUE
645 R(J)=ZERZER
646 S(J)=ZERZER
647 SOS(J)=ZERZER
648 GO TO 33
649 4 CONTINUE
650 S(J)=SQRT (SMALLA(J)*ARG)*TWOTWO
651 SOS(J)=TENTEN*S(J)
652 33 CONTINUE
653 N1=1
654 SES = MAX (S(1),S(2),S(3),S(4),S(5),S(6))
655 DO 1131 J=1,6
656 JS = SOS(J) + ONEONE
657 FJS = JS
658 STRUN(J) = FJS - ONEONE
659 IF ( S(J) .GT. ZERZER ) THEN
660 MM = JA-IA(J)
661 SAS = EXP (S(J)-SES)
662 SUS = EXP (-S(J))
663 EYE1(J) = ( TWOTWO * S(J)**2 -SIXSIX * S(J)
664 & + SIXSIX + SUS * ( S(J)**2 - SIXSIX ) )
665 & / ( EIGEIG * SMALLA(J)**2 )
666 IF ( J .EQ. 1 ) THEN
667 EYE0(J) = ( S(J) - ONEONE + SUS ) / ( TWOTWO*SMALLA(J) )
668* Standard calculation
669 R (J) = RMASS(MM)**2 * ALPH(MM) * ( EYE1(J) + BET(MM)
670 & * EYE0(J) ) * SAS
671 ELSE
672 R (J) = CCOUL(J) * RMASS(MM)**2 * EYE1(J) * SAS
673 END IF
674 R (J) = MAX ( ZERZER, R (J) )
675 SIGMA = SIGMA + R (J)
676 END IF
677 1131 CONTINUE
678 NCOUNT = 0
679 6202 CONTINUE
680 IF(SIGMA)9,9,10
681 9 CONTINUE
682 DO 6100 J = 1,6
683 IF(JA-IA(J))6100,6101,6100
684 6101 CONTINUE
685 IF(JZ-IZ(J))6100,6102,6100
686 6100 CONTINUE
687 GO TO 6099
688 6102 CONTINUE
689 IF ( U .GT. ANGLGB ) GO TO 1000
690 JEMISS = J
691C*****STORE,RESIDUAL NUC IS OF EMITTED PARTICLE TYPE
692* If we are here this means that the residual nucleus is equal to
693* one of the six emitted particle (the j-th one). So give to it
694* all the energy, score it and return with 0 recoil and excitation
695* energy for the residual nucleus
696 EPS = EREC
697 NPART(JEMISS) = NPART(JEMISS)+1
698 ITEMP=NPART(JEMISS)
699 NRNEEP = NRNEEP + 1
700 SMOM1(JEMISS) = SMOM1(JEMISS) + EPS
701 ITEMP=NPART(JEMISS)
702 EPART(ITEMP,JEMISS)=EPS
703 COSEVP(1,ITEMP,JEMISS) = COSLBR(1)
704 COSEVP(2,ITEMP,JEMISS) = COSLBR(2)
705 COSEVP(3,ITEMP,JEMISS) = COSLBR(3)
706 GO TO 8002
707*
708*-->-->-->-->--> go to return
709 6099 CONTINUE
710 IF(JA-8)72,51,72
711 51 CONTINUE
712 IF(JZ-4)72,1224,72
713* Come here for a "normal" step
714 10 CONTINUE
715 LOPPAR = .FALSE.
716 CALL GRNDM(RNDM,1)
717 URAN=RNDM(1)*SIGMA
718 SUM = ZERZER
719 DO 16 J=1,6
720 K = J
721 SUM = R(J)+SUM
722 IF ( SUM - URAN .GT. 0.D+00 ) GO TO 17
723 16 CONTINUE
724 17 CONTINUE
725 JEMISS=K
726 NPART(JEMISS) = NPART (JEMISS) + 1
727 JS = SOS (JEMISS) + ONEONE
728 IF(JS-1000)4344,4345,4345
729 4345 CONTINUE
730 RATIO2=(S(JEMISS)**3-6.D0*S(JEMISS)**2+15.D0*
731 &S(JEMISS)-15.D0)/((2.D0*S(JEMISS)**2-6.D0*S(JEMISS)+6.D0)*SMALLA
732 &(JEMISS))
733 GO TO 4347
734 4344 CONTINUE
735 RATIO2=(P2(JS)+(P2(JS+1)-P2(JS))*
736 &(SOS(JEMISS)-STRUN(JEMISS)))/SMALLA(JEMISS)
737 4347 CONTINUE
738 EPSAV=RATIO2*TWOTWO
739* +-------------------------------------------------------------------*
740* | Neutron channel selected:
741 IF (JEMISS .EQ. 1) THEN
742 MM=JA-IA(J)
743 EPSAV=(EPSAV+BET(MM))/(ONEONE+BET(MM)*EYE0(JEMISS)
744 & /EYE1(JEMISS))
745* | +----------------------------------------------------------------*
746* | | Compute the fission width relative to the neutron one:
747* | | this part is taken from subroutine EMIT of LAHET
748 IF ( IFISS .GT. 0 .AND. JZ .GE. 71 .AND. .NOT. FISINH ) THEN
749* | | +-------------------------------------------------------------*
750* | | | Compute the correction factor for the fission width:
751 IF ( JZ .GT. 88 ) THEN
752 AGOES = ONEONE
753* | | |
754* | | +-------------------------------------------------------------*
755* | | |
756 ELSE
757 AGOES = MAX ( ONEONE, ( U-SEVSEV ) / ( EPSAV+SEVSEV ) )
758 END IF
759* | | |
760* | | +-------------------------------------------------------------*
761* | | Finally this is the relative fission width:
762* | | This is : Probfs = 1 / ( 1 + G_n / G_f )
763 PROBFS = FPROB ( Z, A, U ) / AGOES
764* | | +-------------------------------------------------------------*
765* | | | Check if it will be fission:
766 CALL GRNDM(RNDM,1)
767 IF ( RNDM (1) .LT. PROBFS ) THEN
768 FISINH = .TRUE.
769 KFISS = 1
770* | | | Undo the counting of the neutron evaporation
771 NPART (JEMISS) = NPART (JEMISS) -1
772 GO TO 9260
773 END IF
774* | | |
775* | | +-------------------------------------------------------------*
776 END IF
777* | |
778* | +----------------------------------------------------------------*
779 END IF
780* |
781* +-------------------------------------------------------------------*
782 20 CONTINUE
783 CALL GRNDM(RNDM,2)
784 E1=-HLFHLF*LOG(RNDM(1))
785 E2=-HLFHLF*LOG(RNDM(2))
786* Eps should be the total kinetic energy in the CMS frame
787* Standard calculation:
788 EPS=(E1+E2)*EPSAV+THRESH(JEMISS)-Q(JEMISS)
789 AR = A - IA(JEMISS)
790 ZR = Z - IZ(JEMISS)
791* The CMS energy is updated
792 IMASS = NINT (AR)
793 IF ( IMASS .EQ. 8 .AND. NINT (ZR) .EQ. 4 ) THEN
794 UNEW = U - EPS - Q(JEMISS)
795 UMAX = U - THRESH(JEMISS)
796 IF ( UNEW .GT. 6.D+00 ) THEN
797 UMIN = 6.D+00
798 ELSE IF ( UNEW .GT. 4.47D+00 .AND. UMAX .GT. 6.D+00 ) THEN
799 UMIN = 4.47D+00
800 UNEW = 6.D+00
801 ELSE IF ( UNEW .GT. 1.47D+00 .AND. UMAX .GT. 2.94D+00 ) THEN
802 UMIN = 1.47D+00
803 UNEW = 2.94D+00
804 ELSE
805 UMIN = -0.1D+00
806 UNEW = ANGLGB * 0.1D+00
807 END IF
808 ELSE IF ( IMASS .LE. 4 ) THEN
809 IPRO = NINT ( ZR )
810 INEU = IMASS - IPRO
811 IF ( IMASS .EQ. 1 ) THEN
812* Be sure that residual neutrons or protons are not left excited
813 UMIN = 0.D+00
814 UNEW = 0.D+00
815 EPS = U - Q(JEMISS)
816 ELSE IF ( IPRO .EQ. 0 .OR. INEU .EQ. 0 ) THEN
817* Ipro protons or ineu neutrons arrived here!
818 UMIN = CORRRR(JEMISS)
819 UNEW = U - EPS - Q(JEMISS)
820 ELSE IF ( IMASS .LE. 2 ) THEN
821* Be sure that residual deuterons are not left excited!
822 UMIN = 0.D+00
823 UNEW = 0.D+00
824 EPS = U - Q(JEMISS)
825 ELSE IF ( ABS ( INEU - IPRO ) .LE. 1 ) THEN
826* For the moment also residual 3-H, 3-He and 4-He are not left
827* excited !
828 UMIN = 0.D+00
829 UNEW = 0.D+00
830 EPS = U - Q(JEMISS)
831 ELSE
832 UMIN = CORRRR(JEMISS)
833 UNEW = U - EPS - Q(JEMISS)
834 END IF
835 ELSE
836 UMIN = CORRRR(JEMISS)
837 UNEW = U - EPS - Q(JEMISS)
838 END IF
839* Standard calculation
840 IF(UNEW-UMIN)6200,6220,6220
841 6220 CONTINUE
842*or RNMASS = AR * AMUMEV + ENERGY (AR,ZR)
843 RNMASS = AR * AMUMEV + ENERGY (AR,ZR) + UNEW
844 UMIN2 = ( RNMASS + ZMASS (JEMISS) )**2
845 IF ( UMIN2 .GE. UMO2 ) THEN
846 GO TO 6200
847 END IF
848 U = UNEW
849* C(i) are the direction cosines of the evaporated particle in the CMS
850* frame, of course - C(i) are the ones of the residual nucleus
851 CALL RACO(C(1),C(2),C(3))
852* Erncm, Eepcm are the total energies of the residual nucleus and
853* of the evaporated particle in the CMS frame
854 ERNCM = 0.5D+00 * ( UMO2 + RNMASS**2 - Z2MASS (JEMISS) ) / UMO
855 EEPCM = UMO - ERNCM
856 PCMS = SQRT ( EEPCM**2 - Z2MASS (JEMISS) )
857* Now we perform the Lorentz transformation back to the original
858* frame (lab frame)
859* First the evaporated particle:
860 ETAX = ETACM * COSLBR (1)
861 ETAY = ETACM * COSLBR (2)
862 ETAZ = ETACM * COSLBR (3)
863 PCMSX = PCMS * C (1)
864 PCMSY = PCMS * C (2)
865 PCMSZ = PCMS * C (3)
866 ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
867 EPS = GAMCM * EEPCM + ETAPCM - ZMASS (JEMISS)
868 PHELP = ETAPCM / (GAMCM + 1.D0) + EEPCM
869 PLBPX = PCMSX + ETAX * PHELP
870 PLBPY = PCMSY + ETAY * PHELP
871 PLBPZ = PCMSZ + ETAZ * PHELP
872 PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
873 COSLBP (1) = PLBPX / PHELP
874 COSLBP (2) = PLBPY / PHELP
875 COSLBP (3) = PLBPZ / PHELP
876* Then the residual nucleus ( for it c (i) --> - c (i) ):
877 EREC = GAMCM * ERNCM - ETAPCM - RNMASS
878 EKRES = 1.D-03 * EREC
879 PHELP = - ETAPCM / (GAMCM + 1.D0) + ERNCM
880 PXRES = 1.D-03 * ( - PCMSX + ETAX * PHELP )
881 PYRES = 1.D-03 * ( - PCMSY + ETAY * PHELP )
882 PZRES = 1.D-03 * ( - PCMSZ + ETAZ * PHELP )
883 P2RES = PXRES * PXRES + PYRES * PYRES + PZRES * PZRES
884 PTRES = SQRT (P2RES)
885 COSLBR (1) = PXRES / PTRES
886 COSLBR (2) = PYRES / PTRES
887 COSLBR (3) = PZRES / PTRES
888* Check energy and momentum conservation !!
889 IF (EREC .LE. 0.D+00) THEN
890 PTRES = 0.D+00
891 EREC = 0.D+00
892 END IF
893* Umo is the invariant mass of the system!!
894 UMO = RNMASS
895 UMO2 = UMO * UMO
896 ELBTOT = RNMASS + EREC
897 GAMCM = ELBTOT / RNMASS
898 ETACM = 1.D+03 * PTRES / RNMASS
899 GO TO 76
900 6200 CONTINUE
901 NCOUNT = NCOUNT + 1
902 IF ( NCOUNT .LE. 10 ) GO TO 20
903 SIGMA = SIGMA - R(JEMISS)
904* if we are here we have sampled for > 10 times a negative energy Unew
905 NPART(JEMISS)=NPART(JEMISS)-1
906 R(JEMISS) = 0.D0
907 NCOUNT = 0
908 GO TO 6202
909 76 CONTINUE
910 JAT=JA-IA(JEMISS)
911 JZT=JZ-IZ(JEMISS)
912 IF(JAT-JZT)9,9,77
913 77 CONTINUE
914 JA=JAT
915 JZ=JZT
916C*****STORE,END OF NORMAL CYCLE
917 SMOM1(JEMISS)=SMOM1(JEMISS)+EPS
918 ITEMP=NPART(JEMISS)
919 EPART(ITEMP,JEMISS)=EPS
920 COSEVP(1,ITEMP,JEMISS)=COSLBP(1)
921 COSEVP(2,ITEMP,JEMISS)=COSLBP(2)
922 COSEVP(3,ITEMP,JEMISS)=COSLBP(3)
923* The following card switch to the rough splitting treatment
924 IF (JA .LE. 2) GO TO 1000
925 IF(JA-8)23,1223,23
926 1223 CONTINUE
927 IF(JZ-4)23,1224,23
928* If we are here the residual nucleus is a 8-Be one, break it into
929* two alphas with all the available energy (U plus the Q of the breakup)
930* , score them and return with 0 recoil and excitation energy
931 1224 CONTINUE
932 LOPPAR = .FALSE.
933 IF(U)1228,1229,1229
934 1228 CONTINUE
935 EPS=0.D0
936 GO TO 1230
937 1229 CONTINUE
938 1230 CONTINUE
939 NBE8=NBE8+1
940* C(i) are the direction cosines of the first alpha in the CMS
941* frame, of course - C(i) are the ones of the other
942 CALL RACO(C(1),C(2),C(3))
943* Ecms is the total energy of the alphas in the CMS frame
944 ECMS = 0.5D+00 * UMO
945 PCMS = SQRT ( ECMS**2 - Z2MASS (6) )
946* Now we perform the Lorentz transformation back to the original
947* frame (lab frame)
948* First alpha:
949 ETAX = ETACM * COSLBR (1)
950 ETAY = ETACM * COSLBR (2)
951 ETAZ = ETACM * COSLBR (3)
952 PCMSX = PCMS * C (1)
953 PCMSY = PCMS * C (2)
954 PCMSZ = PCMS * C (3)
955 ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
956 EPS = GAMCM * ECMS + ETAPCM - ZMASS (6)
957 PHELP = ETAPCM / (GAMCM + 1.D0) + ECMS
958 PLBPX = PCMSX + ETAX * PHELP
959 PLBPY = PCMSY + ETAY * PHELP
960 PLBPZ = PCMSZ + ETAZ * PHELP
961 PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
962* Store the first alpha!!
963 SMOM1(6) = SMOM1(6) + EPS
964 NPART(6) = NPART(6) + 1
965 ITEMP = NPART(6)
966 EPART (ITEMP,6) = EPS
967 COSEVP (1,ITEMP,6) = PLBPX / PHELP
968 COSEVP (2,ITEMP,6) = PLBPY / PHELP
969 COSEVP (3,ITEMP,6) = PLBPZ / PHELP
970* Then the second alpha ( for it c (i) --> - c (i) ):
971 EPS = GAMCM * ECMS - ETAPCM - ZMASS (6)
972 PHELP = - ETAPCM / (GAMCM + 1.D0) + ECMS
973 PLBPX = - PCMSX + ETAX * PHELP
974 PLBPY = - PCMSY + ETAY * PHELP
975 PLBPZ = - PCMSZ + ETAZ * PHELP
976 PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
977* Store the second alpha !!
978 SMOM1(6) = SMOM1(6) + EPS
979 NPART(6) = NPART(6) + 1
980 ITEMP = NPART(6)
981 EPART (ITEMP,6) = EPS
982 COSEVP (1,ITEMP,6) = PLBPX / PHELP
983 COSEVP (2,ITEMP,6) = PLBPY / PHELP
984 COSEVP (3,ITEMP,6) = PLBPZ / PHELP
985 8002 CONTINUE
986 LOPPAR = .FALSE.
987 EREC = ZERZER
988 U = ZERZER
989 EKRES = ZERZER
990 PTRES = ZERZER
991 72 CONTINUE
992 HEVSUM=SMOM1(3)+SMOM1(5)+SMOM1(6)+SMOM1(4)
993 RETURN
994*.......................................................................
995*///// RAL FISSION ROUTINE /////
996 9260 CONTINUE
997* +-------------------------------------------------------------------*
998* | Record the direction cosines of the fissioning nucleus
999 DO 9270 I=1,3
1000 COSLF0 (I) = COSLBR (I)
1001 9270 CONTINUE
1002* |
1003* +-------------------------------------------------------------------*
1004 CALL FISFRA ( JA, JZ, U, EREC, UMO, GAMCM, ETACM )
1005* +-------------------------------------------------------------------*
1006* | Check for fission failures!!
1007 IF ( .NOT. FISINH ) THEN
1008 PENBAR = .FALSE.
1009 GO TO 23
1010 END IF
1011* |
1012* +-------------------------------------------------------------------*
1013* Do not pick up the fission fragments, rather go back to Evevap
1014 HEVSUM=SMOM1(3)+SMOM1(5)+SMOM1(6)+SMOM1(4)
1015 RETURN
1016 END
1017#endif