]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/fluka/fkdres.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 / fkdres.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1995/10/24 10:20:05 cernlib
6* Geant
7*
8*
9#include "geant321/pilot.h"
10*CMZ : 3.21/02 29/03/94 15.41.45 by S.Giani
11*-- Author :
12*=== dres =============================================================*
13* *
14 SUBROUTINE FKDRES ( M2, M3, T1, U, EREC, LOPPAR, JFISS )
15
16#include "geant321/dblprc.inc"
17#include "geant321/dimpar.inc"
18#include "geant321/iounit.inc"
19*
20*----------------------------------------------------------------------*
21* *
22* New version of DRES created by A.Ferrari & P.Sala, INFN - Milan *
23* *
24* Last change on 10-apr-93 by Alfredo Ferrari, INFN - Milan *
25* *
26* Dres93: Dres91 plus the RAL fission model taken from LAHET thanks *
27* to R.E.Prael *
28* Dres91: new version from A. Ferrari and P. Sala, INFN - Milan *
29* This routine has been adapted from the original one of the *
30* Evap-5 module (KFA - Julich). Main modifications concern *
31* with kinematics which is now fully relativistic and with *
32* the treatment of few nucleons nuclei, which are now frag- *
33* mented, even though in a very rough manner. Changes have *
34* been made also to other routines of the Evap-5 package *
35* *
36*----------------------------------------------------------------------*
37*
38*----------------------------------------------------------------------*
39* *
40* Input variables: *
41* M2 = Mass number of the residual nucleus *
42* M3 = Atomic number of the residual nucleus *
43* T1 = Excitation energy of the residual nucleus before evaporation*
44* U = Excitation energy of the residual nucleus after evaporation *
45* Erec = Recoil kinetic energy of the residual nucleus *
46* The recoil direction is given by Coslbr (i) *
47* *
48* Significant variables: *
49* JA = Present mass number of the residual nucleus *
50* JZ = Present atomic number of the residual nucleus *
51* Smom1 = Energy accumulators for the six types of evaporated *
52* particles *
53* *
54* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *
55* !!!! Please note that the following variables concerning !!!! *
56* !!!! with the present residual nucleus must be set before!!!! *
57* !!!! entering DRES91: Ammres, Ptres !!!! *
58* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *
59* *
60*----------------------------------------------------------------------*
61*
62C---------------------------------------------------------------------
63C SUBNAME = DRES --- EVAPORATION
64C EVAPORATION DATA SHOUD BE READ ON INPUT STAGE
65C---------------------------------------------------------------------
66C*****A LA EVAP III(TWA,8-68)
67#include "geant321/eva0.inc"
68#include "geant321/eva1.inc"
69#include "geant321/forcn.inc"
70#include "geant321/higfis.inc"
71#include "geant321/inpflg.inc"
72#include "geant321/labcos.inc"
73#include "geant321/nucdat.inc"
74#include "geant321/resnuc.inc"
75 DIMENSION ZMASS (6), Z2MASS(6), C(3), Q(0:6), FLKCOU(6), CCOUL(6),
76 & THRESH(6), SMALLA(6), R(6), S (6), SOS (6), STRUN(6),
77 & EYE1 (6), EYE0 (6), SMOM1 (6), BNMASS(6)
78 DIMENSION CORRRR(6)
79 REAL RNDM(2)
80 LOGICAL LOPPAR, PENBAR, LFIRST
81 PARAMETER (EXPMIN=-88,EXPMAX=88)
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) = FKENER ( ONEONE, ONEONE )
96 EXMASS(3) = FKENER ( TWOTWO, ONEONE )
97 EXMASS(4) = FKENER ( THRTHR, ONEONE )
98 EXMASS(5) = FKENER ( THRTHR, TWOTWO )
99 EXMASS(6) = FKENER ( 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 = FKENER ( EIGEIG, FOUFOU ) - TWOTWO * EXMASS (6)
118 EMN = 1.D+03 * AMNEUT
119 EMH = ZMASS (2)
120 TMP16 = 16.D+00
121 UM = AMUMEV + FKENER ( 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 = FKENER (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 - FKENER ( 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 + FKENER ( 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)=FKENER(A-FLA(K),Z-FLZ(K))-FKENER(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 EXPSAS=MIN(EXPMAX,MAX(EXPMIN,S(J)-SES))
662 SAS = EXP (EXPSAS)
663 EXPSUS=MIN(EXPMAX,MAX(EXPMIN,-S(J)))
664 SUS = EXP (EXPSUS)
665 EYE1(J) = ( TWOTWO * S(J)**2 -SIXSIX * S(J)
666 & + SIXSIX + SUS * ( S(J)**2 - SIXSIX ) )
667 & / ( EIGEIG * SMALLA(J)**2 )
668 IF ( J .EQ. 1 ) THEN
669 EYE0(J) = ( S(J) - ONEONE + SUS ) / ( TWOTWO*SMALLA(J) )
670* Standard calculation
671 R (J) = RMASS(MM)**2 * ALPH(MM) * ( EYE1(J) + BET(MM)
672 & * EYE0(J) ) * SAS
673 ELSE
674 R (J) = CCOUL(J) * RMASS(MM)**2 * EYE1(J) * SAS
675 END IF
676 R (J) = MAX ( ZERZER, R (J) )
677 SIGMA = SIGMA + R (J)
678 END IF
679 1131 CONTINUE
680 NCOUNT = 0
681 6202 CONTINUE
682 IF(SIGMA)9,9,10
683 9 CONTINUE
684 DO 6100 J = 1,6
685 IF(JA-IA(J))6100,6101,6100
686 6101 CONTINUE
687 IF(JZ-IZ(J))6100,6102,6100
688 6100 CONTINUE
689 GO TO 6099
690 6102 CONTINUE
691 IF ( U .GT. ANGLGB ) GO TO 1000
692 JEMISS = J
693C*****STORE,RESIDUAL NUC IS OF EMITTED PARTICLE TYPE
694* If we are here this means that the residual nucleus is equal to
695* one of the six emitted particle (the j-th one). So give to it
696* all the energy, score it and return with 0 recoil and excitation
697* energy for the residual nucleus
698 EPS = EREC
699 NPART(JEMISS) = NPART(JEMISS)+1
700 ITEMP=NPART(JEMISS)
701 NRNEEP = NRNEEP + 1
702 SMOM1(JEMISS) = SMOM1(JEMISS) + EPS
703 ITEMP=NPART(JEMISS)
704 EPART(ITEMP,JEMISS)=EPS
705 COSEVP(1,ITEMP,JEMISS) = COSLBR(1)
706 COSEVP(2,ITEMP,JEMISS) = COSLBR(2)
707 COSEVP(3,ITEMP,JEMISS) = COSLBR(3)
708 GO TO 8002
709*
710*-->-->-->-->--> go to return
711 6099 CONTINUE
712 IF(JA-8)72,51,72
713 51 CONTINUE
714 IF(JZ-4)72,1224,72
715* Come here for a "normal" step
716 10 CONTINUE
717 LOPPAR = .FALSE.
718 CALL GRNDM(RNDM,1)
719 URAN=RNDM(1)*SIGMA
720 SUM = ZERZER
721 DO 16 J=1,6
722 K = J
723 SUM = R(J)+SUM
724 IF ( SUM - URAN .GT. 0.D+00 ) GO TO 17
725 16 CONTINUE
726 17 CONTINUE
727 JEMISS=K
728 NPART(JEMISS) = NPART (JEMISS) + 1
729 JS = SOS (JEMISS) + ONEONE
730 IF(JS-1000)4344,4345,4345
731 4345 CONTINUE
732 RATIO2=(S(JEMISS)**3-6.D0*S(JEMISS)**2+15.D0*
733 &S(JEMISS)-15.D0)/((2.D0*S(JEMISS)**2-6.D0*S(JEMISS)+6.D0)*SMALLA
734 &(JEMISS))
735 GO TO 4347
736 4344 CONTINUE
737 RATIO2=(P2(JS)+(P2(JS+1)-P2(JS))*
738 &(SOS(JEMISS)-STRUN(JEMISS)))/SMALLA(JEMISS)
739 4347 CONTINUE
740 EPSAV=RATIO2*TWOTWO
741* +-------------------------------------------------------------------*
742* | Neutron channel selected:
743 IF (JEMISS .EQ. 1) THEN
744 MM=JA-IA(J)
745 EPSAV=(EPSAV+BET(MM))/(ONEONE+BET(MM)*EYE0(JEMISS)
746 & /EYE1(JEMISS))
747* | +----------------------------------------------------------------*
748* | | Compute the fission width relative to the neutron one:
749* | | this part is taken from subroutine EMIT of LAHET
750 IF ( IFISS .GT. 0 .AND. JZ .GE. 71 .AND. .NOT. FISINH ) THEN
751* | | +-------------------------------------------------------------*
752* | | | Compute the correction factor for the fission width:
753 IF ( JZ .GT. 88 ) THEN
754 AGOES = ONEONE
755* | | |
756* | | +-------------------------------------------------------------*
757* | | |
758 ELSE
759 AGOES = MAX ( ONEONE, ( U-SEVSEV ) / ( EPSAV+SEVSEV ) )
760 END IF
761* | | |
762* | | +-------------------------------------------------------------*
763* | | Finally this is the relative fission width:
764* | | This is : Probfs = 1 / ( 1 + G_n / G_f )
765 PROBFS = FPROB ( Z, A, U ) / AGOES
766* | | +-------------------------------------------------------------*
767* | | | Check if it will be fission:
768 CALL GRNDM(RNDM,1)
769 IF ( RNDM (1) .LT. PROBFS ) THEN
770 FISINH = .TRUE.
771 KFISS = 1
772* | | | Undo the counting of the neutron evaporation
773 NPART (JEMISS) = NPART (JEMISS) -1
774 GO TO 9260
775 END IF
776* | | |
777* | | +-------------------------------------------------------------*
778 END IF
779* | |
780* | +----------------------------------------------------------------*
781 END IF
782* |
783* +-------------------------------------------------------------------*
784 20 CONTINUE
785 CALL GRNDM(RNDM,2)
786 E1=-HLFHLF*LOG(RNDM(1))
787 E2=-HLFHLF*LOG(RNDM(2))
788* Eps should be the total kinetic energy in the CMS frame
789* Standard calculation:
790 EPS=(E1+E2)*EPSAV+THRESH(JEMISS)-Q(JEMISS)
791 AR = A - IA(JEMISS)
792 ZR = Z - IZ(JEMISS)
793* The CMS energy is updated
794 IMASS = NINT (AR)
795 IF ( IMASS .EQ. 8 .AND. NINT (ZR) .EQ. 4 ) THEN
796 UNEW = U - EPS - Q(JEMISS)
797 UMAX = U - THRESH(JEMISS)
798 IF ( UNEW .GT. 6.D+00 ) THEN
799 UMIN = 6.D+00
800 ELSE IF ( UNEW .GT. 4.47D+00 .AND. UMAX .GT. 6.D+00 ) THEN
801 UMIN = 4.47D+00
802 UNEW = 6.D+00
803 ELSE IF ( UNEW .GT. 1.47D+00 .AND. UMAX .GT. 2.94D+00 ) THEN
804 UMIN = 1.47D+00
805 UNEW = 2.94D+00
806 ELSE
807 UMIN = -0.1D+00
808 UNEW = ANGLGB * 0.1D+00
809 END IF
810 ELSE IF ( IMASS .LE. 4 ) THEN
811 IPRO = NINT ( ZR )
812 INEU = IMASS - IPRO
813 IF ( IMASS .EQ. 1 ) THEN
814* Be sure that residual neutrons or protons are not left excited
815 UMIN = 0.D+00
816 UNEW = 0.D+00
817 EPS = U - Q(JEMISS)
818 ELSE IF ( IPRO .EQ. 0 .OR. INEU .EQ. 0 ) THEN
819* Ipro protons or ineu neutrons arrived here!
820 UMIN = CORRRR(JEMISS)
821 UNEW = U - EPS - Q(JEMISS)
822 ELSE IF ( IMASS .LE. 2 ) THEN
823* Be sure that residual deuterons are not left excited!
824 UMIN = 0.D+00
825 UNEW = 0.D+00
826 EPS = U - Q(JEMISS)
827 ELSE IF ( ABS ( INEU - IPRO ) .LE. 1 ) THEN
828* For the moment also residual 3-H, 3-He and 4-He are not left
829* excited !
830 UMIN = 0.D+00
831 UNEW = 0.D+00
832 EPS = U - Q(JEMISS)
833 ELSE
834 UMIN = CORRRR(JEMISS)
835 UNEW = U - EPS - Q(JEMISS)
836 END IF
837 ELSE
838 UMIN = CORRRR(JEMISS)
839 UNEW = U - EPS - Q(JEMISS)
840 END IF
841* Standard calculation
842 IF(UNEW-UMIN)6200,6220,6220
843 6220 CONTINUE
844*or RNMASS = AR * AMUMEV + FKENER (AR,ZR)
845 RNMASS = AR * AMUMEV + FKENER (AR,ZR) + UNEW
846 UMIN2 = ( RNMASS + ZMASS (JEMISS) )**2
847 IF ( UMIN2 .GE. UMO2 ) THEN
848 GO TO 6200
849 END IF
850 U = UNEW
851* C(i) are the direction cosines of the evaporated particle in the CMS
852* frame, of course - C(i) are the ones of the residual nucleus
853 CALL RACO(C(1),C(2),C(3))
854* Erncm, Eepcm are the total energies of the residual nucleus and
855* of the evaporated particle in the CMS frame
856 ERNCM = 0.5D+00 * ( UMO2 + RNMASS**2 - Z2MASS (JEMISS) ) / UMO
857 EEPCM = UMO - ERNCM
858 PCMS = SQRT ( EEPCM**2 - Z2MASS (JEMISS) )
859* Now we perform the Lorentz transformation back to the original
860* frame (lab frame)
861* First the evaporated particle:
862 ETAX = ETACM * COSLBR (1)
863 ETAY = ETACM * COSLBR (2)
864 ETAZ = ETACM * COSLBR (3)
865 PCMSX = PCMS * C (1)
866 PCMSY = PCMS * C (2)
867 PCMSZ = PCMS * C (3)
868 ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
869 EPS = GAMCM * EEPCM + ETAPCM - ZMASS (JEMISS)
870 PHELP = ETAPCM / (GAMCM + 1.D0) + EEPCM
871 PLBPX = PCMSX + ETAX * PHELP
872 PLBPY = PCMSY + ETAY * PHELP
873 PLBPZ = PCMSZ + ETAZ * PHELP
874 PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
875 COSLBP (1) = PLBPX / PHELP
876 COSLBP (2) = PLBPY / PHELP
877 COSLBP (3) = PLBPZ / PHELP
878* Then the residual nucleus ( for it c (i) --> - c (i) ):
879 EREC = GAMCM * ERNCM - ETAPCM - RNMASS
880 EKRES = 1.D-03 * EREC
881 PHELP = - ETAPCM / (GAMCM + 1.D0) + ERNCM
882 PXRES = 1.D-03 * ( - PCMSX + ETAX * PHELP )
883 PYRES = 1.D-03 * ( - PCMSY + ETAY * PHELP )
884 PZRES = 1.D-03 * ( - PCMSZ + ETAZ * PHELP )
885 P2RES = PXRES * PXRES + PYRES * PYRES + PZRES * PZRES
886 PTRES = SQRT (P2RES)
887 COSLBR (1) = PXRES / PTRES
888 COSLBR (2) = PYRES / PTRES
889 COSLBR (3) = PZRES / PTRES
890* Check energy and momentum conservation !!
891 IF (EREC .LE. 0.D+00) THEN
892 PTRES = 0.D+00
893 EREC = 0.D+00
894 END IF
895* Umo is the invariant mass of the system!!
896 UMO = RNMASS
897 UMO2 = UMO * UMO
898 ELBTOT = RNMASS + EREC
899 GAMCM = ELBTOT / RNMASS
900 ETACM = 1.D+03 * PTRES / RNMASS
901 GO TO 76
902 6200 CONTINUE
903 NCOUNT = NCOUNT + 1
904 IF ( NCOUNT .LE. 10 ) GO TO 20
905 SIGMA = SIGMA - R(JEMISS)
906* if we are here we have sampled for > 10 times a negative energy Unew
907 NPART(JEMISS)=NPART(JEMISS)-1
908 R(JEMISS) = 0.D0
909 NCOUNT = 0
910 GO TO 6202
911 76 CONTINUE
912 JAT=JA-IA(JEMISS)
913 JZT=JZ-IZ(JEMISS)
914 IF(JAT-JZT)9,9,77
915 77 CONTINUE
916 JA=JAT
917 JZ=JZT
918C*****STORE,END OF NORMAL CYCLE
919 SMOM1(JEMISS)=SMOM1(JEMISS)+EPS
920 ITEMP=NPART(JEMISS)
921 EPART(ITEMP,JEMISS)=EPS
922 COSEVP(1,ITEMP,JEMISS)=COSLBP(1)
923 COSEVP(2,ITEMP,JEMISS)=COSLBP(2)
924 COSEVP(3,ITEMP,JEMISS)=COSLBP(3)
925* The following card switch to the rough splitting treatment
926 IF (JA .LE. 2) GO TO 1000
927 IF(JA-8)23,1223,23
928 1223 CONTINUE
929 IF(JZ-4)23,1224,23
930* If we are here the residual nucleus is a 8-Be one, break it into
931* two alphas with all the available energy (U plus the Q of the breakup)
932* , score them and return with 0 recoil and excitation energy
933 1224 CONTINUE
934 LOPPAR = .FALSE.
935 IF(U)1228,1229,1229
936 1228 CONTINUE
937 EPS=0.D0
938 GO TO 1230
939 1229 CONTINUE
940 1230 CONTINUE
941 NBE8=NBE8+1
942* C(i) are the direction cosines of the first alpha in the CMS
943* frame, of course - C(i) are the ones of the other
944 CALL RACO(C(1),C(2),C(3))
945* Ecms is the total energy of the alphas in the CMS frame
946 ECMS = 0.5D+00 * UMO
947 PCMS = SQRT ( ECMS**2 - Z2MASS (6) )
948* Now we perform the Lorentz transformation back to the original
949* frame (lab frame)
950* First alpha:
951 ETAX = ETACM * COSLBR (1)
952 ETAY = ETACM * COSLBR (2)
953 ETAZ = ETACM * COSLBR (3)
954 PCMSX = PCMS * C (1)
955 PCMSY = PCMS * C (2)
956 PCMSZ = PCMS * C (3)
957 ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
958 EPS = GAMCM * ECMS + ETAPCM - ZMASS (6)
959 PHELP = ETAPCM / (GAMCM + 1.D0) + ECMS
960 PLBPX = PCMSX + ETAX * PHELP
961 PLBPY = PCMSY + ETAY * PHELP
962 PLBPZ = PCMSZ + ETAZ * PHELP
963 PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
964* Store the first alpha!!
965 SMOM1(6) = SMOM1(6) + EPS
966 NPART(6) = NPART(6) + 1
967 ITEMP = NPART(6)
968 EPART (ITEMP,6) = EPS
969 COSEVP (1,ITEMP,6) = PLBPX / PHELP
970 COSEVP (2,ITEMP,6) = PLBPY / PHELP
971 COSEVP (3,ITEMP,6) = PLBPZ / PHELP
972* Then the second alpha ( for it c (i) --> - c (i) ):
973 EPS = GAMCM * ECMS - ETAPCM - ZMASS (6)
974 PHELP = - ETAPCM / (GAMCM + 1.D0) + ECMS
975 PLBPX = - PCMSX + ETAX * PHELP
976 PLBPY = - PCMSY + ETAY * PHELP
977 PLBPZ = - PCMSZ + ETAZ * PHELP
978 PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
979* Store the second alpha !!
980 SMOM1(6) = SMOM1(6) + EPS
981 NPART(6) = NPART(6) + 1
982 ITEMP = NPART(6)
983 EPART (ITEMP,6) = EPS
984 COSEVP (1,ITEMP,6) = PLBPX / PHELP
985 COSEVP (2,ITEMP,6) = PLBPY / PHELP
986 COSEVP (3,ITEMP,6) = PLBPZ / PHELP
987 8002 CONTINUE
988 LOPPAR = .FALSE.
989 EREC = ZERZER
990 U = ZERZER
991 EKRES = ZERZER
992 PTRES = ZERZER
993 72 CONTINUE
994 HEVSUM=SMOM1(3)+SMOM1(5)+SMOM1(6)+SMOM1(4)
995 RETURN
996*.......................................................................
997*///// RAL FISSION ROUTINE /////
998 9260 CONTINUE
999* +-------------------------------------------------------------------*
1000* | Record the direction cosines of the fissioning nucleus
1001 DO 9270 I=1,3
1002 COSLF0 (I) = COSLBR (I)
1003 9270 CONTINUE
1004* |
1005* +-------------------------------------------------------------------*
1006 CALL FISFRA ( JA, JZ, U, EREC, UMO, GAMCM, ETACM )
1007* +-------------------------------------------------------------------*
1008* | Check for fission failures!!
1009 IF ( .NOT. FISINH ) THEN
1010 PENBAR = .FALSE.
1011 GO TO 23
1012 END IF
1013* |
1014* +-------------------------------------------------------------------*
1015* Do not pick up the fission fragments, rather go back to Evevap
1016 HEVSUM=SMOM1(3)+SMOM1(5)+SMOM1(6)+SMOM1(4)
1017 RETURN
1018 END