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