source: trunk/source/processes/hadronic/models/incl/src/G4AblaFissionSimfis18.cc@ 1199

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

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