]>
Commit | Line | Data |
---|---|---|
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 | * | |
63 | C--------------------------------------------------------------------- | |
64 | C SUBNAME = DRES --- EVAPORATION | |
65 | C EVAPORATION DATA SHOUD BE READ ON INPUT STAGE | |
66 | C--------------------------------------------------------------------- | |
67 | C*****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 | |
88 | C-------------------------------------- 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) | |
114 | 1234 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 | * +-------------------------------------------------------------------* | |
127 | C --------------------------------- 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) | |
584 | C 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 | |
691 | C*****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 | |
916 | C*****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 |