Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4AblaFission Class Reference

#include <G4AblaFission.hh>

Inheritance diagram for G4AblaFission:
G4AblaFissionBase

Public Member Functions

 G4AblaFission ()
 
 ~G4AblaFission ()
 
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 spdef (G4int a, G4int z, G4int optxfis)
 
G4double fissility (G4int a, G4int z, G4int optxfis)
 
void even_odd (G4double r_origin, G4double r_even_odd, G4int &i_out)
 
G4double umass (G4double z, G4double n, G4double beta)
 
G4double ecoul (G4double z1, G4double n1, G4double beta1, G4double z2, G4double n2, G4double beta2, G4double d)
 
void fissionDistri (G4double &a, G4double &z, G4double &e, G4double &a1, G4double &z1, G4double &e1, G4double &v1, G4double &a2, G4double &z2, G4double &e2, G4double &v2)
 
void standardRandom (G4double *rndm, G4long *seed)
 
G4double haz (G4int k)
 
G4double gausshaz (int k, double xmoy, double sig)
 
G4int min (G4int a, G4int b)
 
G4double min (G4double a, G4double b)
 
G4int max (G4int a, G4int b)
 
G4double max (G4double a, G4double b)
 
G4int nint (G4double number)
 
G4int secnds (G4int x)
 
G4int mod (G4int a, G4int b)
 
G4double dmod (G4double a, G4double b)
 
G4double dint (G4double a)
 
G4int idint (G4double a)
 
G4double utilabs (G4double a)
 
G4double dmin1 (G4double a, G4double b, G4double c)
 
- Public Member Functions inherited from G4AblaFissionBase
 G4AblaFissionBase ()
 
virtual ~G4AblaFissionBase ()
 
void setVerboseLevel (G4int level)
 
void about ()
 
void setAboutString (G4String anAbout)
 

Detailed Description

Definition at line 42 of file G4AblaFission.hh.

Constructor & Destructor Documentation

G4AblaFission::G4AblaFission ( )

Definition at line 40 of file G4AblaFission.cc.

41 {
42 }
G4AblaFission::~G4AblaFission ( )

Definition at line 44 of file G4AblaFission.cc.

45 {
46 }

Member Function Documentation

G4double G4AblaFission::dint ( G4double  a)

Definition at line 1209 of file G4AblaFission.cc.

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 }
value
const XML_Char int const XML_Char * value
Definition: include/expat.h:331
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AblaFission::dmin1 ( G4double  a,
G4double  b,
G4double  c 
)

Definition at line 1237 of file G4AblaFission.cc.

References test::a, test::b, and test::c.

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 }
G4double G4AblaFission::dmod ( G4double  a,
G4double  b 
)

Definition at line 1199 of file G4AblaFission.cc.

1200 {
1201  if(b != 0) {
1202  return (a - (a/b)*b);
1203  }
1204  else {
1205  return 0.0;
1206  }
1207 }
void G4AblaFission::doFission ( G4double &  A,
G4double &  Z,
G4double &  E,
G4double &  A1,
G4double &  Z1,
G4double &  E1,
G4double &  K1,
G4double &  A2,
G4double &  Z2,
G4double &  E2,
G4double &  K2 
)
virtual

Implements G4AblaFissionBase.

Definition at line 48 of file G4AblaFission.cc.

References fissionDistri().

51 {
52  fissionDistri(A,Z,E,A1,Z1,E1,K1,A2,Z2,E2,K2);
53 }
G4AblaFission::fissionDistri
void fissionDistri(G4double &a, G4double &z, G4double &e, G4double &a1, G4double &z1, G4double &e1, G4double &v1, G4double &a2, G4double &z2, G4double &e2, G4double &v2)
Definition: G4AblaFission.cc:171
G4double G4AblaFission::ecoul ( G4double  z1,
G4double  n1,
G4double  beta1,
G4double  z2,
G4double  n2,
G4double  beta2,
G4double  d 
)

Definition at line 144 of file G4AblaFission.cc.

Referenced by fissionDistri().

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 }
G4double
double G4double
Definition: G4Types.hh:76
void G4AblaFission::even_odd ( G4double  r_origin,
G4double  r_even_odd,
G4int &  i_out 
)

Definition at line 55 of file G4AblaFission.cc.

References int().

Referenced by fissionDistri().

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 }
int
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
G4int
int G4int
Definition: G4Types.hh:78
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AblaFission::fissility ( G4int  a,
G4int  z,
G4int  optxfis 
)
void G4AblaFission::fissionDistri ( G4double &  a,
G4double &  z,
G4double &  e,
G4double &  a1,
G4double &  z1,
G4double &  e1,
G4double &  v1,
G4double &  a2,
G4double &  z2,
G4double &  e2,
G4double &  v2 
)

Definition at line 171 of file G4AblaFission.cc.

References test::a, ecoul(), even_odd(), gausshaz(), haz(), int(), max(), n, umass(), and z.

