Geant4  10.00.p03
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 
387  // /* average Z of asymmetric and symmetric components: */
388  n = a - z; /* neutron number of the fissioning nucleus */
389 
390  k = 0;
391  icz = 0;
392  if ( (std::pow(z,2)/a < 25.0) || (n < nheavy2) || (e > 500.0) ) {
393  icz = -1;
394  // GOTO 1002;
395  goto milledeux;
396  }
397 
398  // /* Polarisation assumed for standard I and standard II:
399  // Z - Zucd = cpol (for A = const);
400  // from this we get (see Armbruster)
401  // Z - Zucd = Acn/Ncn * cpol (for N = const) */
402 
403  // zheavy1_shell = ((nheavy1/n) * z) - ((a/n) * cpol1); // Simfis18 PK:::
404  zheavy1_shell = ((nheavy1/n) * z) - 0.8; // Simfis3 PK:::
405  //zheavy2_shell = ((nheavy2/n) * z) - ((a/n) * cpol2); // Simfis18 PK:::
406  zheavy2_shell = ((nheavy2/n) * z) - 0.8; // Simfis3 PK:::
407 
408  // p(zheavy1_shell, zheavy2_shell);
409 
410  // e_saddle_scission =
411  // (-24.0 + 0.02227 * (std::pow(z,2))/(std::pow(a,0.33333)) ) * friction_factor;
412  e_saddle_scission = (3.535 * std::pow(z,2)/a - 121.1) * friction_factor; // PK:::
413 
414  // /* Energy dissipated from saddle to scission */
415  // /* F. Rejmund et al., Nucl. Phys. A 678 (2000) 215, fig. 4 b */
416  // E_saddle_scission = DMAX1(0.,E_saddle_scission);
417  // Heavy Ion Induced Reactions, Schroeder W. ed., Harwood, 1986, 101
418  if (e_saddle_scission < 0.0) {
419  e_saddle_scission = 0.0;
420  }
421 
422  // /* Semiempirical fission model: */
423 
424  // /* Fit to experimental result on curvature of potential at saddle */
425  // /* reference: */
426  // /* IF Z**2/A < 33.15E0 THEN
427  // MassCurv = 30.5438538E0 - 4.00212049E0 * Z**2/A
428  // + 0.11983384E0 * Z**4 / (A**2) ;
429  // ELSE
430  // MassCurv = 10.E0 ** (7.16993332E0 - 0.26602401E0 * Z**2/A
431  // + 0.00283802E0 * Z**4 / (A**2)) ; */
432  // /* New parametrization of T. Enqvist according to Mulgin et al. 1998 NPA 640(1998) 375 */
433  if ( (std::pow(z,2))/a < 33.9186) {
434  masscurv = std::pow( 10.0,(-1.093364 + 0.082933 * (std::pow(z,2)/a)
435  - 0.0002602 * (std::pow(z,4)/std::pow(a,2))) );
436  } else {
437  masscurv = std::pow( 10.0,(3.053536 - 0.056477 * (std::pow(z,2)/a)
438  + 0.0002454 * (std::pow(z,4)/std::pow(a,2))) );
439  }
440 
441  cz_symm = (8.0/std::pow(z,2)) * masscurv;
442 
443  if(itest == 1) {
444  // // G4cout << "cz_symmetry= " << cz_symm << G4endl;
445  }
446 
447  icz = 0;
448  if (cz_symm < 0) {
449  icz = -1;
450  // GOTO 1002;
451  goto milledeux;
452  }
453 
454  // /* proton number in symmetric fission (centre) */
455  zsymm = z/2.0;
456  nsymm = n/2.0;
457  asymm = nsymm + zsymm;
458 
459  zheavy1 = (cz_symm*zsymm + cz_asymm1_shell*zheavy1_shell)/(cz_symm + cz_asymm1_shell);
460  zheavy2 = (cz_symm*zsymm + cz_asymm2_shell*zheavy2_shell)/(cz_symm + cz_asymm2_shell);
461 
462  // /* position of valley due to influence of liquid-drop potential */
463  nheavy1_eff = (zheavy1 + 0.8)*(n/z);
464  nheavy2_eff = (zheavy2 + 0.8)*(n/z);
465  aheavy1 = nheavy1_eff + zheavy1;
466  aheavy2 = nheavy2_eff + zheavy2;
467  // Eheavy1 = E * Aheavy1 / A
468  // Eheavy2 = E * Aheavy2 / A
469  // Elight1 = E * Alight1 / A
470  // Elight2 = E * Alight2 / A
471  a_levdens = a / 8.0;
472  a_levdens_heavy1 = aheavy1 / 8.0;
473  a_levdens_heavy2 = aheavy2 / 8.0;
474  gamma = a_levdens / (0.4 * (std::pow(a,1.3333)) );
475  gamma_heavy1 = ( a_levdens_heavy1 / (0.4 * (std::pow(aheavy1,1.3333)) ) ) * fgamma1;
476  gamma_heavy2 = a_levdens_heavy2 / (0.4 * (std::pow(aheavy2,1.3333)) );
477 
478  // Up to here: Ok! Checked CS 10/10/05
479 
480  cn = umass(zsymm,(nsymm+1.),0.0) + umass(zsymm,(nsymm-1.),0.0)
481  + 1.44 * (std::pow(zsymm,2))/
482  ( (std::pow(r_null,2)) *
483  ( std::pow((asymm+1.0),1.0/3.0) + std::pow((asymm-1.0),1.0/3.0) ) *
484  ( std::pow((asymm+1.0),1.0/3.0) + std::pow((asymm-1.0),1.0/3.0) ) )
485  - 2.0 * umass(zsymm,nsymm,0.0)
486  - 1.44 * (std::pow(zsymm,2))/
487  ( ( 2.0 * r_null * (std::pow(asymm,1.0/3.0)) ) *
488  ( 2.0 * r_null * (std::pow(asymm,1.0/3.0)) ) );
489 
490  // /* shell effect in valley of mode 1 */
491  delta_u1 = delta_u1_shell + (std::pow((zheavy1_shell-zheavy1),2))*cz_asymm1_shell;
492  // /* shell effect in valley of mode 2 */
493  delta_u2 = delta_u2_shell + (std::pow((zheavy2_shell-zheavy2),2))*cz_asymm2_shell;
494 
495  Bsym = 0.0;
496  Basym_1 = Bsym + std::pow((zheavy1-zsymm), 2) * cz_symm + delta_u1;
497  Basym_2 = Bsym + std::pow((zheavy2-zsymm), 2) * cz_symm + delta_u2;
498  if(Bsym < Basym_1 && Bsym < Basym_2) {
499  // Excitation energies at the saddle point
500  // without and with shell effect
501  epsilon0_1_saddle = (e + e_zero_point - std::pow((zheavy1 - zsymm), 2) * cz_symm);
502  epsilon0_2_saddle = (e + e_zero_point - std::pow((zheavy2 - zsymm), 2) * cz_symm);
503 
504  epsilon_1_saddle = epsilon0_1_saddle - delta_u1;
505  epsilon_2_saddle = epsilon0_2_saddle - delta_u2;
506 
507  epsilon_symm_saddle = e + e_zero_point;
508  eexc1_saddle = epsilon_1_saddle;
509  eexc2_saddle = epsilon_2_saddle;
510 
511  // Excitation energies at the scission point
512  // heavy fragment without and with shell effect
513  epsilon0_1_scission = (e + e_saddle_scission - std::pow((zheavy1 - zsymm), 2) * cz_symm) * aheavy1/a;
514  epsilon_1_scission = epsilon0_1_scission - delta_u1*(aheavy1/a);
515 
516  epsilon0_2_scission = (e + e_saddle_scission - std::pow((zheavy2 - zsymm), 2) * cz_symm) * aheavy2/a;
517  epsilon_2_scission = epsilon0_2_scission - delta_u2*(aheavy2/a);
518 
519  epsilon_symm_scission = e + e_saddle_scission;
520  } else if (Basym_1 < Bsym && Basym_1 < Basym_2) {
521  // Excitation energies at the saddle point
522  // without and with shell effect
523  epsilon_symm_saddle = (e + e_zero_point + delta_u1 + std::pow((zheavy1-zsymm), 2) * cz_symm);
524  epsilon0_2_saddle = (epsilon_symm_saddle - std::pow((zheavy2-zsymm), 2) * cz_symm);
525  epsilon_2_saddle = epsilon0_2_saddle - delta_u2;
526  epsilon0_1_saddle = e + e_zero_point + delta_u1;
527  epsilon_1_saddle = e + e_zero_point;
528  eexc1_saddle = epsilon_1_saddle;
529  eexc2_saddle = epsilon_2_saddle;
530 
531  // Excitation energies at the scission point
532  // heavy fragment without and with shell effect
533  epsilon_symm_scission = (e + e_saddle_scission + std::pow((zheavy1-zsymm), 2) * cz_symm + delta_u1);
534  epsilon0_2_scission = (epsilon_symm_scission - std::pow((zheavy2-zsymm), 2) * cz_symm) * aheavy2/a;
535  epsilon_2_scission = epsilon0_2_scission - delta_u2*aheavy2/a;
536  epsilon0_1_scission = (e + e_saddle_scission + delta_u1) * aheavy1/a;
537  epsilon_1_scission = (e + e_saddle_scission) * aheavy1/a;
538  } else if (Basym_2 < Bsym && Basym_2 < Basym_1) {
539  // Excitation energies at the saddle point
540  // without and with shell effect
541  epsilon_symm_saddle = (e + e_zero_point + delta_u2 + std::pow((zheavy2-zsymm), 2) * cz_symm);
542  epsilon0_1_saddle = (epsilon_symm_saddle - std::pow((zheavy1-zsymm), 2) * cz_symm);
543  epsilon_1_saddle = epsilon0_1_saddle - delta_u1;
544  epsilon0_2_saddle = e + e_zero_point + delta_u2;
545  epsilon_2_saddle = e + e_zero_point;
546  eexc1_saddle = epsilon_1_saddle;
547  eexc2_saddle = epsilon_2_saddle;
548 
549  // Excitation energies at the scission point
550  // heavy fragment without and with shell effect
551  epsilon_symm_scission = (e + e_saddle_scission + std::pow((zheavy2-zsymm), 2) * cz_symm + delta_u2);
552  epsilon0_1_scission = (epsilon_symm_scission - std::pow((zheavy1-zsymm), 2) * cz_symm) * aheavy1/a;
553  epsilon_1_scission = epsilon0_1_scission - delta_u1*aheavy1/a;
554  epsilon0_2_scission = (e + e_saddle_scission + delta_u2) * aheavy2/a;
555  epsilon_2_scission = (e + e_saddle_scission) * aheavy2/a;
556 
557  } else {
558  // G4cout <<"G4AblaFission: " << G4endl;
559  }
560  if(epsilon_1_saddle < 0.0) epsilon_1_saddle = 0.0;
561  if(epsilon_2_saddle < 0.0) epsilon_2_saddle = 0.0;
562  if(epsilon0_1_saddle < 0.0) epsilon0_1_saddle = 0.0;
563  if(epsilon0_2_saddle < 0.0) epsilon0_2_saddle = 0.0;
564  if(epsilon_symm_saddle < 0.0) epsilon_symm_saddle = 0.0;
565 
566  if(epsilon_1_scission < 0.0) epsilon_1_scission = 0.0;
567  if(epsilon_2_scission < 0.0) epsilon_2_scission = 0.0;
568  if(epsilon0_1_scission < 0.0) epsilon0_1_scission = 0.0;
569  if(epsilon0_2_scission < 0.0) epsilon0_2_scission = 0.0;
570  if(epsilon_symm_scission < 0.0) epsilon_symm_scission = 0.0;
571 
572  if(itest == 1) {
573  // G4cout <<"E, E1, E2, Es" << e << epsilon_1_saddle << epsilon_2_saddle << epsilon_symm_saddle << G4endl;
574  }
575 
576  e_eff1_saddle = epsilon0_1_saddle - delta_u1 * (std::exp((-epsilon_1_saddle*gamma)));
577 
578  if (e_eff1_saddle > 0.0) {
579  wzasymm1_saddle = std::sqrt( (0.5) *
580  (std::sqrt(1.0/a_levdens*e_eff1_saddle)) /
581  (cz_asymm1_shell * std::exp((-epsilon_1_saddle*gamma)) + cz_symm) );
582  } else {
583  wzasymm1_saddle = 1.0;
584  }
585 
586  e_eff2_saddle = epsilon0_2_saddle - delta_u2 * std::exp((-epsilon_2_saddle*gamma));
587  if (e_eff2_saddle > 0.0) {
588  wzasymm2_saddle = std::sqrt( (0.5 *
589  (std::sqrt(1.0/a_levdens*e_eff2_saddle)) /
590  (cz_asymm2_shell * std::exp((-epsilon_2_saddle*gamma)) + cz_symm) ) );
591  } else {
592  wzasymm2_saddle = 1.0;
593  }
594 
595  if(e - e_zero_point > 0.0) {
596  wzsymm_saddle = std::sqrt( (0.5 *
597  (std::sqrt(1.0/a_levdens*(e+epsilon_symm_saddle))) / cz_symm ) );
598  } else {
599  wzsymm_saddle = 1.0;
600  }
601 
602 // if (itest == 1) {
603 // // G4cout << "wz1(saddle) = " << wzasymm1_saddle << G4endl;
604 // // G4cout << "wz2(saddle) = " << wzasymm2_saddle << G4endl;
605 // // G4cout << "wzsymm(saddle) = " << wzsymm_saddle << G4endl;
606 // }
607 
608  // /* Calculate widhts at the scission point: */
609  // /* fits of ref. Beizin 1991 (Plots brought to GSI by Sergei Zhdanov) */
610 
611  wzsymm_scission = wzsymm_saddle;
612 
613  if (e_saddle_scission == 0.0) {
614  wzasymm1_scission = wzasymm1_saddle;
615  wzasymm2_scission = wzasymm2_saddle;
616  } else {
617  if (nheavy1_eff > 75.0) {
618  wzasymm1_scission = (std::sqrt(21.0)) * z/a;
619  double RR = (70.0-28.0)/3.0*(z*z/a-35.0)+28.0;
620  if(RR > 0.0) {
621  wzasymm2_scission = std::sqrt(RR)*(z/a);
622  } else {
623  wzasymm2_scission = 0.0;
624  }
625  wzasymm2_scission = (std::sqrt (max( (70.0-28.0)/3.0*(z*z/a-35.0)+28.,0.0 )) ) * z/a;
626  } else {
627  wzasymm1_scission = wzasymm1_saddle;
628  wzasymm2_scission = wzasymm2_saddle;
629  }
630  }
631 
632  wzasymm1_scission = max(wzasymm1_scission,wzasymm1_saddle);
633  wzasymm2_scission = max(wzasymm2_scission,wzasymm2_saddle);
634 
635  wzasymm1 = wzasymm1_scission * fwidth_asymm1;
636  wzasymm2 = wzasymm2_scission * fwidth_asymm2;
637  wzsymm = wzsymm_scission;
638 
639 // /* if (ITEST == 1) {
640 // // G4cout << "WZ1(scission) = " << WZasymm1_scission << G4endl;
641 // // G4cout << "WZ2(scission) = " << WZasymm2_scission << G4endl;
642 // // G4cout << "WZsymm(scission) = " << WZsymm_scission << G4endl;
643 // }
644 // if (ITEST == 1) {
645 // // G4cout << "WZ1(scission) final= " << WZasymm1 << G4endl;
646 // // G4cout << "WZ2(scission) final= " << WZasymm2 << G4endl;
647 // // G4cout << "WZsymm(scission) final= " << WZsymm << G4endl;
648 // } */
649 
650 
651  // // G4cout <<"al, e, es, cn " << a_levdens << e << e_saddle_scission << cn << G4endl;
652 
653  // if (itest == 1) {
654  // // G4cout << "wasymm = " << wzsymm << G4endl;
655  // // G4cout << "waheavy1 = " << waheavy1 << G4endl;
656  // // G4cout << "waheavy2 = " << waheavy2 << G4endl;
657  // }
658 
659  // sig_0 = quantum fluctuation = 0.45 z units for A=cte
660  // 0.45*2.58 = 1.16 n units for Z=cte
661  // sig_0^2 = 1.16*2 = 1.35 n units for Z=cte
662  n1width = std::sqrt(0.5 * std::sqrt(1.0/a_levdens*(e + e_saddle_scission)) / cn + 1.35);
663  if ( (epsilon0_1_saddle - delta_u1*std::exp((-epsilon_1_saddle*gamma_heavy1))) < 0.0) {
664  sasymm1 = -10.0;
665  } else {
666  sasymm1 = 2.0 * std::sqrt( a_levdens * (epsilon0_1_saddle -
667  delta_u1*(std::exp((-epsilon_1_saddle*gamma_heavy1))) ) );
668  }
669 
670  if ( (epsilon0_2_saddle - delta_u2*std::exp((-epsilon_2_saddle*gamma_heavy2))) < 0.0) {
671  sasymm2 = -10.0;
672  } else {
673  sasymm2 = 2.0 * std::sqrt( a_levdens * (epsilon0_2_saddle -
674  delta_u2*(std::exp((-epsilon_2_saddle*gamma_heavy2))) ) );
675  }
676 
677  if (epsilon_symm_saddle > 0.0) {
678  ssymm = 2.0 * std::sqrt( a_levdens*(epsilon_symm_saddle) );
679  } else {
680  ssymm = -10.0;
681  }
682 
683  if (ssymm > -10.0) {
684  ysymm = 1.0;
685  if (epsilon0_1_saddle < 0.0) { // /* low energy */
686  yasymm1 = std::exp((sasymm1-ssymm)) * wzasymm1_saddle / wzsymm_saddle * 2.0;
687  // /* factor of 2 for symmetry classes */
688  } else { // /* high energy */
689  ssymm_mode1 = 2.0 * std::sqrt( a_levdens*(epsilon0_1_saddle) );
690  yasymm1 = ( std::exp((sasymm1-ssymm)) - std::exp((ssymm_mode1 - ssymm)) )
691  * wzasymm1_saddle / wzsymm_saddle * 2.0;
692  }
693 
694  if (epsilon0_2_saddle < 0.0) { // /* low energy */
695  yasymm2 = std::exp((sasymm2-ssymm)) * wzasymm2_saddle / wzsymm_saddle * 2.0;
696  // /* factor of 2 for symmetry classes */
697  } else { // /* high energy */
698  ssymm_mode2 = 2.0 * std::sqrt( a_levdens*(epsilon0_2_saddle) );
699  yasymm2 = ( std::exp((sasymm2-ssymm)) - std::exp((ssymm_mode2 - ssymm)) )
700  * wzasymm2_saddle / wzsymm_saddle * 2.0;
701  }
702  // /* difference in the exponent in order */
703  // /* to avoid numerical overflow */
704  }
705  else {
706  if ( (sasymm1 > -10.0) && (sasymm2 > -10.0) ) {
707  ysymm = 0.0;
708  yasymm1 = std::exp(sasymm1) * wzasymm1_saddle * 2.0;
709  yasymm2 = std::exp(sasymm2) * wzasymm2_saddle * 2.0;
710  }
711  }
712 
713  // /* normalize */
714  ysum = ysymm + yasymm1 + yasymm2;
715  if (ysum > 0.0) {
716  ysymm = ysymm / ysum;
717  yasymm1 = yasymm1 / ysum;
718  yasymm2 = yasymm2 / ysum;
719  } else {
720  ysymm = 0.0;
721  yasymm1 = 0.0;
722  yasymm2 = 0.0;
723  // /* search minimum threshold and attribute all events to this mode */
724  if ( (epsilon_symm_saddle < epsilon_1_saddle) && (epsilon_symm_saddle < epsilon_2_saddle) ) {
725  ysymm = 1.0;
726  } else {
727  if (epsilon_1_saddle < epsilon_2_saddle) {
728  yasymm1 = 1.0;
729  } else {
730  yasymm2 = 1.0;
731  }
732  }
733  }
734 
735 // if (itest == 1) {
736 // // G4cout << "ysymm normalized= " << ysymm << G4endl;
737 // // G4cout << "yasymm1 normalized= " << yasymm1 << G4endl;
738 // // G4cout << "yasymm2 normalized= " << yasymm2 << G4endl;
739 // }
740 
741  // /* even-odd effect */
742  // /* simple parametrization KHS, Nov. 2000. From Rejmund et al. */
743  eexc_max = max(eexc1_saddle, eexc2_saddle);
744  eexc_max = max(eexc_max, e);
745  // // G4cout << "mod(z, 2)" << iz%2 << G4endl;
746  if ((G4int)(z) % 2 == 0) {
747  r_e_o_max = 0.3 * (1.0 - 0.2 * (std::pow(z, 2)/a - std::pow(92.0, 2)/238.0));
748  e_pair = 2.0 * 12.0 / std::sqrt(a);
749  if(eexc_max > (e_crit + e_pair)) {
750  r_e_o = 0.0;
751  } else {
752  if(eexc_max < e_pair) {
753  r_e_o = r_e_o_max;
754  } else {
755  r_e_o = std::pow((eexc_max - e_crit - e_pair)/e_crit, 2) * r_e_o_max;
756  }
757  }
758  } else {
759  r_e_o = 0.0;
760  }
761 
762  // // G4cout <<"rmax " << r_e_o_max << G4endl;
763  // if(r_e_o > 0.0) // G4cout <<"e_crit, r_e_o" << e_crit << r_e_o << G4endl;
764  // $LOOP; /* event loop */
765  // I_COUNT = I_COUNT + 1;
766 
767  /* random decision: symmetric or asymmetric */
768  /* IMODE = 3 means asymmetric fission, mode 1,
769  IMODE = 2 means asymmetric fission, mode 2,
770  IMODE = 1 means symmetric */
771  // RMODE = dble(HAZ(k));
772  // rmode = rnd.rndm();
773 // // Safety check added to make sure we always select well defined
774 // // fission mode.
775  rmode = haz(k);
776  // Cast for test CS 11/10/05
777  // RMODE = 0.54;
778  // rmode = 0.54;
779  if (rmode < ysymm) {
780  imode = 1;
781  } else if (rmode < (ysymm + yasymm1)) {
782  imode = 2;
783  } else {
784  imode = 3;
785  }
786  // /* determine parameters of the Z distribution */
787  // force imode (for testing, PK)
788  // imode = 3;
789 
790  if (imode == 1) {
791  z1mean = zsymm;
792  z1width = wzsymm;
793  } else if (imode == 2) {
794  z1mean = zheavy1;
795  z1width = wzasymm1;
796  } else if (imode == 3) {
797  z1mean = zheavy2;
798  z1width = wzasymm2;
799  }
800 
801  if (itest == 1) {
802  // G4cout << "nbre aleatoire tire " << rmode << G4endl;
803  // G4cout << "fission mode " << imode << G4endl;
804  // G4cout << "z1mean= " << z1mean << G4endl;
805  // G4cout << "z1width= " << z1width << G4endl;
806  }
807 
808  // /* random decision: Z1 and Z2 at scission: */
809  z1 = 1.0;
810  z2 = 1.0;
811 
812  while ( (z1<5.0) || (z2<5.0) ) {
813  // Z1 = dble(GAUSSHAZ(K,sngl(Z1mean),sngl(Z1width)));
814  // z1 = rnd.gaus(z1mean,z1width);
815  // z1 = 48.26; // gausshaz(k, z1mean, z1width);
816  z1 = gausshaz(k, z1mean, z1width);
817  even_odd(z1, r_e_o, i_help);
818  z1 = double(i_help);
819  z2 = z - z1;
820  }
821 
822  if (itest == 1) {
823  // G4cout << "ff charge sample " << G4endl;
824  // G4cout << "z1 = " << z1 << G4endl;
825  // G4cout << "z2 = " << z2 << G4endl;
826  }
827 
828 // // CALL EVEN_ODD(Z1,R_E_O,I_HELP);
829 // // /* Integer proton number with even-odd effect */
830 // // Z1 = REAL(I_HELP)
831 // // /* Z1 = INT(Z1+0.5E0); */
832 // z2 = z - z1;
833 
834  // /* average N of both fragments: */
835  if (imode == 1) {
836  n1ucd = z1 * n/z;
837  n2ucd = z2 * n/z;
838  re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) + ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,2.0); // umass == massdef
839  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);
840  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);
841  reps2 = (re1-2.0*re2+re3) / 2.0;
842  reps1 = re2 - re1 - reps2;
843  rn1_pol = - reps1 / (2.0 * reps2);
844  n1mean = n1ucd + rn1_pol;
845  } else {
846  n1mean = (z1 + 0.5) * n/z;
847  }
848 
849 // n1mean nsymm + (z1 - zsymm) * 1.6 from 238 U(nth, f)
850 // n1width = 0.9 + E * 0.002 KHS
851 
852 // random decision: N1R and N2R at scission, before evaporation
853  n1r = 1.0;
854  n2r = 1.0;
855  while (n1r < 5 || n2r < 5) {
856  // n1r = 76.93; gausshaz(kkk,n1mean,n1width);
857  n1r = gausshaz(kkk,n1mean,n1width);
858  // modification to have n1r as integer, and n=n1r+n2r rigorously a.b. 19/4/2001
859  i_inter = int(n1r + 0.5);
860  n1r = double(i_inter);
861  n2r = n - n1r;
862  }
863 
864  // neutron evaporation from fragments
865  if (i_eva > 0) {
866  // treatment sz
867  ne_min = 0.095e0 * a - 20.4e0;
868  if (ne_min < 0) ne_min = 0.0;
869  ne_min = ne_min + e / 8.e0; // 1 neutron per 8 mev */
870  }
871 
872  // excitation energy due to deformation
873 
874  a1 = z1 + n1r; // mass of first fragment */
875  a2 = z2 + n2r; // mass of second fragment */
876  if (a1 < 80) {
877  ed1_low = 0.0;
878  } else if (a1 >= 80 && a1 < 110) {
879  ed1_low = (a1-80.)*20./30.;
880  } else if (a1 >= 110 && a1 < 130) {
881  ed1_low = -(a1-110.)*20./20. + 20.;
882  } else if (a1 >= 130) {
883  ed1_low = (a1-130.)*20./30.;
884  }
885 
886  if (a2 < 80) {
887  ed2_low = 0.0;
888  } else if (a2 >= 80 && a2 < 110) {
889  ed2_low = (a2-80.)*20./30.;
890  } else if (a2 >= 110 && a2 < 130) {
891  ed2_low = -(a2-110.)*20./20. + 20.;
892  } else if (a2 >= 130) {
893  ed2_low = (a2-130.)*20./30.;
894  }
895 
896  ed1_high = 20.0*a1/(a1+a2);
897  ed2_high = 20.0 - ed1_high;
898  ed1 = ed1_low*std::exp(-e/el) + ed1_high*(1-std::exp(-e/el));
899  ed2 = ed2_low*std::exp(-e/el) + ed2_high*(1-std::exp(-e/el));
900 
901  // write(6,101)e,a1,a2,ed1,ed2,ed1+ed2
902  // write(6,102)ed1_low,ed1_high,ed2_low,ed2_high
903  e1 = e*a1/(a1+a2) + ed1;
904  e2 = e - e*a1/(a1+a2) + ed2;
905  atot = a1+a2;
906  if (atot > a+1) {
907  // write(6,*)'a,,a1,a2,atot',a,a1,a2,atot
908  // write(6,*)'n,n1r,n2r',n,n1r,n2r
909  // write(6,*)'z,z1,z2',z,z1,z2
910  }
911 
912  milledeux:
913  // only symmetric fission
914  // Symmetric fission: Ok! Checked CS 10/10/05
915  if ( (icz == -1) || (a1 < 0.0) || (a2 < 0.0) ) {
916  // IF (z.eq.92) THEN
917  // write(6,*)'symmetric fission'
918  // write(6,*)'Z,A,E,A1,A2,icz,Atot',Z,A,E,A1,A2,icz,Atot
919  // END IF
920 
921  if (itest == 1) {
922  // G4cout << "milledeux: liquid-drop option " << G4endl;
923  }
924 
925  n = a-z;
926  // proton number in symmetric fission (centre) *
927  zsymm = z / 2.0;
928  nsymm = n / 2.0;
929  asymm = nsymm + zsymm;
930 
931  a_levdens = a / 8.0;
932 
933  masscurv = 2.0;
934  cz_symm = 8.0 / std::pow(z,2) * masscurv;
935 
936  wzsymm = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cz_symm) ) ;
937 
938  if (itest == 1) {
939  // G4cout << " symmetric high energy fission " << G4endl;
940  // G4cout << "wzsymm " << wzsymm << G4endl;
941  }
942 
943  z1mean = zsymm;
944  z1width = wzsymm;
945 
946  // random decision: Z1 and Z2 at scission: */
947  z1 = 1.0;
948  z2 = 1.0;
949  while ( (z1 < 5.0) || (z2 < 5.0) ) {
950  // z1 = dble(gausshaz(kkk,sngl(z1mean),sngl(z1width)));
951  // z1 = rnd.gaus(z1mean,z1width);
952  // z1 = 24.8205585; //gausshaz(kkk, z1mean, z1width);
953  z1 = gausshaz(kkk, z1mean, z1width);
954  z2 = z - z1;
955  }
956 
957  if (itest == 1) {
958  // G4cout << " z1 " << z1 << G4endl;
959  // G4cout << " z2 " << z2 << G4endl;
960  }
961  if (itest == 1) {
962  // G4cout << " zsymm " << zsymm << G4endl;
963  // G4cout << " nsymm " << nsymm << G4endl;
964  // G4cout << " asymm " << asymm << G4endl;
965  }
966 
967  cn = umass(zsymm, nsymm+1.0, 0.0) + umass(zsymm, nsymm-1.0, 0.0)
968  + 1.44 * std::pow(zsymm, 2)/
969  (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))
970  - 2.0 * umass(zsymm, nsymm, 0.0) - 1.44 * std::pow(zsymm, 2) /
971  std::pow(r_null * 2.0 *std::pow(asymm, 1.0/3.0), 2);
972  // This is an approximation! Coulomb energy is neglected.
973 
974  n1width = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cn) + 1.35);
975  if (itest == 1) {
976  // G4cout << " cn " << cn << G4endl;
977  // G4cout << " n1width " << n1width << G4endl;
978  }
979 
980  // /* average N of both fragments: */
981  n1ucd = z1 * n/z;
982  n2ucd = z2 * n/z;
983  re1 = umass(z1,n1ucd, 0.6) + umass(z2,n2ucd, 0.6) + ecoul(z1,n1ucd, 0.6,z2,n2ucd, 0.6,2.0);
984  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);
985  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);
986  reps2 = (re1-2.0*re2+re3) / 2.0;
987  reps1 = re2 - re1 - reps2;
988  rn1_pol = - reps1 / (2.0 * reps2);
989  n1mean = n1ucd + rn1_pol;
990 
991  // random decision: N1R and N2R at scission, before evaporation: */
992  // N1R = dfloat(NINT(GAUSSHAZ(KKK,sngl(N1mean),sngl(N1width))));
993  // n1r = (float)( (int)(rnd.gaus(n1mean,n1width)) );
994  // n1r = 34.0; //(float)( (int)(gausshaz(k, n1mean,n1width)) );
995  n1r = (float)( (int)(gausshaz(k, n1mean,n1width)) );
996  n2r = n - n1r;
997  // Mass of first and second fragment */
998  a1 = z1 + n1r;
999  a2 = z2 + n2r;
1000 
1001  e1 = e*a1/(a1+a2);
1002  e2 = e - e*a1/(a1+a2);
1003  }
1004  v1 = 0.0; // These are not calculated in SimFis3.
1005  v2 = 0.0;
1006  if (itest == 1) {
1007  // G4cout << " n1r " << n1r << G4endl;
1008  // G4cout << " n2r " << n2r << G4endl;
1009  }
1010 
1011  if (itest == 1) {
1012  // G4cout << " a1 " << a1 << G4endl;
1013  // G4cout << " z1 " << z1 << G4endl;
1014  // G4cout << " a2 " << a2 << G4endl;
1015  // G4cout << " z2 " << z2 << G4endl;
1016  // G4cout << " e1 " << e1 << G4endl;
1017  // G4cout << " e2 " << e << G4endl;
1018  }
1019 }
1020 
1022 {
1023  const G4int pSize = 110;
1024  static G4ThreadLocal G4double p[pSize];
1025  static G4ThreadLocal G4long ix = 0, i = 0;
1026  static G4ThreadLocal G4double x = 0.0, y = 0.0, a = 0.0, fhaz = 0.0;
1027  // k =< -1 on initialise
1028  // k = -1 c'est reproductible
1029  // k < -1 || k > -1 ce n'est pas reproductible
1030 
1031  // Zero is invalid random seed. Set proper value from our random seed collection:
1032  if(ix == 0) {
1033  // ix = hazard->ial;
1034  }
1035 
1036  if (k <= -1) { //then
1037  if(k == -1) { //then
1038  ix = 0;
1039  }
1040  else {
1041  x = 0.0;
1042  y = secnds(int(x));
1043  ix = int(y * 100 + 43543000);
1044  if(mod(ix,2) == 0) {
1045  ix = ix + 1;
1046  }
1047  }
1048 
1049  // Here we are using random number generator copied from INCL code
1050  // instead of the CERNLIB one! This causes difficulties for
1051  // automatic testing since the random number generators, and thus
1052  // the behavior of the routines in C++ and FORTRAN versions is no
1053  // longer exactly the same!
1054  x = G4AblaRandom::flat();
1055  // standardRandom(&x, &ix);
1056  for(G4int iRandom = 0; iRandom < pSize; iRandom++) { //do i=1,110
1057  p[iRandom] = G4AblaRandom::flat();
1058  // standardRandom(&(p[i]), &ix);
1059  }
1060  a = G4AblaRandom::flat();
1061  //standardRandom(&a, &ix);
1062  k = 0;
1063  }
1064 
1065  i = nint(100*a)+1;
1066  fhaz = p[i];
1067  a = G4AblaRandom::flat();
1068  // standardRandom(&a, &ix);
1069  p[i] = a;
1070 
1071  // hazard->ial = ix;
1072  // haz=0.4;
1073  return fhaz;
1074 }
1075 
1076 G4double G4AblaFission::gausshaz(int k, double xmoy, double sig)
1077 {
1078  // Gaussian random numbers:
1079 
1080  // 1005 C*** TIRAGE ALEATOIRE DANS UNE GAUSSIENNE DE LARGEUR SIG ET MOYENNE XMOY
1081  static G4ThreadLocal G4int iset = 0;
1082  static G4ThreadLocal G4double v1,v2,r,fac,gset,fgausshaz;
1083 
1084  if(iset == 0) { //then
1085  do {
1086  v1 = 2.0*haz(k) - 1.0;
1087  v2 = 2.0*haz(k) - 1.0;
1088  r = std::pow(v1,2) + std::pow(v2,2);
1089  } while(r >= 1);
1090 
1091  fac = std::sqrt(-2.*std::log(r)/r);
1092  gset = v1*fac;
1093  fgausshaz = v2*fac*sig+xmoy;
1094  iset = 1;
1095  }
1096  else {
1097  fgausshaz=gset*sig+xmoy;
1098  iset=0;
1099  }
1100  return fgausshaz;
1101 }
1102 
1103 // Utilities
1104 
1106 {
1107  if(a < b) {
1108  return a;
1109  }
1110  else {
1111  return b;
1112  }
1113 }
1114 
1116 {
1117  if(a < b) {
1118  return a;
1119  }
1120  else {
1121  return b;
1122  }
1123 }
1124 
1126 {
1127  if(a > b) {
1128  return a;
1129  }
1130  else {
1131  return b;
1132  }
1133 }
1134 
1136 {
1137  if(a > b) {
1138  return a;
1139  }
1140  else {
1141  return b;
1142  }
1143 }
1144 
1146 {
1147  G4double intpart = 0.0;
1148  G4double fractpart = 0.0;
1149  fractpart = std::modf(number, &intpart);
1150  if(number == 0) {
1151  return 0;
1152  }
1153  if(number > 0) {
1154  if(fractpart < 0.5) {
1155  return int(std::floor(number));
1156  }
1157  else {
1158  return int(std::ceil(number));
1159  }
1160  }
1161  if(number < 0) {
1162  if(fractpart < -0.5) {
1163  return int(std::floor(number));
1164  }
1165  else {
1166  return int(std::ceil(number));
1167  }
1168  }
1169 
1170  return int(std::floor(number));
1171 }
1172 
1174 {
1175  time_t mytime;
1176  tm *mylocaltime;
1177 
1178  time(&mytime);
1179  mylocaltime = localtime(&mytime);
1180 
1181  if(x == 0) {
1182  return(mylocaltime->tm_hour*60*60 + mylocaltime->tm_min*60 + mylocaltime->tm_sec);
1183  }
1184  else {
1185  return(mytime - x);
1186  }
1187 }
1188 
1190 {
1191  if(b != 0) {
1192  return (a - (a/b)*b);
1193  }
1194  else {
1195  return 0;
1196  }
1197 }
1198 
1200 {
1201  if(b != 0) {
1202  return (a - (a/b)*b);
1203  }
1204  else {
1205  return 0.0;
1206  }
1207 }
1208 
1210 {
1211  G4double value = 0.0;
1212 
1213  if(a < 0.0) {
1214  value = double(std::ceil(a));
1215  }
1216  else {
1217  value = double(std::floor(a));
1218  }
1219 
1220  return value;
1221 }
1222 
1224 {
1225  G4int value = 0;
1226 
1227  if(a < 0) {
1228  value = int(std::ceil(a));
1229  }
1230  else {
1231  value = int(std::floor(a));
1232  }
1233 
1234  return value;
1235 }
1236 
1238 {
1239  if(a < b && a < c) {
1240  return a;
1241  }
1242  if(b < a && b < c) {
1243  return b;
1244  }
1245  if(c < a && c < b) {
1246  return c;
1247  }
1248  return a;
1249 }
1250 
1252 {
1253  if(a > 0) {
1254  return a;
1255  }
1256  if(a < 0) {
1257  return (-1*a);
1258  }
1259  if(a == 0) {
1260  return a;
1261  }
1262 
1263  return a;
1264 }
1265 
G4int mod(G4int a, G4int b)
static const G4double fac
G4double z
Definition: TRTMaterials.hh:39
G4double dint(G4double a)
static const G4double a1
G4int secnds(G4int x)
const G4double pi
long G4long
Definition: G4Types.hh:80
static const G4double e2
G4double a
Definition: TRTMaterials.hh:39
#define G4ThreadLocal
Definition: tls.hh:52
int G4int
Definition: G4Types.hh:78
G4double umass(G4double z, G4double n, G4double beta)
G4int idint(G4double a)
G4int max(G4int a, G4int b)
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 A[nN]
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)
G4double haz(G4int k)
G4double dmod(G4double a, G4double b)
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