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