Referenced by doFission().

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 }
int
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
z
G4double z
Definition: TRTMaterials.hh:39
G4int
int G4int
Definition: G4Types.hh:78
G4AblaFission::umass
G4double umass(G4double z, G4double n, G4double beta)
Definition: G4AblaFission.cc:115
G4AblaFission::max
G4int max(G4int a, G4int b)
Definition: G4AblaFission.cc:1135
G4AblaFission::gausshaz
G4double gausshaz(int k, double xmoy, double sig)
Definition: G4AblaFission.cc:1076
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4AblaFission::ecoul
G4double ecoul(G4double z1, G4double n1, G4double beta1, G4double z2, G4double n2, G4double beta2, G4double d)
Definition: G4AblaFission.cc:144
G4AblaFission::even_odd
void even_odd(G4double r_origin, G4double r_even_odd, G4int &i_out)
Definition: G4AblaFission.cc:55
G4AblaFission::haz
G4double haz(G4int k)
Definition: G4AblaFission.cc:1021
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AblaFission::gausshaz ( int  k,
double  xmoy,
double  sig 
)

Definition at line 1076 of file G4AblaFission.cc.

References G4ThreadLocal, gset(), and haz().

Referenced by fissionDistri().

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 }
gset
subroutine gset(AX, BX, NX, Z, W)
Definition: dpm25nulib.f:689
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:52
G4int
int G4int
Definition: G4Types.hh:78
G4AblaFission::haz
G4double haz(G4int k)
Definition: G4AblaFission.cc:1021
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AblaFission::haz ( G4int  k)

Definition at line 1021 of file G4AblaFission.cc.

References test::a, G4AblaRandom::flat(), G4ThreadLocal, int(), mod(), nint(), secnds(), and test::x.

Referenced by fissionDistri(), and gausshaz().

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 }
G4AblaFission::mod
G4int mod(G4int a, G4int b)
Definition: G4AblaFission.cc:1189
int
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
p
const char * p
Definition: xmltok.h:285
G4AblaFission::secnds
G4int secnds(G4int x)
Definition: G4AblaFission.cc:1173
G4long
long G4long
Definition: G4Types.hh:80
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:52
G4int
int G4int
Definition: G4Types.hh:78
G4AblaRandom::flat
double flat()
Definition: G4AblaRandom.cc:47
G4double
double G4double
Definition: G4Types.hh:76
G4AblaFission::nint
G4int nint(G4double number)
Definition: G4AblaFission.cc:1145
G4int G4AblaFission::idint ( G4double  a)

Definition at line 1223 of file G4AblaFission.cc.

References int().

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 }
int
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
G4int
int G4int
Definition: G4Types.hh:78
value
const XML_Char int const XML_Char * value
Definition: include/expat.h:331
G4int G4AblaFission::max ( G4int  a,
G4int  b 
)

Definition at line 1135 of file G4AblaFission.cc.

References test::a, and test::b.

Referenced by fissionDistri().

1136 {
1137  if(a > b) {
1138  return a;
1139  }
1140  else {
1141  return b;
1142  }
1143 }
G4double G4AblaFission::max ( G4double  a,
G4double  b 
)

Definition at line 1125 of file G4AblaFission.cc.

References test::a, and test::b.

1126 {
1127  if(a > b) {
1128  return a;
1129  }
1130  else {
1131  return b;
1132  }
1133 }
G4int G4AblaFission::min ( G4int  a,
G4int  b 
)

Definition at line 1115 of file G4AblaFission.cc.

References test::a, and test::b.

1116 {
1117  if(a < b) {
1118  return a;
1119  }
1120  else {
1121  return b;
1122  }
1123 }
G4double G4AblaFission::min ( G4double  a,
G4double  b 
)

Definition at line 1105 of file G4AblaFission.cc.

References test::a, and test::b.

1106 {
1107  if(a < b) {
1108  return a;
1109  }
1110  else {
1111  return b;
1112  }
1113 }
G4int G4AblaFission::mod ( G4int  a,
G4int  b 
)

Definition at line 1189 of file G4AblaFission.cc.

Referenced by haz().

1190 {
1191  if(b != 0) {
1192  return (a - (a/b)*b);
1193  }
1194  else {
1195  return 0;
1196  }
1197 }
G4int G4AblaFission::nint ( G4double  number)

Definition at line 1145 of file G4AblaFission.cc.

References int().

Referenced by haz().

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 }
int
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
G4double
double G4double
Definition: G4Types.hh:76
G4int G4AblaFission::secnds ( G4int  x)

Definition at line 1173 of file G4AblaFission.cc.

References G4InuclParticleNames::tm.

Referenced by haz().

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 }
G4double G4AblaFission::spdef ( G4int  a,
G4int  z,
G4int  optxfis 
)
void G4AblaFission::standardRandom ( G4double *  rndm,
G4long *  seed 
)
G4double G4AblaFission::umass ( G4double  z,
G4double  n,
G4double  beta 
)

Definition at line 115 of file G4AblaFission.cc.

References test::a, python.hepunit::pi, and z.

Referenced by fissionDistri().

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 }
z
G4double z
Definition: TRTMaterials.hh:39
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4double
double G4double
Definition: G4Types.hh:76
G4double G4AblaFission::utilabs ( G4double  a)

Definition at line 1251 of file G4AblaFission.cc.

References test::a.

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 }
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:10 for Geant4.10 by   doxygen 1.8.7