source: trunk/source/processes/hadronic/models/incl/src/G4AblaFission.cc@ 1201

Last change on this file since 1201 was 968, checked in by garnier, 17 years ago

fichier ajoutes
File size: 46.2 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: G4AblaFission.cc,v 1.3 2008/11/06 08:42:00 gcosmo Exp $
27// Translation of INCL4.2/ABLA V3
28// Pekka Kaitaniemi, HIP (translation)
29// Christelle Schmidt, IPNL (fission code)
30// Alain Boudard, CEA (contact person INCL/ABLA)
31// Aatos Heikkinen, HIP (project coordination)
32
33#include "G4AblaFission.hh"
34#include <time.h>
35
36G4AblaFission::G4AblaFission()
37{
38}
39
40G4AblaFission::G4AblaFission(G4Hazard *hzr, G4InclRandomInterface *rndm)
41{
42 hazard = hzr;
43 randomGenerator = rndm;
44 setAboutString("Fission model: Based on ABLA V3 and SimFis3");
45}
46
47G4AblaFission::~G4AblaFission()
48{
49}
50
51void G4AblaFission::doFission(G4double &A, G4double &Z, G4double &E,
52 G4double &A1, G4double &Z1, G4double &E1, G4double &K1,
53 G4double &A2, G4double &Z2, G4double &E2, G4double &K2)
54{
55 fissionDistri(A,Z,E,A1,Z1,E1,K1,A2,Z2,E2,K2);
56}
57
58void G4AblaFission::even_odd(G4double r_origin,G4double r_even_odd,G4int &i_out)
59{
60 // Procedure to calculate I_OUT from R_IN in a way that
61 // on the average a flat distribution in R_IN results in a
62 // fluctuating distribution in I_OUT with an even-odd effect as
63 // given by R_EVEN_ODD
64
65 // /* ------------------------------------------------------------ */
66 // /* EXAMPLES : */
67 // /* ------------------------------------------------------------ */
68 // /* If R_EVEN_ODD = 0 : */
69 // /* CEIL(R_IN) ---- */
70 // /* */
71 // /* R_IN -> */
72 // /* (somewhere in between CEIL(R_IN) and FLOOR(R_IN)) */ */
73 // /* */
74 // /* FLOOR(R_IN) ---- --> I_OUT */
75 // /* ------------------------------------------------------------ */
76 // /* If R_EVEN_ODD > 0 : */
77 // /* The interval for the above treatment is */
78 // /* larger for FLOOR(R_IN) = even and */
79 // /* smaller for FLOOR(R_IN) = odd */
80 // /* For R_EVEN_ODD < 0 : just opposite treatment */
81 // /* ------------------------------------------------------------ */
82
83 // /* ------------------------------------------------------------ */
84 // /* On input: R_ORIGIN nuclear charge (real number) */
85 // /* R_EVEN_ODD requested even-odd effect */
86 // /* Intermediate quantity: R_IN = R_ORIGIN + 0.5 */
87 // /* On output: I_OUT nuclear charge (integer) */
88 // /* ------------------------------------------------------------ */
89
90 // G4double R_ORIGIN,R_IN,R_EVEN_ODD,R_REST,R_HELP;
91 G4double r_in = 0.0, r_rest = 0.0, r_help = 0.0;
92 G4double r_floor = 0.0;
93 G4double r_middle = 0.0;
94 // G4int I_OUT,N_FLOOR;
95 G4int n_floor = 0;
96
97 r_in = r_origin + 0.5;
98 r_floor = (float)((int)(r_in));
99 if (r_even_odd < 0.001) {
100 i_out = (int)(r_floor);
101 }
102 else {
103 r_rest = r_in - r_floor;
104 r_middle = r_floor + 0.5;
105 n_floor = (int)(r_floor);
106 if (n_floor%2 == 0) {
107 // even before modif.
108 r_help = r_middle + (r_rest - 0.5) * (1.0 - r_even_odd);
109 }
110 else {
111 // odd before modification
112 r_help = r_middle + (r_rest - 0.5) * (1.0 + r_even_odd);
113 }
114 i_out = (int)(r_help);
115 }
116}
117
118G4double G4AblaFission::umass(G4double z,G4double n,G4double beta)
119{
120 // liquid-drop mass, Myers & Swiatecki, Lysekil, 1967
121 // pure liquid drop, without pairing and shell effects
122
123 // On input: Z nuclear charge of nucleus
124 // N number of neutrons in nucleus
125 // beta deformation of nucleus
126 // On output: binding energy of nucleus
127
128 G4double a = 0.0, umass = 0.0;
129 G4double alpha = 0.0;
130 G4double xcom = 0.0, xvs = 0.0, xe = 0.0;
131 const G4double pi = 3.1416;
132
133 a = n + z;
134 alpha = ( std::sqrt(5.0/(4.0*pi)) ) * beta;
135
136 xcom = 1.0 - 1.7826 * ((a - 2.0*z)/a)*((a - 2.0*z)/a);
137 // factor for asymmetry dependence of surface and volume term
138 xvs = - xcom * ( 15.4941 * a -
139 17.9439 * std::pow(a,2.0/3.0) * (1.0+0.4*alpha*alpha) );
140 // sum of volume and surface energy
141 xe = z*z * (0.7053/(std::pow(a,1.0/3.0)) * (1.0-0.2*alpha*alpha) - 1.1529/a);
142 umass = xvs + xe;
143
144 return umass;
145}
146
147G4double G4AblaFission::ecoul(G4double z1,G4double n1,G4double beta1,G4double z2,G4double n2,G4double beta2,G4double d)
148{
149 // Coulomb potential between two nuclei
150 // surfaces are in a distance of d
151 // in a tip to tip configuration
152
153 // approximate formulation
154 // On input: Z1 nuclear charge of first nucleus
155 // N1 number of neutrons in first nucleus
156 // beta1 deformation of first nucleus
157 // Z2 nuclear charge of second nucleus
158 // N2 number of neutrons in second nucleus
159 // beta2 deformation of second nucleus
160 // d distance of surfaces of the nuclei
161
162 // G4double Z1,N1,beta1,Z2,N2,beta2,d,ecoul;
163 G4double ecoul = 0;
164 G4double dtot = 0;
165 const G4double r0 = 1.16;
166
167 dtot = r0 * ( std::pow((z1+n1),1.0/3.0) * (1.0+0.6666*beta1)
168 + std::pow((z2+n2),1.0/3.0) * (1.0+0.6666*beta2) ) + d;
169 ecoul = z1 * z2 * 1.44 / dtot;
170
171 return ecoul;
172}
173
174void G4AblaFission::fissionDistri(G4double &a,G4double &z,G4double &e,
175 G4double &a1,G4double &z1,G4double &e1,G4double &v1,
176 G4double &a2,G4double &z2,G4double &e2,G4double &v2)
177{
178 // G4cout <<"Fission: a = " << a << " z = " << z << " e = " << e << G4endl;
179 // On input: A, Z, E (mass, atomic number and exc. energy of compound nucleus
180 // before fission)
181 // On output: Ai, Zi, Ei (mass, atomic number and exc. energy of fragment 1 and 2
182 // after fission)
183
184 // Additionally calculated but not put in the parameter list:
185 // Kinetic energy of prefragments EkinR1, EkinR2
186
187 // Translation of SIMFIS18.PLI (KHS, 2.1.2001)
188
189 // This program calculates isotopic distributions of fission fragments
190 // with a semiempirical model
191 // Copy from SIMFIS3, KHS, 8. February 1995
192 // Modifications made by Jose Benlliure and KHS in August 1996
193 // Energy counted from lowest barrier (J. Benlliure, KHS 1997)
194 // Some bugs corrected (J. Benlliure, KHS 1997)
195 // Version used for thesis S. Steinhaueser (August 1997)
196 // (Curvature of LD potential increased by factor of 2!)
197
198 // Weiter veraendert mit der Absicht, eine Version zu erhalten, die
199 // derjenigen entspricht, die von J. Benlliure et al.
200 // in Nucl. Phys. A 628 (1998) 458 verwendet wurde,
201 // allerdings ohne volle Neutronenabdampfung.
202
203 // The excitation energy was calculate now for each fission channel
204 // separately. The dissipation from saddle to scission was taken from
205 // systematics, the deformation energy at scission considers the shell
206 // effects in a simplified way, and the fluctuation is included.
207 // KHS, April 1999
208
209 // The width in N/Z was carefully adapted to values given by Lang et al.
210
211 // The width and eventually a shift in N/Z (polarization) follows the
212 // following rules:
213
214 // The line N/Z following UCD has an angle of std::atan(Zcn/Ncn)
215 // to the horizontal axis on a chart of nuclides.
216 // (For 238U the angle is 32.2 deg.)
217
218 // The following relations hold: (from Armbruster)
219 //
220 // sigma(N) (A=const) = sigma(Z) (A=const)
221 // sigma(A) (N=const) = sigma(Z) (N=const)
222 // sigma(A) (Z=const) = sigma(N) (Z=const)
223 //
224 // From this we get:
225 // sigma(Z) (N=const) * N = sigma(N) (Z=const) * Z
226 // sigma(A) (Z=const) = sigma(Z) (A=const) * A/Z
227 // sigma(N) (Z=const) = sigma(Z) (A=const) * A/Z
228 // Z*sigma(N) (Z=const) = N*sigma(Z) (N=const) = A*sigma(Z) (A=const)
229
230 // Excitation energy now calculated above the lowest potential point
231 // Inclusion of a distribution of excitation energies
232
233 // Several modifications, starting from SIMFIS12: KHS November 2000
234 // This version seems to work quite well for 238U.
235 // The transition from symmetric to asymmetric fission around 226Th
236 // is reasonably well reproduced, although St. I is too strong and St. II
237 // is too weak. St. I and St. II are also weakly seen for 208Pb.
238
239 // Extensions for an event generator of fission events (21.11.2000,KHS)
240
241 // Defalt parameters (IPARS) rather carefully adjusted to
242 // pre-neutron mass distributions of Vives et al. (238U + n)
243 // Die Parameter Fgamma1 und Fgamma2 sind kleiner als die resultierenden
244 // Breiten der Massenverteilungen!!!
245 // Fgamma1 und Fgamma2 wurden angepaￜ, so daￜ
246 // Sigma-A(ST-I) = 3.3, Sigma-A(St-II) = 5.8 (nach Vives)
247
248 // Parameters of the model carefully adjusted by KHS (2.2.2001) to
249 // 238U + 208Pb, 1000 A MeV, Timo Enqvist et al.
250
251 G4double n = 0.0;
252 G4double nlight1 = 0.0, nlight2 = 0.0;
253 G4double aheavy1 = 0.0,alight1 = 0.0, aheavy2 = 0.0, alight2 = 0.0;
254 // G4double eheavy1 = 0.0, elight1 = 0.0, eheavy2 = 0.0, elight2 = 0.0;
255 G4double zheavy1_shell = 0.0, zheavy2_shell = 0.0;
256 G4double zlight1 = 0.0, zlight2 = 0.0;
257 G4double masscurv = 0.0;
258 G4double sasymm1 = 0.0, sasymm2 = 0.0, ssymm = 0.0, ysum = 0.0, yasymm = 0.0;
259 G4double ssymm_mode1 = 0.0, ssymm_mode2 = 0.0;
260 G4double cz_asymm1_saddle = 0.0, cz_asymm2_saddle = 0.0;
261 // Curvature at saddle, modified by ld-potential
262 G4double wzasymm1_saddle, wzasymm2_saddle, wzsymm_saddle = 0.0;
263 G4double wzasymm1_scission = 0.0, wzasymm2_scission = 0.0, wzsymm_scission = 0.0;
264 G4double wzasymm1 = 0.0, wzasymm2 = 0.0, wzsymm = 0.0;
265 G4double nlight1_eff = 0.0, nlight2_eff = 0.0;
266 G4int imode = 0;
267 G4double rmode = 0.0;
268 G4double z1mean = 0.0, z1width = 0.0;
269 // G4double Z1,Z2,N1R,N2R,A1R,A2R,N1,N2,A1,A2;
270 G4double n1r = 0.0, n2r = 0.0, n1 = 0.0, n2 = 0.0;
271
272 G4double zsymm = 0.0, nsymm = 0.0, asymm = 0.0;
273 G4double n1mean = 0.0, n1width = 0.0;
274 // effective shell effect at lowest barrier
275 // Excitation energy with respect to ld barrier
276 G4double re1 = 0.0, re2 = 0.0, re3 = 0.0;
277 G4double n1ucd = 0.0, n2ucd = 0.0;
278
279 // shift of most probable neutron number for given Z,
280 // according to polarization
281 G4int i_help = 0;
282
283 // /* Parameters of the semiempirical fission model */
284 G4double a_levdens = 0.0;
285 // /* level-density parameter */
286 G4double a_levdens_light1 = 0.0, a_levdens_light2 = 0.0;
287 G4double a_levdens_heavy1 = 0.0, a_levdens_heavy2 = 0.0;
288 const G4double r_null = 1.16;
289 // /* radius parameter */
290 G4double epsilon_1_saddle = 0.0, epsilon0_1_saddle = 0.0;
291 G4double epsilon_2_saddle = 0.0, epsilon0_2_saddle = 0.0, epsilon_symm_saddle = 0.0;
292 G4double epsilon_1_scission = 0.0, epsilon0_1_scission = 0.0;
293 G4double epsilon_2_scission = 0.0, epsilon0_2_scission = 0.0;
294 G4double epsilon_symm_scission = 0.0;
295 // /* modified energy */
296 G4double e_eff1_saddle = 0.0, e_eff2_saddle = 0.0;
297 G4double a1r;
298 G4int icz = 0, k = 0;
299
300 G4int i_inter = 0;
301 G4double ne_min = 0, ne_m1 = 0, ne_m2 = 0;
302 G4double ed1_low = 0.0, ed2_low = 0.0, ed1_high = 0.0, ed2_high = 0.0, ed1 = 0.0, ed2 = 0.0;
303 G4double atot = 0.0;
304
305 // Input parameters:
306 //OMMENT(Nuclear charge number);
307 // G4double Z;
308 //OMMENT(Nuclear mass number);
309 // G4double A;
310 //OMMENT(Excitation energy above fission barrier);
311 // G4double E;
312
313 // Model parameters:
314 //OMMENT(position of heavy peak valley 1);
315 const G4double nheavy1 = 82.0;
316 const G4double nheavy2 = 89.0;
317 //OMMENT(position of heavy peak valley 2);
318 const G4double e_crit = 5; // Critical pairing energy :::PK
319 //OMMENT(Shell effect for valley 1);
320 // Parameter (Delta_U2_shell = -3.2)
321 //OMMENT(I: used shell effect);
322 G4double delta_u1 = 0.0;
323 //omment(I: used shell effect);
324 G4double delta_u2 = 0.0;
325 const G4double delta_u1_shell = -2.5;
326 // Parameter (Delta_U1_shell = -2)
327 //OMMENT(Shell effect for valley 2);
328 const G4double delta_u2_shell = -5.5;
329 const G4double el = 30.0;
330 //OMMENT(Curvature of asymmetric valley 1);
331 const G4double cz_asymm1_shell = 0.7;
332 //OMMENT(Curvature of asymmetric valley 2);
333 const G4double cz_asymm2_shell = 0.08;
334 //OMMENT(Factor for width of distr. valley 1);
335 const G4double fwidth_asymm1 = 1.2;
336 //OMMENT(Factor for width of distr. valley 2);
337 const G4double fwidth_asymm2 = 1.0;
338 // Parameter (CZ_asymm2_scission = 0.12)
339 //OMMENT(Factor to gamma_heavy1);
340 const G4double fgamma1 = 1.0;
341 //OMMENT(I: fading of shells (general));
342 G4double gamma = 0.0;
343 //OMMENT(I: fading of shell 1);
344 G4double gamma_heavy1 = 0.0;
345 //OMMENT(I: fading of shell 2);
346 G4double gamma_heavy2 = 0.0;
347 //OMMENT(Zero-point energy at saddle);
348 const G4double e_zero_point = 0.5;
349 G4int i_eva = 0; // Calculate A = 1 or Aprime = 0 :::PK
350 //OMMENT(I: friction from saddle to scission);
351 G4double e_saddle_scission = 10.0;
352 //OMMENT(Friction factor);
353 const G4double friction_factor = 1.0;
354 //OMMENT(I: Internal counter for different modes); INIT(0,0,0)
355 // Integer*4 I_MODE(3)
356 //OMMENT(I: Yield of symmetric mode);
357 G4double ysymm = 0.0;
358 //OMMENT(I: Yield of asymmetric mode 1);
359 G4double yasymm1 = 0.0;
360 //OMMENT(I: Yield of asymmetric mode 2);
361 G4double yasymm2 = 0.0;
362 //OMMENT(I: Effective position of valley 1);
363 G4double nheavy1_eff = 0.0;
364 //OMMENT(I: position of heavy peak valley 1);
365 G4double zheavy1 = 0.0;
366 //omment(I: Effective position of valley 2);
367 G4double nheavy2_eff = 0.0;
368 //OMMENT(I: position of heavy peak valley 2);
369 G4double zheavy2 = 0.0;
370 //omment(I: Excitation energy above saddle 1);
371 G4double eexc1_saddle = 0.0;
372 //omment(I: Excitation energy above saddle 2);
373 G4double eexc2_saddle = 0.0;
374 //omment(I: Excitation energy above lowest saddle);
375 G4double eexc_max = 0.0;
376 //omment(I: Effective mass mode 1);
377 G4double aheavy1_mean = 0.0;
378 //omment(I: Effective mass mode 2);
379 G4double aheavy2_mean = 0.0;
380 //omment(I: Width of symmetric mode);
381 G4double wasymm_saddle = 0.0;
382 //OMMENT(I: Width of asymmetric mode 1);
383 G4double waheavy1_saddle = 0.0;
384 //OMMENT(I: Width of asymmetric mode 2);
385 G4double waheavy2_saddle = 0.0;
386 //omment(I: Width of symmetric mode);
387 G4double wasymm = 0.0;
388 //OMMENT(I: Width of asymmetric mode 1);
389 G4double waheavy1 = 0.0;
390 //OMMENT(I: Width of asymmetric mode 2);
391 G4double waheavy2 = 0.0;
392 //OMMENT(I: Even-odd effect in Z);
393 G4double r_e_o = 0.0;
394 G4double r_e_o_max = 0.0;
395 G4double e_pair = 0.0;
396 //OMMENT(I: Curveture of symmetric valley);
397 G4double cz_symm = 0.0;
398 //OMMENT(I: Curvature of mass distribution for fixed Z);
399 G4double cn = 0.0;
400 //OMMENT(=1: test output, =0: no test output);
401 const G4int itest = 0;
402 // G4double UMASS, ECOUL, reps1, reps2, rn1_pol;
403 G4double reps1 = 0.0, reps2 = 0.0, rn1_pol = 0.0;
404 // Float_t HAZ,GAUSSHAZ;
405 //G4int kkk = 0;
406 G4int kkk = 10;
407 G4double Bsym = 0.0;
408 G4double Basym_1 = 0.0;
409 G4double Basym_2 = 0.0;
410 G4int iz = 0;
411 // I_MODE = 0;
412
413 if(itest == 1) {
414 G4cout << " cn mass " << a << G4endl;
415 G4cout << " cn charge " << z << G4endl;
416 G4cout << " cn energy " << e << G4endl;
417 }
418
419 // /* average Z of asymmetric and symmetric components: */
420 n = a - z; /* neutron number of the fissioning nucleus */
421
422 k = 0;
423 icz = 0;
424 if ( (std::pow(z,2)/a < 25.0) || (n < nheavy2) || (e > 500.0) ) {
425 icz = -1;
426 // GOTO 1002;
427 goto milledeux;
428 }
429
430 nlight1 = n - nheavy1;
431 nlight2 = n - nheavy2;
432
433 // /* Polarisation assumed for standard I and standard II:
434 // Z - Zucd = cpol (for A = const);
435 // from this we get (see Armbruster)
436 // Z - Zucd = Acn/Ncn * cpol (for N = const) */
437
438 // zheavy1_shell = ((nheavy1/n) * z) - ((a/n) * cpol1); // Simfis18 PK:::
439 zheavy1_shell = ((nheavy1/n) * z) - 0.8; // Simfis3 PK:::
440 //zheavy2_shell = ((nheavy2/n) * z) - ((a/n) * cpol2); // Simfis18 PK:::
441 zheavy2_shell = ((nheavy2/n) * z) - 0.8; // Simfis3 PK:::
442
443 // p(zheavy1_shell, zheavy2_shell);
444
445 // e_saddle_scission =
446 // (-24.0 + 0.02227 * (std::pow(z,2))/(std::pow(a,0.33333)) ) * friction_factor;
447 e_saddle_scission = (3.535 * std::pow(z,2)/a - 121.1) * friction_factor; // PK:::
448
449 // /* Energy dissipated from saddle to scission */
450 // /* F. Rejmund et al., Nucl. Phys. A 678 (2000) 215, fig. 4 b */
451 // E_saddle_scission = DMAX1(0.,E_saddle_scission);
452 // Heavy Ion Induced Reactions, Schroeder W. ed., Harwood, 1986, 101
453 if (e_saddle_scission < 0.0) {
454 e_saddle_scission = 0.0;
455 }
456
457 // /* Semiempirical fission model: */
458
459 // /* Fit to experimental result on curvature of potential at saddle */
460 // /* reference: */
461 // /* IF Z**2/A < 33.15E0 THEN
462 // MassCurv = 30.5438538E0 - 4.00212049E0 * Z**2/A
463 // + 0.11983384E0 * Z**4 / (A**2) ;
464 // ELSE
465 // MassCurv = 10.E0 ** (7.16993332E0 - 0.26602401E0 * Z**2/A
466 // + 0.00283802E0 * Z**4 / (A**2)) ; */
467 // /* New parametrization of T. Enqvist according to Mulgin et al. 1998 NPA 640(1998) 375 */
468 if ( (std::pow(z,2))/a < 33.9186) {
469 masscurv = std::pow( 10.0,(-1.093364 + 0.082933 * (std::pow(z,2)/a)
470 - 0.0002602 * (std::pow(z,4)/std::pow(a,2))) );
471 } else {
472 masscurv = std::pow( 10.0,(3.053536 - 0.056477 * (std::pow(z,2)/a)
473 + 0.0002454 * (std::pow(z,4)/std::pow(a,2))) );
474 }
475
476 cz_symm = (8.0/std::pow(z,2)) * masscurv;
477
478 if(itest == 1) {
479 G4cout << "cz_symmetry= " << cz_symm << G4endl;
480 }
481
482 icz = 0;
483 if (cz_symm < 0) {
484 icz = -1;
485 // GOTO 1002;
486 goto milledeux;
487 }
488
489 // /* proton number in symmetric fission (centre) */
490 zsymm = z/2.0;
491 nsymm = n/2.0;
492 asymm = nsymm + zsymm;
493
494 zheavy1 = (cz_symm*zsymm + cz_asymm1_shell*zheavy1_shell)/(cz_symm + cz_asymm1_shell);
495 zheavy2 = (cz_symm*zsymm + cz_asymm2_shell*zheavy2_shell)/(cz_symm + cz_asymm2_shell);
496
497 // /* position of valley due to influence of liquid-drop potential */
498 nheavy1_eff = (zheavy1 + 0.8)*(n/z);
499 nheavy2_eff = (zheavy2 + 0.8)*(n/z);
500 nlight1_eff = n - nheavy1_eff;
501 nlight2_eff = n - nheavy2_eff;
502 // /* proton number of light fragments (centre) */
503 zlight1 = z - zheavy1;
504 // /* proton number of light fragments (centre) */
505 zlight2 = z - zheavy2;
506 aheavy1 = nheavy1_eff + zheavy1;
507 aheavy2 = nheavy2_eff + zheavy2;
508 aheavy1_mean = aheavy1;
509 aheavy2_mean = aheavy2;
510 alight1 = nlight1_eff + zlight1;
511 alight2 = nlight2_eff + zlight2;
512 // Eheavy1 = E * Aheavy1 / A
513 // Eheavy2 = E * Aheavy2 / A
514 // Elight1 = E * Alight1 / A
515 // Elight2 = E * Alight2 / A
516 a_levdens = a / 8.0;
517 a_levdens_heavy1 = aheavy1 / 8.0;
518 a_levdens_heavy2 = aheavy2 / 8.0;
519 a_levdens_light1 = alight1 / 8.0;
520 a_levdens_light2 = alight2 / 8.0;
521 gamma = a_levdens / (0.4 * (std::pow(a,1.3333)) );
522 gamma_heavy1 = ( a_levdens_heavy1 / (0.4 * (std::pow(aheavy1,1.3333)) ) ) * fgamma1;
523 gamma_heavy2 = a_levdens_heavy2 / (0.4 * (std::pow(aheavy2,1.3333)) );
524
525 cz_asymm1_saddle = cz_asymm1_shell + cz_symm;
526 cz_asymm2_saddle = cz_asymm2_shell + cz_symm;
527
528 // Up to here: Ok! Checked CS 10/10/05
529
530 cn = umass(zsymm,(nsymm+1.),0.0) + umass(zsymm,(nsymm-1.),0.0)
531 + 1.44 * (std::pow(zsymm,2))/
532 ( (std::pow(r_null,2)) *
533 ( std::pow((asymm+1.0),1.0/3.0) + std::pow((asymm-1.0),1.0/3.0) ) *
534 ( std::pow((asymm+1.0),1.0/3.0) + std::pow((asymm-1.0),1.0/3.0) ) )
535 - 2.0 * umass(zsymm,nsymm,0.0)
536 - 1.44 * (std::pow(zsymm,2))/
537 ( ( 2.0 * r_null * (std::pow(asymm,1.0/3.0)) ) *
538 ( 2.0 * r_null * (std::pow(asymm,1.0/3.0)) ) );
539
540 // /* shell effect in valley of mode 1 */
541 delta_u1 = delta_u1_shell + (std::pow((zheavy1_shell-zheavy1),2))*cz_asymm1_shell;
542 // /* shell effect in valley of mode 2 */
543 delta_u2 = delta_u2_shell + (std::pow((zheavy2_shell-zheavy2),2))*cz_asymm2_shell;
544
545 Bsym = 0.0;
546 Basym_1 = Bsym + std::pow((zheavy1-zsymm), 2) * cz_symm + delta_u1;
547 Basym_2 = Bsym + std::pow((zheavy2-zsymm), 2) * cz_symm + delta_u2;
548 if(Bsym < Basym_1 && Bsym < Basym_2) {
549 // Excitation energies at the saddle point
550 // without and with shell effect
551 epsilon0_1_saddle = (e + e_zero_point - std::pow((zheavy1 - zsymm), 2) * cz_symm);
552 epsilon0_2_saddle = (e + e_zero_point - std::pow((zheavy2 - zsymm), 2) * cz_symm);
553
554 epsilon_1_saddle = epsilon0_1_saddle - delta_u1;
555 epsilon_2_saddle = epsilon0_2_saddle - delta_u2;
556
557 epsilon_symm_saddle = e + e_zero_point;
558 eexc1_saddle = epsilon_1_saddle;
559 eexc2_saddle = epsilon_2_saddle;
560
561 // Excitation energies at the scission point
562 // heavy fragment without and with shell effect
563 epsilon0_1_scission = (e + e_saddle_scission - std::pow((zheavy1 - zsymm), 2) * cz_symm) * aheavy1/a;
564 epsilon_1_scission = epsilon0_1_scission - delta_u1*(aheavy1/a);
565
566 epsilon0_2_scission = (e + e_saddle_scission - std::pow((zheavy2 - zsymm), 2) * cz_symm) * aheavy2/a;
567 epsilon_2_scission = epsilon0_2_scission - delta_u2*(aheavy2/a);
568
569 epsilon_symm_scission = e + e_saddle_scission;
570 } else if (Basym_1 < Bsym && Basym_1 < Basym_2) {
571 // Excitation energies at the saddle point
572 // without and with shell effect
573 epsilon_symm_saddle = (e + e_zero_point + delta_u1 + std::pow((zheavy1-zsymm), 2) * cz_symm);
574 epsilon0_2_saddle = (epsilon_symm_saddle - std::pow((zheavy2-zsymm), 2) * cz_symm);
575 epsilon_2_saddle = epsilon0_2_saddle - delta_u2;
576 epsilon0_1_saddle = e + e_zero_point + delta_u1;
577 epsilon_1_saddle = e + e_zero_point;
578 eexc1_saddle = epsilon_1_saddle;
579 eexc2_saddle = epsilon_2_saddle;
580
581 // Excitation energies at the scission point
582 // heavy fragment without and with shell effect
583 epsilon_symm_scission = (e + e_saddle_scission + std::pow((zheavy1-zsymm), 2) * cz_symm + delta_u1);
584 epsilon0_2_scission = (epsilon_symm_scission - std::pow((zheavy2-zsymm), 2) * cz_symm) * aheavy2/a;
585 epsilon_2_scission = epsilon0_2_scission - delta_u2*aheavy2/a;
586 epsilon0_1_scission = (e + e_saddle_scission + delta_u1) * aheavy1/a;
587 epsilon_1_scission = (e + e_saddle_scission) * aheavy1/a;
588 } else if (Basym_2 < Bsym && Basym_2 < Basym_1) {
589 // Excitation energies at the saddle point
590 // without and with shell effect
591 epsilon_symm_saddle = (e + e_zero_point + delta_u2 + std::pow((zheavy2-zsymm), 2) * cz_symm);
592 epsilon0_1_saddle = (epsilon_symm_saddle - std::pow((zheavy1-zsymm), 2) * cz_symm);
593 epsilon_1_saddle = epsilon0_1_saddle - delta_u1;
594 epsilon0_2_saddle = e + e_zero_point + delta_u2;
595 epsilon_2_saddle = e + e_zero_point;
596 eexc1_saddle = epsilon_1_saddle;
597 eexc2_saddle = epsilon_2_saddle;
598
599 // Excitation energies at the scission point
600 // heavy fragment without and with shell effect
601 epsilon_symm_scission = (e + e_saddle_scission + std::pow((zheavy2-zsymm), 2) * cz_symm + delta_u2);
602 epsilon0_1_scission = (epsilon_symm_scission - std::pow((zheavy1-zsymm), 2) * cz_symm) * aheavy1/a;
603 epsilon_1_scission = epsilon0_1_scission - delta_u1*aheavy1/a;
604 epsilon0_2_scission = (e + e_saddle_scission + delta_u2) * aheavy2/a;
605 epsilon_2_scission = (e + e_saddle_scission) * aheavy2/a;
606
607 } else {
608 G4cout <<"G4AblaFission: " << G4endl;
609 }
610 if(epsilon_1_saddle < 0.0) epsilon_1_saddle = 0.0;
611 if(epsilon_2_saddle < 0.0) epsilon_2_saddle = 0.0;
612 if(epsilon0_1_saddle < 0.0) epsilon0_1_saddle = 0.0;
613 if(epsilon0_2_saddle < 0.0) epsilon0_2_saddle = 0.0;
614 if(epsilon_symm_saddle < 0.0) epsilon_symm_saddle = 0.0;
615
616 if(epsilon_1_scission < 0.0) epsilon_1_scission = 0.0;
617 if(epsilon_2_scission < 0.0) epsilon_2_scission = 0.0;
618 if(epsilon0_1_scission < 0.0) epsilon0_1_scission = 0.0;
619 if(epsilon0_2_scission < 0.0) epsilon0_2_scission = 0.0;
620 if(epsilon_symm_scission < 0.0) epsilon_symm_scission = 0.0;
621
622 if(itest == 1) {
623 G4cout <<"E, E1, E2, Es" << e << epsilon_1_saddle << epsilon_2_saddle << epsilon_symm_saddle << G4endl;
624 }
625
626 e_eff1_saddle = epsilon0_1_saddle - delta_u1 * (std::exp((-epsilon_1_saddle*gamma)));
627
628 if (e_eff1_saddle > 0.0) {
629 wzasymm1_saddle = std::sqrt( (0.5) *
630 (std::sqrt(1.0/a_levdens*e_eff1_saddle)) /
631 (cz_asymm1_shell * std::exp((-epsilon_1_saddle*gamma)) + cz_symm) );
632 } else {
633 wzasymm1_saddle = 1.0;
634 }
635
636 e_eff2_saddle = epsilon0_2_saddle - delta_u2 * std::exp((-epsilon_2_saddle*gamma));
637 if (e_eff2_saddle > 0.0) {
638 wzasymm2_saddle = std::sqrt( (0.5 *
639 (std::sqrt(1.0/a_levdens*e_eff2_saddle)) /
640 (cz_asymm2_shell * std::exp((-epsilon_2_saddle*gamma)) + cz_symm) ) );
641 } else {
642 wzasymm2_saddle = 1.0;
643 }
644
645 if(e - e_zero_point > 0.0) {
646 wzsymm_saddle = std::sqrt( (0.5 *
647 (std::sqrt(1.0/a_levdens*(e+epsilon_symm_saddle))) / cz_symm ) );
648 } else {
649 wzsymm_saddle = 1.0;
650 }
651
652// if (itest == 1) {
653// G4cout << "wz1(saddle) = " << wzasymm1_saddle << G4endl;
654// G4cout << "wz2(saddle) = " << wzasymm2_saddle << G4endl;
655// G4cout << "wzsymm(saddle) = " << wzsymm_saddle << G4endl;
656// }
657
658 // /* Calculate widhts at the scission point: */
659 // /* fits of ref. Beizin 1991 (Plots brought to GSI by Sergei Zhdanov) */
660
661 wzsymm_scission = wzsymm_saddle;
662
663 if (e_saddle_scission == 0.0) {
664 wzasymm1_scission = wzasymm1_saddle;
665 wzasymm2_scission = wzasymm2_saddle;
666 } else {
667 if (nheavy1_eff > 75.0) {
668 wzasymm1_scission = (std::sqrt(21.0)) * z/a;
669 double RR = (70.0-28.0)/3.0*(z*z/a-35.0)+28.0;
670 if(RR > 0.0) {
671 wzasymm2_scission = std::sqrt(RR)*(z/a);
672 } else {
673 wzasymm2_scission = 0.0;
674 }
675 wzasymm2_scission = (std::sqrt (max( (70.0-28.0)/3.0*(z*z/a-35.0)+28.,0.0 )) ) * z/a;
676 } else {
677 wzasymm1_scission = wzasymm1_saddle;
678 wzasymm2_scission = wzasymm2_saddle;
679 }
680 }
681
682 wzasymm1_scission = max(wzasymm1_scission,wzasymm1_saddle);
683 wzasymm2_scission = max(wzasymm2_scission,wzasymm2_saddle);
684
685 wzasymm1 = wzasymm1_scission * fwidth_asymm1;
686 wzasymm2 = wzasymm2_scission * fwidth_asymm2;
687 wzsymm = wzsymm_scission;
688
689// /* if (ITEST == 1) {
690// G4cout << "WZ1(scission) = " << WZasymm1_scission << G4endl;
691// G4cout << "WZ2(scission) = " << WZasymm2_scission << G4endl;
692// G4cout << "WZsymm(scission) = " << WZsymm_scission << G4endl;
693// }
694// if (ITEST == 1) {
695// G4cout << "WZ1(scission) final= " << WZasymm1 << G4endl;
696// G4cout << "WZ2(scission) final= " << WZasymm2 << G4endl;
697// G4cout << "WZsymm(scission) final= " << WZsymm << G4endl;
698// } */
699
700 wasymm = wzsymm * a/z;
701 waheavy1 = wzasymm1 * a/z;
702 waheavy2 = wzasymm2 * a/z;
703
704 // G4cout <<"al, e, es, cn " << a_levdens << e << e_saddle_scission << cn << G4endl;
705
706 wasymm_saddle = wzsymm_saddle * a/z;
707 waheavy1_saddle = wzasymm1_saddle * a/z;
708 waheavy2_saddle = wzasymm2_saddle * a/z;
709
710 // if (itest == 1) {
711 // G4cout << "wasymm = " << wzsymm << G4endl;
712 // G4cout << "waheavy1 = " << waheavy1 << G4endl;
713 // G4cout << "waheavy2 = " << waheavy2 << G4endl;
714 // }
715
716 // sig_0 = quantum fluctuation = 0.45 z units for A=cte
717 // 0.45*2.58 = 1.16 n units for Z=cte
718 // sig_0^2 = 1.16*2 = 1.35 n units for Z=cte
719 n1width = std::sqrt(0.5 * std::sqrt(1.0/a_levdens*(e + e_saddle_scission)) / cn + 1.35);
720 if ( (epsilon0_1_saddle - delta_u1*std::exp((-epsilon_1_saddle*gamma_heavy1))) < 0.0) {
721 sasymm1 = -10.0;
722 } else {
723 sasymm1 = 2.0 * std::sqrt( a_levdens * (epsilon0_1_saddle -
724 delta_u1*(std::exp((-epsilon_1_saddle*gamma_heavy1))) ) );
725 }
726
727 if ( (epsilon0_2_saddle - delta_u2*std::exp((-epsilon_2_saddle*gamma_heavy2))) < 0.0) {
728 sasymm2 = -10.0;
729 } else {
730 sasymm2 = 2.0 * std::sqrt( a_levdens * (epsilon0_2_saddle -
731 delta_u2*(std::exp((-epsilon_2_saddle*gamma_heavy2))) ) );
732 }
733
734 if (epsilon_symm_saddle > 0.0) {
735 ssymm = 2.0 * std::sqrt( a_levdens*(epsilon_symm_saddle) );
736 } else {
737 ssymm = -10.0;
738 }
739
740 if (ssymm > -10.0) {
741 ysymm = 1.0;
742 if (epsilon0_1_saddle < 0.0) { // /* low energy */
743 yasymm1 = std::exp((sasymm1-ssymm)) * wzasymm1_saddle / wzsymm_saddle * 2.0;
744 // /* factor of 2 for symmetry classes */
745 } else { // /* high energy */
746 ssymm_mode1 = 2.0 * std::sqrt( a_levdens*(epsilon0_1_saddle) );
747 yasymm1 = ( std::exp((sasymm1-ssymm)) - std::exp((ssymm_mode1 - ssymm)) )
748 * wzasymm1_saddle / wzsymm_saddle * 2.0;
749 }
750
751 if (epsilon0_2_saddle < 0.0) { // /* low energy */
752 yasymm2 = std::exp((sasymm2-ssymm)) * wzasymm2_saddle / wzsymm_saddle * 2.0;
753 // /* factor of 2 for symmetry classes */
754 } else { // /* high energy */
755 ssymm_mode2 = 2.0 * std::sqrt( a_levdens*(epsilon0_2_saddle) );
756 yasymm2 = ( std::exp((sasymm2-ssymm)) - std::exp((ssymm_mode2 - ssymm)) )
757 * wzasymm2_saddle / wzsymm_saddle * 2.0;
758 }
759 // /* difference in the exponent in order */
760 // /* to avoid numerical overflow */
761 }
762 else {
763 if ( (sasymm1 > -10.0) && (sasymm2 > -10.0) ) {
764 ysymm = 0.0;
765 yasymm1 = std::exp(sasymm1) * wzasymm1_saddle * 2.0;
766 yasymm2 = std::exp(sasymm2) * wzasymm2_saddle * 2.0;
767 }
768 }
769
770 // /* normalize */
771 ysum = ysymm + yasymm1 + yasymm2;
772 if (ysum > 0.0) {
773 ysymm = ysymm / ysum;
774 yasymm1 = yasymm1 / ysum;
775 yasymm2 = yasymm2 / ysum;
776 yasymm = yasymm1 + yasymm2;
777 } else {
778 ysymm = 0.0;
779 yasymm1 = 0.0;
780 yasymm2 = 0.0;
781 // /* search minimum threshold and attribute all events to this mode */
782 if ( (epsilon_symm_saddle < epsilon_1_saddle) && (epsilon_symm_saddle < epsilon_2_saddle) ) {
783 ysymm = 1.0;
784 } else {
785 if (epsilon_1_saddle < epsilon_2_saddle) {
786 yasymm1 = 1.0;
787 } else {
788 yasymm2 = 1.0;
789 }
790 }
791 }
792
793// if (itest == 1) {
794// G4cout << "ysymm normalized= " << ysymm << G4endl;
795// G4cout << "yasymm1 normalized= " << yasymm1 << G4endl;
796// G4cout << "yasymm2 normalized= " << yasymm2 << G4endl;
797// }
798
799 // /* even-odd effect */
800 // /* simple parametrization KHS, Nov. 2000. From Rejmund et al. */
801 eexc_max = max(eexc1_saddle, eexc2_saddle);
802 eexc_max = max(eexc_max, e);
803 iz = (G4int)z;
804 // G4cout << "mod(z, 2)" << iz%2 << G4endl;
805 if ((G4int)(z) % 2 == 0) {
806 r_e_o_max = 0.3 * (1.0 - 0.2 * (std::pow(z, 2)/a - std::pow(92.0, 2)/238.0));
807 e_pair = 2.0 * 12.0 / std::sqrt(a);
808 if(eexc_max > (e_crit + e_pair)) {
809 r_e_o = 0.0;
810 } else {
811 if(eexc_max < e_pair) {
812 r_e_o = r_e_o_max;
813 } else {
814 r_e_o = std::pow((eexc_max - e_crit - e_pair)/e_crit, 2) * r_e_o_max;
815 }
816 }
817 } else {
818 r_e_o = 0.0;
819 }
820
821 // G4cout <<"rmax " << r_e_o_max << G4endl;
822 // if(r_e_o > 0.0) G4cout <<"e_crit, r_e_o" << e_crit << r_e_o << G4endl;
823 // $LOOP; /* event loop */
824 // I_COUNT = I_COUNT + 1;
825
826 /* random decision: symmetric or asymmetric */
827 /* IMODE = 3 means asymmetric fission, mode 1,
828 IMODE = 2 means asymmetric fission, mode 2,
829 IMODE = 1 means symmetric */
830 // RMODE = dble(HAZ(k));
831 // rmode = rnd.rndm();
832// // Safety check added to make sure we always select well defined
833// // fission mode.
834 rmode = haz(k);
835 // Cast for test CS 11/10/05
836 // RMODE = 0.54;
837 // rmode = 0.54;
838 if (rmode < ysymm) {
839 imode = 1;
840 } else if (rmode < (ysymm + yasymm1)) {
841 imode = 2;
842 } else {
843 imode = 3;
844 }
845 // /* determine parameters of the Z distribution */
846 // force imode (for testing, PK)
847 // imode = 3;
848
849 if (imode == 1) {
850 z1mean = zsymm;
851 z1width = wzsymm;
852 } else if (imode == 2) {
853 z1mean = zheavy1;
854 z1width = wzasymm1;
855 } else if (imode == 3) {
856 z1mean = zheavy2;
857 z1width = wzasymm2;
858 }
859
860 if (itest == 1) {
861 G4cout << "nbre aleatoire tire " << rmode << G4endl;
862 G4cout << "fission mode " << imode << G4endl;
863 G4cout << "z1mean= " << z1mean << G4endl;
864 G4cout << "z1width= " << z1width << G4endl;
865 }
866
867 // /* random decision: Z1 and Z2 at scission: */
868 z1 = 1.0;
869 z2 = 1.0;
870
871 while ( (z1<5.0) || (z2<5.0) ) {
872 // Z1 = dble(GAUSSHAZ(K,sngl(Z1mean),sngl(Z1width)));
873 // z1 = rnd.gaus(z1mean,z1width);
874 // z1 = 48.26; // gausshaz(k, z1mean, z1width);
875 z1 = gausshaz(k, z1mean, z1width);
876 even_odd(z1, r_e_o, i_help);
877 z1 = double(i_help);
878 z2 = z - z1;
879 }
880
881 if (itest == 1) {
882 G4cout << "ff charge sample " << G4endl;
883 G4cout << "z1 = " << z1 << G4endl;
884 G4cout << "z2 = " << z2 << G4endl;
885 }
886
887// // CALL EVEN_ODD(Z1,R_E_O,I_HELP);
888// // /* Integer proton number with even-odd effect */
889// // Z1 = REAL(I_HELP)
890// // /* Z1 = INT(Z1+0.5E0); */
891// z2 = z - z1;
892
893 // /* average N of both fragments: */
894 if (imode == 1) {
895 n1ucd = z1 * n/z;
896 n2ucd = z2 * n/z;
897 re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) + ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,2.0); // umass == massdef
898 re2 = umass(z1,n1ucd+1.,0.6) + umass(z2,n2ucd-1.,0.6) + ecoul(z1,n1ucd+1.,0.6,z2,n2ucd-1.,0.6,2.0);
899 re3 = umass(z1,n1ucd+2.,0.6) + umass(z2,n2ucd-2.,0.6) + ecoul(z1,n1ucd+2.,0.6,z2,n2ucd-2.,0.6,2.0);
900 reps2 = (re1-2.0*re2+re3) / 2.0;
901 reps1 = re2 - re1 - reps2;
902 rn1_pol = - reps1 / (2.0 * reps2);
903 n1mean = n1ucd + rn1_pol;
904 } else {
905 n1mean = (z1 + 0.5) * n/z;
906 }
907
908// n1mean nsymm + (z1 - zsymm) * 1.6 from 238 U(nth, f)
909// n1width = 0.9 + E * 0.002 KHS
910
911// random decision: N1R and N2R at scission, before evaporation
912 n1r = 1.0;
913 n2r = 1.0;
914 while (n1r < 5 || n2r < 5) {
915 // n1r = 76.93; gausshaz(kkk,n1mean,n1width);
916 n1r = gausshaz(kkk,n1mean,n1width);
917 // modification to have n1r as integer, and n=n1r+n2r rigorously a.b. 19/4/2001
918 i_inter = int(n1r + 0.5);
919 n1r = double(i_inter);
920 n2r = n - n1r;
921 }
922
923 // neutron evaporation from fragments
924 if (i_eva > 0) {
925 // treatment sz
926 ne_min = 0.095e0 * a - 20.4e0;
927 if (ne_min < 0) ne_min = 0.0;
928 ne_min = ne_min + e / 8.e0; // 1 neutron per 8 mev */
929 a1r = z1 + n1r; // mass of first fragment */
930 ne_m1 = a1r / a * ne_min; // devide nbr. of neutrons acc. mass
931 ne_m2 = ne_min - ne_m1; // nmbr. of neutrons of 2. fragment
932 n1 = n1r - ne_m1; // final neutron number 1. fragment
933 n2 = n2r - ne_m2; // ! final neutron number 2. fragment
934 } else {
935 n1 = n1r;
936 n2 = n2r;
937 }
938
939 // excitation energy due to deformation
940
941 a1 = z1 + n1r; // mass of first fragment */
942 a2 = z2 + n2r; // mass of second fragment */
943 if (a1 < 80) {
944 ed1_low = 0.0;
945 } else if (a1 >= 80 && a1 < 110) {
946 ed1_low = (a1-80.)*20./30.;
947 } else if (a1 >= 110 && a1 < 130) {
948 ed1_low = -(a1-110.)*20./20. + 20.;
949 } else if (a1 >= 130) {
950 ed1_low = (a1-130.)*20./30.;
951 }
952
953 if (a2 < 80) {
954 ed2_low = 0.0;
955 } else if (a2 >= 80 && a2 < 110) {
956 ed2_low = (a2-80.)*20./30.;
957 } else if (a2 >= 110 && a2 < 130) {
958 ed2_low = -(a2-110.)*20./20. + 20.;
959 } else if (a2 >= 130) {
960 ed2_low = (a2-130.)*20./30.;
961 }
962
963 ed1_high = 20.0*a1/(a1+a2);
964 ed2_high = 20.0 - ed1_high;
965 ed1 = ed1_low*std::exp(-e/el) + ed1_high*(1-std::exp(-e/el));
966 ed2 = ed2_low*std::exp(-e/el) + ed2_high*(1-std::exp(-e/el));
967
968 // write(6,101)e,a1,a2,ed1,ed2,ed1+ed2
969 // write(6,102)ed1_low,ed1_high,ed2_low,ed2_high
970 e1 = e*a1/(a1+a2) + ed1;
971 e2 = e - e*a1/(a1+a2) + ed2;
972 atot = a1+a2;
973 if (atot > a+1) {
974 // write(6,*)'a,,a1,a2,atot',a,a1,a2,atot
975 // write(6,*)'n,n1r,n2r',n,n1r,n2r
976 // write(6,*)'z,z1,z2',z,z1,z2
977 }
978
979 milledeux:
980 // only symmetric fission
981 // Symmetric fission: Ok! Checked CS 10/10/05
982 if ( (icz == -1) || (a1 < 0.0) || (a2 < 0.0) ) {
983 // IF (z.eq.92) THEN
984 // write(6,*)'symmetric fission'
985 // write(6,*)'Z,A,E,A1,A2,icz,Atot',Z,A,E,A1,A2,icz,Atot
986 // END IF
987
988 if (itest == 1) {
989 G4cout << "milledeux: liquid-drop option " << G4endl;
990 }
991
992 n = a-z;
993 // proton number in symmetric fission (centre) *
994 zsymm = z / 2.0;
995 nsymm = n / 2.0;
996 asymm = nsymm + zsymm;
997
998 a_levdens = a / 8.0;
999
1000 masscurv = 2.0;
1001 cz_symm = 8.0 / std::pow(z,2) * masscurv;
1002
1003 wzsymm = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cz_symm) ) ;
1004
1005 if (itest == 1) {
1006 G4cout << " symmetric high energy fission " << G4endl;
1007 G4cout << "wzsymm " << wzsymm << G4endl;
1008 }
1009
1010 z1mean = zsymm;
1011 z1width = wzsymm;
1012
1013 // random decision: Z1 and Z2 at scission: */
1014 z1 = 1.0;
1015 z2 = 1.0;
1016 while ( (z1 < 5.0) || (z2 < 5.0) ) {
1017 // z1 = dble(gausshaz(kkk,sngl(z1mean),sngl(z1width)));
1018 // z1 = rnd.gaus(z1mean,z1width);
1019 // z1 = 24.8205585; //gausshaz(kkk, z1mean, z1width);
1020 z1 = gausshaz(kkk, z1mean, z1width);
1021 z2 = z - z1;
1022 }
1023
1024 if (itest == 1) {
1025 G4cout << " z1 " << z1 << G4endl;
1026 G4cout << " z2 " << z2 << G4endl;
1027 }
1028 if (itest == 1) {
1029 G4cout << " zsymm " << zsymm << G4endl;
1030 G4cout << " nsymm " << nsymm << G4endl;
1031 G4cout << " asymm " << asymm << G4endl;
1032 }
1033
1034 cn = umass(zsymm, nsymm+1.0, 0.0) + umass(zsymm, nsymm-1.0, 0.0)
1035 + 1.44 * std::pow(zsymm, 2)/
1036 (std::pow(r_null, 2) * std::pow(std::pow(asymm+1.0, 1.0/3.0) + std::pow(asymm-1.0, 1.0/3.0), 2))
1037 - 2.0 * umass(zsymm, nsymm, 0.0) - 1.44 * std::pow(zsymm, 2) /
1038 std::pow(r_null * 2.0 *std::pow(asymm, 1.0/3.0), 2);
1039 // This is an approximation! Coulomb energy is neglected.
1040
1041 n1width = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cn) + 1.35);
1042 if (itest == 1) {
1043 G4cout << " cn " << cn << G4endl;
1044 G4cout << " n1width " << n1width << G4endl;
1045 }
1046
1047 // /* average N of both fragments: */
1048 n1ucd = z1 * n/z;
1049 n2ucd = z2 * n/z;
1050 re1 = umass(z1,n1ucd, 0.6) + umass(z2,n2ucd, 0.6) + ecoul(z1,n1ucd, 0.6,z2,n2ucd, 0.6,2.0);
1051 re2 = umass(z1,n1ucd+1.,0.6) + umass(z2,n2ucd-1.,0.6) + ecoul(z1,n1ucd+1.,0.6,z2,n2ucd-1.,0.6,2.0);
1052 re3 = umass(z1,n1ucd+2.,0.6) + umass(z2,n2ucd-2.,0.6) + ecoul(z1,n1ucd+2.,0.6,z2,n2ucd-2.,0.6,2.0);
1053 reps2 = (re1-2.0*re2+re3) / 2.0;
1054 reps1 = re2 - re1 - reps2;
1055 rn1_pol = - reps1 / (2.0 * reps2);
1056 n1mean = n1ucd + rn1_pol;
1057
1058 // random decision: N1R and N2R at scission, before evaporation: */
1059 // N1R = dfloat(NINT(GAUSSHAZ(KKK,sngl(N1mean),sngl(N1width))));
1060 // n1r = (float)( (int)(rnd.gaus(n1mean,n1width)) );
1061 // n1r = 34.0; //(float)( (int)(gausshaz(k, n1mean,n1width)) );
1062 n1r = (float)( (int)(gausshaz(k, n1mean,n1width)) );
1063 n2r = n - n1r;
1064 // Mass of first and second fragment */
1065 a1 = z1 + n1r;
1066 a2 = z2 + n2r;
1067
1068 e1 = e*a1/(a1+a2);
1069 e2 = e - e*a1/(a1+a2);
1070 }
1071 v1 = 0.0; // These are not calculated in SimFis3.
1072 v2 = 0.0;
1073 if (itest == 1) {
1074 G4cout << " n1r " << n1r << G4endl;
1075 G4cout << " n2r " << n2r << G4endl;
1076 }
1077
1078 if (itest == 1) {
1079 G4cout << " a1 " << a1 << G4endl;
1080 G4cout << " z1 " << z1 << G4endl;
1081 G4cout << " a2 " << a2 << G4endl;
1082 G4cout << " z2 " << z2 << G4endl;
1083 G4cout << " e1 " << e1 << G4endl;
1084 G4cout << " e2 " << e << G4endl;
1085 }
1086}
1087
1088void G4AblaFission::standardRandom(G4double *rndm, G4long *seed)
1089{
1090 (*seed) = (*seed); // Avoid warning during compilation.
1091 // Use Geant4 G4UniformRand
1092 (*rndm) = randomGenerator->getRandom();
1093}
1094
1095G4double G4AblaFission::haz(G4int k)
1096{
1097 const G4int pSize = 110;
1098 static G4double p[pSize];
1099 static G4long ix = 0, i = 0;
1100 static G4double x = 0.0, y = 0.0, a = 0.0, haz = 0.0;
1101 // k =< -1 on initialise
1102 // k = -1 c'est reproductible
1103 // k < -1 || k > -1 ce n'est pas reproductible
1104
1105 // Zero is invalid random seed. Set proper value from our random seed collection:
1106 if(ix == 0) {
1107 ix = hazard->ial;
1108 }
1109
1110 if (k <= -1) { //then
1111 if(k == -1) { //then
1112 ix = 0;
1113 }
1114 else {
1115 x = 0.0;
1116 y = secnds(int(x));
1117 ix = int(y * 100 + 43543000);
1118 if(mod(ix,2) == 0) {
1119 ix = ix + 1;
1120 }
1121 }
1122
1123 // Here we are using random number generator copied from INCL code
1124 // instead of the CERNLIB one! This causes difficulties for
1125 // automatic testing since the random number generators, and thus
1126 // the behavior of the routines in C++ and FORTRAN versions is no
1127 // longer exactly the same!
1128 x = randomGenerator->getRandom();
1129 // standardRandom(&x, &ix);
1130 for(G4int i = 0; i < pSize; i++) { //do i=1,110
1131 p[i] = randomGenerator->getRandom();
1132 // standardRandom(&(p[i]), &ix);
1133 }
1134 a = randomGenerator->getRandom();
1135 //standardRandom(&a, &ix);
1136 k = 0;
1137 }
1138
1139 i = nint(100*a)+1;
1140 haz = p[i];
1141 a = randomGenerator->getRandom();
1142 // standardRandom(&a, &ix);
1143 p[i] = a;
1144
1145 hazard->ial = ix;
1146 // haz=0.4;
1147 return haz;
1148}
1149
1150G4double G4AblaFission::gausshaz(int k, double xmoy, double sig)
1151{
1152 // Gaussian random numbers:
1153
1154 // 1005 C*** TIRAGE ALEATOIRE DANS UNE GAUSSIENNE DE LARGEUR SIG ET MOYENNE XMOY
1155 static G4int iset = 0;
1156 static G4double v1,v2,r,fac,gset,gausshaz;
1157
1158 if(iset == 0) { //then
1159 do {
1160 v1 = 2.0*haz(k) - 1.0;
1161 v2 = 2.0*haz(k) - 1.0;
1162 r = std::pow(v1,2) + std::pow(v2,2);
1163 } while(r >= 1);
1164
1165 fac = std::sqrt(-2.*std::log(r)/r);
1166 gset = v1*fac;
1167 gausshaz = v2*fac*sig+xmoy;
1168 iset = 1;
1169 }
1170 else {
1171 gausshaz=gset*sig+xmoy;
1172 iset=0;
1173 }
1174 return gausshaz;
1175}
1176
1177// Utilities
1178
1179G4double G4AblaFission::min(G4double a, G4double b)
1180{
1181 if(a < b) {
1182 return a;
1183 }
1184 else {
1185 return b;
1186 }
1187}
1188
1189G4int G4AblaFission::min(G4int a, G4int b)
1190{
1191 if(a < b) {
1192 return a;
1193 }
1194 else {
1195 return b;
1196 }
1197}
1198
1199G4double G4AblaFission::max(G4double a, G4double b)
1200{
1201 if(a > b) {
1202 return a;
1203 }
1204 else {
1205 return b;
1206 }
1207}
1208
1209G4int G4AblaFission::max(G4int a, G4int b)
1210{
1211 if(a > b) {
1212 return a;
1213 }
1214 else {
1215 return b;
1216 }
1217}
1218
1219G4int G4AblaFission::nint(G4double number)
1220{
1221 G4double intpart = 0.0;
1222 G4double fractpart = 0.0;
1223 fractpart = std::modf(number, &intpart);
1224 if(number == 0) {
1225 return 0;
1226 }
1227 if(number > 0) {
1228 if(fractpart < 0.5) {
1229 return int(std::floor(number));
1230 }
1231 else {
1232 return int(std::ceil(number));
1233 }
1234 }
1235 if(number < 0) {
1236 if(fractpart < -0.5) {
1237 return int(std::floor(number));
1238 }
1239 else {
1240 return int(std::ceil(number));
1241 }
1242 }
1243
1244 return int(std::floor(number));
1245}
1246
1247G4int G4AblaFission::secnds(G4int x)
1248{
1249 time_t mytime;
1250 tm *mylocaltime;
1251
1252 time(&mytime);
1253 mylocaltime = localtime(&mytime);
1254
1255 if(x == 0) {
1256 return(mylocaltime->tm_hour*60*60 + mylocaltime->tm_min*60 + mylocaltime->tm_sec);
1257 }
1258 else {
1259 return(mytime - x);
1260 }
1261}
1262
1263G4int G4AblaFission::mod(G4int a, G4int b)
1264{
1265 if(b != 0) {
1266 return (a - (a/b)*b);
1267 }
1268 else {
1269 return 0;
1270 }
1271}
1272
1273G4double G4AblaFission::dmod(G4double a, G4double b)
1274{
1275 if(b != 0) {
1276 return (a - (a/b)*b);
1277 }
1278 else {
1279 return 0.0;
1280 }
1281}
1282
1283G4double G4AblaFission::dint(G4double a)
1284{
1285 G4double value = 0.0;
1286
1287 if(a < 0.0) {
1288 value = double(std::ceil(a));
1289 }
1290 else {
1291 value = double(std::floor(a));
1292 }
1293
1294 return value;
1295}
1296
1297G4int G4AblaFission::idint(G4double a)
1298{
1299 G4int value = 0;
1300
1301 if(a < 0) {
1302 value = int(std::ceil(a));
1303 }
1304 else {
1305 value = int(std::floor(a));
1306 }
1307
1308 return value;
1309}
1310
1311G4int G4AblaFission::idnint(G4double value)
1312{
1313 G4double valueCeil = int(std::ceil(value));
1314 G4double valueFloor = int(std::floor(value));
1315
1316 if(std::fabs(value - valueCeil) < std::fabs(value - valueFloor)) {
1317 return int(valueCeil);
1318 }
1319 else {
1320 return int(valueFloor);
1321 }
1322}
1323
1324G4double G4AblaFission::dmin1(G4double a, G4double b, G4double c)
1325{
1326 if(a < b && a < c) {
1327 return a;
1328 }
1329 if(b < a && b < c) {
1330 return b;
1331 }
1332 if(c < a && c < b) {
1333 return c;
1334 }
1335 return a;
1336}
1337
1338G4double G4AblaFission::utilabs(G4double a)
1339{
1340 if(a > 0) {
1341 return a;
1342 }
1343 if(a < 0) {
1344 return (-1*a);
1345 }
1346 if(a == 0) {
1347 return a;
1348 }
1349
1350 return a;
1351}
1352
1353void G4AblaFission::p(G4int a, G4int b) {
1354 G4cout << a << std::setw(8) << b << G4endl;
1355}
1356
1357void G4AblaFission::p(G4double a, G4double b) {
1358 G4cout << " " << a << std::setw(12) << b << G4endl;
1359}
1360
1361void G4AblaFission::p(G4String msg, G4double a, G4double b) {
1362 G4cout << " " << msg << " " << a << std::setw(12) << b << G4endl;
1363}
Note: See TracBrowser for help on using the repository browser.