// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // $Id: G4AblaFissionSimfis18.cc,v 1.3 2008/11/06 08:42:00 gcosmo Exp $ // Translation of INCL4.2/ABLA V3 // Pekka Kaitaniemi, HIP (translation) // Christelle Schmidt, IPNL (fission code) // Alain Boudard, CEA (contact person INCL/ABLA) // Aatos Heikkinen, HIP (project coordination) #include "G4AblaFissionSimfis18.hh" #include G4AblaFissionSimfis18::G4AblaFissionSimfis18() { } G4AblaFissionSimfis18::G4AblaFissionSimfis18(G4Hazard *hzr, G4InclRandomInterface *rndm) { hazard = hzr; randomGenerator = rndm; setAboutString("Fission model: Based on ABLA with SimFis18"); } G4AblaFissionSimfis18::~G4AblaFissionSimfis18() { } void G4AblaFissionSimfis18::doFission(G4double &A, G4double &Z, G4double &E, G4double &A1, G4double &Z1, G4double &E1, G4double &K1, G4double &A2, G4double &Z2, G4double &E2, G4double &K2) { fissionDistri(A,Z,E,A1,Z1,E1,K1,A2,Z2,E2,K2); } void G4AblaFissionSimfis18::even_odd(G4double r_origin,G4double r_even_odd,G4int &i_out) { // Procedure to calculate I_OUT from R_IN in a way that // on the average a flat distribution in R_IN results in a // fluctuating distribution in I_OUT with an even-odd effect as // given by R_EVEN_ODD // /* ------------------------------------------------------------ */ // /* EXAMPLES : */ // /* ------------------------------------------------------------ */ // /* If R_EVEN_ODD = 0 : */ // /* CEIL(R_IN) ---- */ // /* */ // /* R_IN -> */ // /* (somewhere in between CEIL(R_IN) and FLOOR(R_IN)) */ */ // /* */ // /* FLOOR(R_IN) ---- --> I_OUT */ // /* ------------------------------------------------------------ */ // /* If R_EVEN_ODD > 0 : */ // /* The interval for the above treatment is */ // /* larger for FLOOR(R_IN) = even and */ // /* smaller for FLOOR(R_IN) = odd */ // /* For R_EVEN_ODD < 0 : just opposite treatment */ // /* ------------------------------------------------------------ */ // /* ------------------------------------------------------------ */ // /* On input: R_ORIGIN nuclear charge (real number) */ // /* R_EVEN_ODD requested even-odd effect */ // /* Intermediate quantity: R_IN = R_ORIGIN + 0.5 */ // /* On output: I_OUT nuclear charge (integer) */ // /* ------------------------------------------------------------ */ // G4double R_ORIGIN,R_IN,R_EVEN_ODD,R_REST,R_HELP; G4double r_in = 0.0, r_rest = 0.0, r_help = 0.0; G4double r_floor = 0.0; G4double r_middle = 0.0; // G4int I_OUT,N_FLOOR; G4int n_floor = 0; r_in = r_origin + 0.5; r_floor = (float)((int)(r_in)); if (r_even_odd < 0.001) { i_out = (int)(r_floor); } else { r_rest = r_in - r_floor; r_middle = r_floor + 0.5; n_floor = (int)(r_floor); if (n_floor%2 == 0) { // even before modif. r_help = r_middle + (r_rest - 0.5) * (1.0 - r_even_odd); } else { // odd before modification r_help = r_middle + (r_rest - 0.5) * (1.0 + r_even_odd); } i_out = (int)(r_help); } } G4double G4AblaFissionSimfis18::umass(G4double z,G4double n,G4double beta) { // liquid-drop mass, Myers & Swiatecki, Lysekil, 1967 // pure liquid drop, without pairing and shell effects // On input: Z nuclear charge of nucleus // N number of neutrons in nucleus // beta deformation of nucleus // On output: binding energy of nucleus G4double a = 0.0, umass = 0.0; G4double alpha = 0.0; G4double xcom = 0.0, xvs = 0.0, xe = 0.0; const G4double pi = 3.1416; a = n + z; alpha = ( std::sqrt(5.0/(4.0*pi)) ) * beta; xcom = 1.0 - 1.7826 * ((a - 2.0*z)/a)*((a - 2.0*z)/a); // factor for asymmetry dependence of surface and volume term xvs = - xcom * ( 15.4941 * a - 17.9439 * std::pow(a,0.66667) * (1.0+0.4*alpha*alpha) ); // sum of volume and surface energy xe = z*z * (0.7053/(std::pow(a,0.33333)) * (1.0-0.2*alpha*alpha) - 1.1529/a); umass = xvs + xe; return umass; } G4double G4AblaFissionSimfis18::ecoul(G4double z1,G4double n1,G4double beta1,G4double z2,G4double n2,G4double beta2,G4double d) { // Coulomb potential between two nuclei // surfaces are in a distance of d // in a tip to tip configuration // approximate formulation // On input: Z1 nuclear charge of first nucleus // N1 number of neutrons in first nucleus // beta1 deformation of first nucleus // Z2 nuclear charge of second nucleus // N2 number of neutrons in second nucleus // beta2 deformation of second nucleus // d distance of surfaces of the nuclei // G4double Z1,N1,beta1,Z2,N2,beta2,d,ecoul; G4double ecoul = 0; G4double dtot = 0; const G4double r0 = 1.16; dtot = r0 * ( std::pow((z1+n1),0.33333) * (1.0+(2.0/3.0)*beta1) + std::pow((z2+n2),0.33333) * (1.0+(2.0/3.0)*beta2) ) + d; ecoul = z1 * z2 * 1.44 / dtot; return ecoul; } void G4AblaFissionSimfis18::fissionDistri(G4double &a,G4double &z,G4double &e, G4double &a1,G4double &z1,G4double &e1,G4double &v1, G4double &a2,G4double &z2,G4double &e2,G4double &v2) { // On input: A, Z, E (mass, atomic number and exc. energy of compound nucleus // before fission) // On output: Ai, Zi, Ei (mass, atomic number and exc. energy of fragment 1 and 2 // after fission) // Additionally calculated but not put in the parameter list: // Kinetic energy of prefragments EkinR1, EkinR2 // Translation of SIMFIS18.PLI (KHS, 2.1.2001) // This program calculates isotopic distributions of fission fragments // with a semiempirical model // Copy from SIMFIS3, KHS, 8. February 1995 // Modifications made by Jose Benlliure and KHS in August 1996 // Energy counted from lowest barrier (J. Benlliure, KHS 1997) // Some bugs corrected (J. Benlliure, KHS 1997) // Version used for thesis S. Steinhaueser (August 1997) // (Curvature of LD potential increased by factor of 2!) // Weiter veraendert mit der Absicht, eine Version zu erhalten, die // derjenigen entspricht, die von J. Benlliure et al. // in Nucl. Phys. A 628 (1998) 458 verwendet wurde, // allerdings ohne volle Neutronenabdampfung. // The excitation energy was calculate now for each fission channel // separately. The dissipation from saddle to scission was taken from // systematics, the deformation energy at scission considers the shell // effects in a simplified way, and the fluctuation is included. // KHS, April 1999 // The width in N/Z was carefully adapted to values given by Lang et al. // The width and eventually a shift in N/Z (polarization) follows the // following rules: // The line N/Z following UCD has an angle of std::atan(Zcn/Ncn) // to the horizontal axis on a chart of nuclides. // (For 238U the angle is 32.2 deg.) // The following relations hold: (from Armbruster) // // sigma(N) (A=const) = sigma(Z) (A=const) // sigma(A) (N=const) = sigma(Z) (N=const) // sigma(A) (Z=const) = sigma(N) (Z=const) // // From this we get: // sigma(Z) (N=const) * N = sigma(N) (Z=const) * Z // sigma(A) (Z=const) = sigma(Z) (A=const) * A/Z // sigma(N) (Z=const) = sigma(Z) (A=const) * A/Z // Z*sigma(N) (Z=const) = N*sigma(Z) (N=const) = A*sigma(Z) (A=const) // Excitation energy now calculated above the lowest potential point // Inclusion of a distribution of excitation energies // Several modifications, starting from SIMFIS12: KHS November 2000 // This version seems to work quite well for 238U. // The transition from symmetric to asymmetric fission around 226Th // is reasonably well reproduced, although St. I is too strong and St. II // is too weak. St. I and St. II are also weakly seen for 208Pb. // Extensions for an event generator of fission events (21.11.2000,KHS) // Defalt parameters (IPARS) rather carefully adjusted to // pre-neutron mass distributions of Vives et al. (238U + n) // Die Parameter Fgamma1 und Fgamma2 sind kleiner als die resultierenden // Breiten der Massenverteilungen!!! // Fgamma1 und Fgamma2 wurden angepa�, so da� // Sigma-A(ST-I) = 3.3, Sigma-A(St-II) = 5.8 (nach Vives) // Parameters of the model carefully adjusted by KHS (2.2.2001) to // 238U + 208Pb, 1000 A MeV, Timo Enqvist et al. G4double n = 0.0; G4double nlight1 = 0.0, nlight2 = 0.0; G4double aheavy1 = 0.0,alight1 = 0.0, aheavy2 = 0.0, alight2 = 0.0; G4double eheavy1 = 0.0, elight1 = 0.0, eheavy2 = 0.0, elight2 = 0.0; G4double zheavy1_shell = 0.0, zheavy2_shell = 0.0; G4double zlight1 = 0.0, zlight2 = 0.0; G4double masscurv = 0.0; G4double sasymm1 = 0.0, sasymm2 = 0.0, ssymm = 0.0, ysum = 0.0, yasymm = 0.0; G4double ssymm_mode1 = 0.0, ssymm_mode2 = 0.0; G4double cz_asymm1_saddle = 0.0, cz_asymm2_saddle = 0.0; // Curvature at saddle, modified by ld-potential G4double wzasymm1_saddle, wzasymm2_saddle, wzsymm_saddle = 0.0; G4double wzasymm1_scission = 0.0, wzasymm2_scission = 0.0, wzsymm_scission = 0.0; G4double wzasymm1 = 0.0, wzasymm2 = 0.0, wzsymm = 0.0; G4double nlight1_eff = 0.0, nlight2_eff = 0.0; G4int imode = 0; G4double rmode = 0.0; G4double z1mean = 0.0, z2mean = 0.0, z1width = 0.0, za1width = 0.0; // G4double Z1,Z2,N1R,N2R,A1R,A2R,N1,N2,A1,A2; G4double n1r = 0.0, n2r = 0.0, a1r = 0.0, a2r = 0.0, n1 = 0.0, n2 = 0.0; G4double zsymm = 0.0, nsymm = 0.0, asymm = 0.0; G4double n1mean = 0.0, n2mean, n1width; G4double dueff = 0.0; // effective shell effect at lowest barrier G4double eld = 0.0; // Excitation energy with respect to ld barrier G4double re1 = 0.0, re2 = 0.0, re3 = 0.0; G4double eps1 = 0.0, eps2 = 0.0; G4double n1ucd = 0.0, n2ucd = 0.0, z1ucd = 0.0, z2ucd = 0.0; G4double beta = 0.0, beta1 = 0.0, beta2 = 0.0; G4double dn1_pol = 0.0; // shift of most probable neutron number for given Z, // according to polarization G4int i_help = 0; // /* Parameters of the semiempirical fission model */ G4double a_levdens = 0.0; // /* level-density parameter */ G4double a_levdens_light1 = 0.0, a_levdens_light2 = 0.0; G4double a_levdens_heavy1 = 0.0, a_levdens_heavy2 = 0.0; const G4double r_null = 1.16; // /* radius parameter */ G4double epsilon_1_saddle = 0.0, epsilon0_1_saddle = 0.0; G4double epsilon_2_saddle = 0.0, epsilon0_2_saddle = 0.0, epsilon_symm_saddle = 0.0; G4double epsilon_1_scission = 0.0, epsilon0_1_scission = 0.0; G4double epsilon_2_scission = 0.0, epsilon0_2_scission = 0.0; G4double epsilon_symm_scission = 0.0; // /* modified energy */ G4double e_eff1_saddle = 0.0, e_eff2_saddle = 0.0; G4double epot0_mode1_saddle = 0.0, epot0_mode2_saddle = 0.0, epot0_symm_saddle = 0.0; G4double epot_mode1_saddle = 0.0, epot_mode2_saddle = 0.0, epot_symm_saddle = 0.0; G4double e_defo = 0.0, e_defo1 = 0.0, e_defo2 = 0.0, e_scission = 0.0, e_asym = 0.0; G4double e1exc = 0.0, e2exc = 0.0; G4double e1exc_sigma = 0.0, e2exc_sigma = 0.0; G4double e1final = 0.0, e2final = 0.0; const G4double r0 = 1.16; G4double tker = 0.0; G4double ekin1 = 0.0, ekin2 = 0.0; // G4double EkinR1,EkinR2,E1,E2,V1,V2; G4double ekinr1 = 0.0, ekinr2 = 0.0; G4int icz = 0, k = 0; // Input parameters: //OMMENT(Nuclear charge number); // G4double Z; //OMMENT(Nuclear mass number); // G4double A; //OMMENT(Excitation energy above fission barrier); // G4double E; // Model parameters: //OMMENT(position of heavy peak valley 1); const G4double nheavy1 = 83.0; //OMMENT(position of heavy peak valley 2); const G4double nheavy2 = 90.0; //OMMENT(Shell effect for valley 1); const G4double delta_u1_shell = -2.65; // Parameter (Delta_U1_shell = -2) //OMMENT(Shell effect for valley 2); const G4double delta_u2_shell = -3.8; // Parameter (Delta_U2_shell = -3.2) //OMMENT(I: used shell effect); G4double delta_u1 = 0.0; //omment(I: used shell effect); G4double delta_u2 = 0.0; //OMMENT(Curvature of asymmetric valley 1); const G4double cz_asymm1_shell = 0.7; //OMMENT(Curvature of asymmetric valley 2); const G4double cz_asymm2_shell = 0.15; //OMMENT(Factor for width of distr. valley 1); const G4double fwidth_asymm1 = 0.63; //OMMENT(Factor for width of distr. valley 2); const G4double fwidth_asymm2 = 0.97; // Parameter (CZ_asymm2_scission = 0.12) //OMMENT(Parameter x: a = A/x); const G4double xlevdens = 12.0; //OMMENT(Factor to gamma_heavy1); const G4double fgamma1 = 2.0; //OMMENT(I: fading of shells (general)); G4double gamma = 0.0; //OMMENT(I: fading of shell 1); G4double gamma_heavy1 = 0.0; //OMMENT(I: fading of shell 2); G4double gamma_heavy2 = 0.0; //OMMENT(Zero-point energy at saddle); const G4double e_zero_point = 0.5; //OMMENT(I: friction from saddle to scission); G4double e_saddle_scission = 0.0; //OMMENT(Friction factor); const G4double friction_factor = 1.0; //OMMENT(I: Internal counter for different modes); INIT(0,0,0) // Integer*4 I_MODE(3) //OMMENT(I: Yield of symmetric mode); G4double ysymm = 0.0; //OMMENT(I: Yield of asymmetric mode 1); G4double yasymm1 = 0.0; //OMMENT(I: Yield of asymmetric mode 2); G4double yasymm2 = 0.0; //OMMENT(I: Effective position of valley 1); G4double nheavy1_eff = 0.0; //OMMENT(I: position of heavy peak valley 1); G4double zheavy1 = 0.0; //omment(I: Effective position of valley 2); G4double nheavy2_eff = 0.0; //OMMENT(I: position of heavy peak valley 2); G4double zheavy2 = 0.0; //omment(I: Excitation energy above saddle 1); G4double eexc1_saddle = 0.0; //omment(I: Excitation energy above saddle 2); G4double eexc2_saddle = 0.0; //omment(I: Excitation energy above lowest saddle); G4double eexc_max = 0.0; //omment(I: Effective mass mode 1); G4double aheavy1_mean = 0.0; //omment(I: Effective mass mode 2); G4double aheavy2_mean = 0.0; //omment(I: Width of symmetric mode); G4double wasymm_saddle = 0.0; //OMMENT(I: Width of asymmetric mode 1); G4double waheavy1_saddle = 0.0; //OMMENT(I: Width of asymmetric mode 2); G4double waheavy2_saddle = 0.0; //omment(I: Width of symmetric mode); G4double wasymm = 0.0; //OMMENT(I: Width of asymmetric mode 1); G4double waheavy1 = 0.0; //OMMENT(I: Width of asymmetric mode 2); G4double waheavy2 = 0.0; //OMMENT(I: Even-odd effect in Z); G4double r_e_o = 0.0, r_e_o_exp = 0.0; //OMMENT(I: Curveture of symmetric valley); G4double cz_symm = 0.0; //OMMENT(I: Curvature of mass distribution for fixed Z); G4double cn = 0.0; //OMMENT(I: Curvature of Z distribution for fixed A); G4double cz = 0.0; //OMMENT(Minimum neutron width for constant Z); const G4double sigzmin = 1.16; //OMMENT(Surface distance of scission configuration); const G4double d = 2.0; // /* Charge polarisation from Wagemanns p. 397: */ //OMMENT(Charge polarisation standard I); const G4double cpol1 = 0.65; //OMMENT(Charge polarisation standard II); const G4double cpol2 = 0.55; //OMMENT(=1: Polarisation simult. in N and Z); const G4int nzpol = 1; //OMMENT(=1: test output, =0: no test output); const G4int itest = 0; // G4double UMASS, ECOUL, reps1, reps2, rn1_pol; G4double reps1 = 0.0, reps2 = 0.0, rn1_pol = 0.0; // Float_t HAZ,GAUSSHAZ; G4int kkk = 0; // G4int kkk = 10; // PK // I_MODE = 0; if(itest == 1) { G4cout << " cn mass " << a << G4endl; G4cout << " cn charge " << z << G4endl; G4cout << " cn energy " << e << G4endl; } // /* average Z of asymmetric and symmetric components: */ n = a - z; /* neutron number of the fissioning nucleus */ k = 0; icz = 0; if ( (std::pow(z,2)/a < 25.0) || (n < nheavy2) || (e > 500.0) ) { icz = -1; // GOTO 1002; goto milledeux; } nlight1 = n - nheavy1; nlight2 = n - nheavy2; // /* Polarisation assumed for standard I and standard II: // Z - Zucd = cpol (for A = const); // from this we get (see Armbruster) // Z - Zucd = Acn/Ncn * cpol (for N = const) */ zheavy1_shell = ((nheavy1/n) * z) - ((a/n) * cpol1); zheavy2_shell = ((nheavy2/n) * z) - ((a/n) * cpol2); e_saddle_scission = (-24.0 + 0.02227 * (std::pow(z,2))/(std::pow(a,0.33333)) ) * friction_factor; // /* Energy dissipated from saddle to scission */ // /* F. Rejmund et al., Nucl. Phys. A 678 (2000) 215, fig. 4 b */ // E_saddle_scission = DMAX1(0.,E_saddle_scission); if (e_saddle_scission > 0.) { e_saddle_scission = e_saddle_scission; } else { e_saddle_scission = 0.; } // /* Semiempirical fission model: */ // /* Fit to experimental result on curvature of potential at saddle */ // /* reference: */ // /* IF Z**2/A < 33.15E0 THEN // MassCurv = 30.5438538E0 - 4.00212049E0 * Z**2/A // + 0.11983384E0 * Z**4 / (A**2) ; // ELSE // MassCurv = 10.E0 ** (7.16993332E0 - 0.26602401E0 * Z**2/A // + 0.00283802E0 * Z**4 / (A**2)) ; */ // /* New parametrization of T. Enqvist according to Mulgin et al. 1998 */ if ( (std::pow(z,2))/a < 34.0) { masscurv = std::pow( 10.0,(-1.093364 + 0.082933 * (std::pow(z,2)/a) - 0.0002602 * (std::pow(z,4)/std::pow(a,2))) ); } else { masscurv = std::pow( 10.0,(3.053536 - 0.056477 * (std::pow(z,2)/a) + 0.0002454 * (std::pow(z,4)/std::pow(a,2))) ); } cz_symm = (8.0/std::pow(z,2)) * masscurv; if(itest == 1) { G4cout << "cz_symmetry= " << cz_symm << G4endl; } if (cz_symm < 0) { icz = -1; // GOTO 1002; goto milledeux; } // /* proton number in symmetric fission (centre) */ zsymm = z/2.0; nsymm = n/2.0; asymm = nsymm + zsymm; zheavy1 = (cz_symm*zsymm + cz_asymm1_shell*zheavy1_shell)/(cz_symm + cz_asymm1_shell); zheavy2 = (cz_symm*zsymm + cz_asymm2_shell*zheavy2_shell)/(cz_symm + cz_asymm2_shell); // /* position of valley due to influence of liquid-drop potential */ nheavy1_eff = (zheavy1 + (a/n * cpol1))*(n/z); nheavy2_eff = (zheavy2 + (a/n * cpol2))*(n/z); nlight1_eff = n - nheavy1_eff; nlight2_eff = n - nheavy2_eff; // /* proton number of light fragments (centre) */ zlight1 = z - zheavy1; // /* proton number of light fragments (centre) */ zlight2 = z - zheavy2; aheavy1 = nheavy1_eff + zheavy1; aheavy2 = nheavy2_eff + zheavy2; aheavy1_mean = aheavy1; aheavy2_mean = aheavy2; alight1 = nlight1_eff + zlight1; alight2 = nlight2_eff + zlight2; a_levdens = a / xlevdens; a_levdens_heavy1 = aheavy1 / xlevdens; a_levdens_heavy2 = aheavy2 / xlevdens; a_levdens_light1 = alight1 / xlevdens; a_levdens_light2 = alight2 / xlevdens; gamma = a_levdens / (0.4 * (std::pow(a,1.3333)) ); gamma_heavy1 = ( a_levdens_heavy1 / (0.4 * (std::pow(aheavy1,1.3333)) ) ) * fgamma1; gamma_heavy2 = a_levdens_heavy2 / (0.4 * (std::pow(aheavy2,1.3333)) ); cz_asymm1_saddle = cz_asymm1_shell + cz_symm; cz_asymm2_saddle = cz_asymm2_shell + cz_symm; // Up to here: Ok! Checked CS 10/10/05 cn = umass(zsymm,(nsymm+1.),0.0) + umass(zsymm,(nsymm-1.),0.0) + 1.44 * (std::pow(zsymm,2))/ ( (std::pow(r_null,2)) * ( std::pow((asymm+1.0),0.33333) + std::pow((asymm-1.0),0.33333) ) * ( std::pow((asymm+1.0),0.33333) + std::pow((asymm-1.0),0.33333) ) ) - 2.0 * umass(zsymm,nsymm,0.0) - 1.44 * (std::pow(zsymm,2))/ ( ( 2.0 * r_null * (std::pow(asymm,0.33333)) ) * ( 2.0 * r_null * (std::pow(asymm,0.33333)) ) ); // /* shell effect in valley of mode 1 */ delta_u1 = delta_u1_shell + (std::pow((zheavy1_shell-zheavy1),2))*cz_asymm1_shell; // /* shell effect in valley of mode 2 */ delta_u2 = delta_u2_shell + (std::pow((zheavy2_shell-zheavy2),2))*cz_asymm2_shell; // /* liquid drop energies // at the centres of the different shell effects // with respect to liquid drop at symmetry: */ epot0_mode1_saddle = (std::pow((zheavy1-zsymm),2)) * cz_symm; epot0_mode2_saddle = (std::pow((zheavy2-zsymm),2)) * cz_symm; epot0_symm_saddle = 0.0; if (itest == 1) { G4cout << "check zheavy1 = " << zheavy1 << G4endl; G4cout << "check zheavy2 = " << zheavy2 << G4endl; G4cout << "check zsymm = " << zsymm << G4endl; G4cout << "check czsymm = " << cz_symm << G4endl; G4cout << "check epot0_mode1_saddle = " << epot0_mode1_saddle << G4endl; G4cout << "check epot0_mode2_saddle = " << epot0_mode2_saddle << G4endl; G4cout << "check epot0_symm_saddle = " << epot0_symm_saddle << G4endl; G4cout << "delta_u1 = " << delta_u1 << G4endl; G4cout << "delta_u2 = " << delta_u2 << G4endl; } // /* energies including shell effects // at the centres of the different shell effects // with respect to liquid drop at symmetry: */ epot_mode1_saddle = epot0_mode1_saddle + delta_u1; epot_mode2_saddle = epot0_mode2_saddle + delta_u2; epot_symm_saddle = epot0_symm_saddle; if (itest == 1) { G4cout << "check epot_mode1_saddle = " << epot_mode1_saddle << G4endl; G4cout << "check epot_mode2_saddle = " << epot_mode2_saddle << G4endl; G4cout << "check epot_symm_saddle = " << epot_symm_saddle << G4endl; } // /* Minimum of potential with respect to ld potential at symmetry */ dueff = min(epot_mode1_saddle,epot_mode2_saddle); dueff = min(dueff,epot_symm_saddle); dueff = dueff - epot_symm_saddle; eld = e + dueff + e_zero_point; if (itest == 1) { G4cout << "check dueff = " << dueff << G4endl; G4cout << "check e = " << e << G4endl; G4cout << "check e_zero_point = " << e_zero_point << G4endl; G4cout << "check eld = " << eld << G4endl; } // Up to here: Ok! Checked CS 10/10/05 // /* E = energy above lowest effective barrier */ // /* Eld = energy above liquid-drop barrier */ // /* Due to this treatment the energy E on input means the excitation */ // /* energy above the lowest saddle. */ // /* These energies are not used */ eheavy1 = e * aheavy1 / a; eheavy2 = e * aheavy2 / a; elight1 = e * alight1 / a; elight2 = e * alight2 / a; epsilon0_1_saddle = eld - e_zero_point - epot0_mode1_saddle; // /* excitation energy at saddle mode 1 without shell effect */ epsilon0_2_saddle = eld - e_zero_point - epot0_mode2_saddle; // /* excitation energy at saddle mode 2 without shell effect */ epsilon_1_saddle = eld - e_zero_point - epot_mode1_saddle; // /* excitation energy at saddle mode 1 with shell effect */ epsilon_2_saddle = eld - e_zero_point - epot_mode2_saddle; // /* excitation energy at saddle mode 2 with shell effect */ epsilon_symm_saddle = eld - e_zero_point - epot_symm_saddle; // /* global parameters */ eexc1_saddle = epsilon_1_saddle; eexc2_saddle = epsilon_2_saddle; eexc_max = max(eexc1_saddle,eexc2_saddle); eexc_max = max(eexc_max,eld); // /* EEXC_MAX is energy above the lowest saddle */ epsilon0_1_scission = eld + e_saddle_scission - epot0_mode1_saddle; // /* excitation energy without shell effect */ epsilon0_2_scission = eld + e_saddle_scission - epot0_mode2_saddle; // /* excitation energy without shell effect */ epsilon_1_scission = eld + e_saddle_scission - epot_mode1_saddle; // /* excitation energy at scission */ epsilon_2_scission = eld+ e_saddle_scission - epot_mode2_saddle; // /* excitation energy at scission */ epsilon_symm_scission = eld + e_saddle_scission - epot_symm_saddle; // /* excitation energy of symmetric fragment at scission */ // /* Calculate widhts at the saddle: */ e_eff1_saddle = epsilon0_1_saddle - delta_u1 * (std::exp((-epsilon_1_saddle*gamma))); if (e_eff1_saddle > 0.0) { wzasymm1_saddle = std::sqrt( (0.5 * (std::sqrt(1.0/a_levdens*e_eff1_saddle)) / (cz_asymm1_shell * std::exp((-epsilon_1_saddle*gamma)) + cz_symm) ) ); } else { wzasymm1_saddle = 1.0; } e_eff2_saddle = epsilon0_2_saddle - delta_u2 * (std::exp((-epsilon_2_saddle*gamma))); if (e_eff2_saddle > 0.0) { wzasymm2_saddle = std::sqrt( (0.5 * (std::sqrt(1.0/a_levdens*e_eff2_saddle)) / (cz_asymm2_shell * std::exp((-epsilon_2_saddle*gamma)) + cz_symm) ) ); } else { wzasymm2_saddle = 1.0; } if (eld > e_zero_point) { if ( (eld + epsilon_symm_saddle) < 0.0) { G4cout << " eld + epsilon_symm_saddle < 0" << G4endl; } wzsymm_saddle = std::sqrt( (0.5 * (std::sqrt(1.0/a_levdens*(eld+epsilon_symm_saddle))) / cz_symm ) ); } else { wzsymm_saddle = 1.0; } if (itest == 1) { G4cout << "wz1(saddle) = " << wzasymm1_saddle << G4endl; G4cout << "wz2(saddle) = " << wzasymm2_saddle << G4endl; G4cout << "wzsymm(saddle) = " << wzsymm_saddle << G4endl; } // /* Calculate widhts at the scission point: */ // /* fits of ref. Beizin 1991 (Plots brought to GSI by Sergei Zhdanov) */ wzsymm_scission = wzsymm_saddle; if (e_saddle_scission == 0.0) { wzasymm1_scission = wzasymm1_saddle; wzasymm2_scission = wzasymm2_saddle; } else { if (nheavy1_eff > 75.0) { wzasymm1_scission = (std::sqrt(21.0)) * z/a; wzasymm2_scission = (std::sqrt (max( (70.0-28.0)/3.0*(z*z/a-35.0)+28.,0.0 )) ) * z/a; } else { wzasymm1_scission = wzasymm1_saddle; wzasymm2_scission = wzasymm2_saddle; } } wzasymm1_scission = max(wzasymm1_scission,wzasymm1_saddle); wzasymm2_scission = max(wzasymm2_scission,wzasymm2_saddle); wzasymm1 = wzasymm1_scission * fwidth_asymm1; wzasymm2 = wzasymm2_scission * fwidth_asymm2; wzsymm = wzsymm_scission; /* if (ITEST == 1) { G4cout << "WZ1(scission) = " << WZasymm1_scission << G4endl; G4cout << "WZ2(scission) = " << WZasymm2_scission << G4endl; G4cout << "WZsymm(scission) = " << WZsymm_scission << G4endl; } if (ITEST == 1) { G4cout << "WZ1(scission) final= " << WZasymm1 << G4endl; G4cout << "WZ2(scission) final= " << WZasymm2 << G4endl; G4cout << "WZsymm(scission) final= " << WZsymm << G4endl; } */ wasymm = wzsymm * a/z; waheavy1 = wzasymm1 * a/z; waheavy2 = wzasymm2 * a/z; wasymm_saddle = wzsymm_saddle * a/z; waheavy1_saddle = wzasymm1_saddle * a/z; waheavy2_saddle = wzasymm2_saddle * a/z; if (itest == 1) { G4cout << "wasymm = " << wzsymm << G4endl; G4cout << "waheavy1 = " << waheavy1 << G4endl; G4cout << "waheavy2 = " << waheavy2 << G4endl; } // Up to here: Ok! Checked CS 11/10/05 if ( (epsilon0_1_saddle - delta_u1*std::exp((-epsilon_1_saddle*gamma_heavy1))) < 0.0) { sasymm1 = -10.0; } else { sasymm1 = 2.0 * std::sqrt( a_levdens * (epsilon0_1_saddle - delta_u1*(std::exp((-epsilon_1_saddle*gamma_heavy1))) ) ); } if ( (epsilon0_2_saddle - delta_u2*std::exp((-epsilon_2_saddle*gamma_heavy2))) < 0.0) { sasymm2 = -10.0; } else { sasymm2 = 2.0 * std::sqrt( a_levdens * (epsilon0_2_saddle - delta_u2*(std::exp((-epsilon_2_saddle*gamma_heavy2))) ) ); } if (epsilon_symm_saddle > 0.0) { ssymm = 2.0 * std::sqrt( a_levdens*(epsilon_symm_saddle) ); } else { ssymm = -10.0; } if (ssymm > -10.0) { ysymm = 1.0; if (epsilon0_1_saddle < 0.0) { // /* low energy */ yasymm1 = std::exp((sasymm1-ssymm)) * wzasymm1_saddle / wzsymm_saddle * 2.0; // /* factor of 2 for symmetry classes */ } else { // /* high energy */ ssymm_mode1 = 2.0 * std::sqrt( a_levdens*(epsilon0_1_saddle) ); yasymm1 = ( std::exp((sasymm1-ssymm)) - std::exp((ssymm_mode1 - ssymm)) ) * wzasymm1_saddle / wzsymm_saddle * 2.0; } if (epsilon0_2_saddle < 0.0) { // /* low energy */ yasymm2 = std::exp((sasymm2-ssymm)) * wzasymm2_saddle / wzsymm_saddle * 2.0; // /* factor of 2 for symmetry classes */ } else { // /* high energy */ ssymm_mode2 = 2.0 * std::sqrt( a_levdens*(epsilon0_2_saddle) ); yasymm2 = ( std::exp((sasymm2-ssymm)) - std::exp((ssymm_mode2 - ssymm)) ) * wzasymm2_saddle / wzsymm_saddle * 2.0; } // /* difference in the exponent in order */ // /* to avoid numerical overflow */ } else { if ( (sasymm1 > -10.0) && (sasymm2 > -10.0) ) { ysymm = 0.0; yasymm1 = std::exp(sasymm1) * wzasymm1_saddle * 2.0; yasymm2 = std::exp(sasymm2) * wzasymm2_saddle * 2.0; } } // /* normalize */ ysum = ysymm + yasymm1 + yasymm2; if (ysum > 0.0) { ysymm = ysymm / ysum; yasymm1 = yasymm1 / ysum; yasymm2 = yasymm2 / ysum; yasymm = yasymm1 + yasymm2; } else { ysymm = 0.0; yasymm1 = 0.0; yasymm2 = 0.0; // /* search minimum threshold and attribute all events to this mode */ if ( (epsilon_symm_saddle < epsilon_1_saddle) && (epsilon_symm_saddle < epsilon_2_saddle) ) { ysymm = 1.0; } else { if (epsilon_1_saddle < epsilon_2_saddle) { yasymm1 = 1.0; } else { yasymm2 = 1.0; } } } if (itest == 1) { G4cout << "ysymm normalized= " << ysymm << G4endl; G4cout << "yasymm1 normalized= " << yasymm1 << G4endl; G4cout << "yasymm2 normalized= " << yasymm2 << G4endl; } // Up to here: Ok! Ckecked CS 11/10/05 // /* even-odd effect */ // /* simple parametrization KHS, Nov. 2000. From Rejmund et al. */ if ((int)(z) % 2 == 0) { r_e_o_exp = -0.017 * (e_saddle_scission + eld) * (e_saddle_scission + eld); if ( r_e_o_exp < -307.0) { r_e_o_exp = -307.0; r_e_o = std::pow(10.0,r_e_o_exp); } else { r_e_o = std::pow(10.0,r_e_o_exp); } } else { r_e_o = 0.0; } // $LOOP; /* event loop */ // I_COUNT = I_COUNT + 1; // /* random decision: symmetric or asymmetric */ // /* IMODE = 1 means asymmetric fission, mode 1, // IMODE = 2 means asymmetric fission, mode 2, // IMODE = 3 means symmetric */ // RMODE = dble(HAZ(k)); // rmode = rnd.rndm(); // Safety check added to make sure we always select well defined // fission mode. do { rmode = haz(k); // Cast for test CS 11/10/05 // RMODE = 0.54; // rmode = 0.54; if (rmode < yasymm1) { imode = 1; } else if ( (rmode > yasymm1) && (rmode < (yasymm1+yasymm2)) ) { imode = 2; } else if ( (rmode > yasymm1) && (rmode > (yasymm1+yasymm2)) ) { imode = 3; } } while(imode == 0); // /* determine parameters of the Z distribution */ // force imode (for testing, PK) // imode = 3; if (imode == 1) { z1mean = zheavy1; z1width = wzasymm1; } if (imode == 2) { z1mean = zheavy2; z1width = wzasymm2; } if (imode == 3) { z1mean = zsymm; z1width = wzsymm; } if (itest == 1) { G4cout << "nbre aleatoire tire " << rmode << G4endl; G4cout << "fission mode " << imode << G4endl; G4cout << "z1mean= " << z1mean << G4endl; G4cout << "z1width= " << z1width << G4endl; } // /* random decision: Z1 and Z2 at scission: */ z1 = 1.0; z2 = 1.0; while ( (z1<5.0) || (z2<5.0) ) { // Z1 = dble(GAUSSHAZ(K,sngl(Z1mean),sngl(Z1width))); // z1 = rnd.gaus(z1mean,z1width); z1 = gausshaz(k, z1mean, z1width); z2 = z - z1; } if (itest == 1) { G4cout << "ff charge sample " << G4endl; G4cout << "z1 = " << z1 << G4endl; G4cout << "z2 = " << z2 << G4endl; } // CALL EVEN_ODD(Z1,R_E_O,I_HELP); // /* Integer proton number with even-odd effect */ // Z1 = REAL(I_HELP) // /* Z1 = INT(Z1+0.5E0); */ z2 = z - z1; // /* average N of both fragments: */ if (imode == 1) { n1mean = (z1 + cpol1 * a/n) * n/z; } if (imode == 2) { n1mean = (z1 + cpol2 * a/n) * n/z; } /* CASE(99) ! only for testing; N1UCD = Z1 * N/Z; N2UCD = Z2 * N/Z; re1 = UMASS(Z1,N1UCD,0.6) +; & UMASS(Z2,N2UCD,0.6) +; & ECOUL(Z1,N1UCD,0.6,Z2,N2UCD,0.6,d); 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); 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); eps2 = (re1-2.0*re2+re3) / 2.0; eps1 = re2 - re1 - eps2; DN1_POL = - eps1 / (2.0 * eps2); N1mean = N1UCD + DN1_POL; */ if (imode == 3) { n1ucd = z1 * n/z; n2ucd = z2 * n/z; re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) + ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,d); 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); 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); eps2 = (re1-2.0*re2+re3) / 2.0; eps1 = re2 - re1 - eps2; dn1_pol = - eps1 / (2.0 * eps2); n1mean = n1ucd + dn1_pol; } // all fission modes features have been checked CS 11/10/05 n2mean = n - n1mean; z2mean = z - z1mean; // /* Excitation energies */ // /* formulated in energies in close consistency with the fission model */ // /* E_defo = UMASS(Z*0.5E0,N*0.5E0,0.6E0) - // UMASS(Z*0.5E0,N*0.5E0,0); */ // /* calculates the deformation energy of the liquid drop for // deformation beta = 0.6 which is most probable at scission */ // /* N1R and N2R provisionaly taken without fluctuations in // polarisation: */ n1r = n1mean; n2r = n2mean; a1r = n1r + z1; a2r = n2r + z2; if (imode == 1) { /* N = 82 */; //! /* Eexc at scission */ e_scission = max(epsilon_1_scission,1.0); if (n1mean > (n * 0.5) ) { //! /* 1. fragment is spherical */ beta1 = 0.0; beta2 = 0.6; e1exc = epsilon_1_scission * a1r / a; e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0); e2exc = epsilon_1_scission * a2r / a + e_defo; } else { //! /* 2. fragment is spherical */ beta1 = 0.6; beta2 = 0.0; e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0); e1exc = epsilon_1_scission * a1r / a + e_defo; e2exc = epsilon_1_scission * a2r / a; } } if (imode == 2) { //! /* N appr. 86 */ e_scission = max(epsilon_2_scission,1.0); if (n1mean > (n * 0.5) ) { //! /* 2. fragment is spherical */ beta1 = (n1r - nheavy2) * 0.034 + 0.3; e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0); e1exc = epsilon_2_scission * a1r / a + e_defo; beta2 = 0.6 - beta1; e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0); e2exc = epsilon_2_scission * a2r / a + e_defo; } else { //! /* 1. fragment is spherical */ beta2 = (n2r - nheavy2) * 0.034 + 0.3; e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0); e2exc = epsilon_2_scission * a2r / a + e_defo; beta1 = 0.6 - beta2; e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0); e1exc = epsilon_2_scission * a1r / a + e_defo; } } if (imode == 3) { // ! /* Symmetric fission channel */ // /* the fit function for beta is the deformation for // optimum energy at the scission point, d = 2 */ // /* beta : deformation of symmetric fragments */ // /* beta1 : deformation of first fragment */ // /* beta2 : deformation of second fragment */ beta = 0.177963 + 0.0153241 * zsymm - 0.000162037 * zsymm*zsymm; beta1 = 0.177963 + 0.0153241 * z1 - 0.000162037 * z1*z1; // beta1 = 0.6 e_defo1 = umass(z1,n1r,beta1) - umass(z1,n1r,0.0); beta2 = 0.177963 + 0.0153241 * z2 - 0.000162037 * z2*z2; // beta2 = 0.6 e_defo2 = umass(z2,n2r,beta2) - umass(z2,n2r,0.0); e_asym = umass(z1 , n1r, beta1) + umass(z2, n2r ,beta2) + ecoul(z1,n1r,beta1,z2,n2r,beta2,2.0) - 2.0 * umass(zsymm,nsymm,beta) - ecoul(zsymm,nsymm,beta,zsymm,nsymm,beta,2.0); // E_asym = CZ_symm * (Z1 - Zsymm)**2 e_scission = max((epsilon_symm_scission - e_asym),1.0); // /* $LIST(Z1,N1R,Z2,N2R,E_asym,E_scission); */ e1exc = e_scission * a1r / a + e_defo1; e2exc = e_scission * a2r / a + e_defo2; } // Energies checked for all the modes CS 11/10/05 // /* random decision: N1R and N2R at scission, before evaporation: */ // /* CN = UMASS(Zsymm , Nsymm + 1.E0,0) + // UMASS(Zsymm, Nsymm - 1.E0,0) // + 1.44E0 * (Zsymm)**2 / // (r_null**2 * ((Asymm+1)**1/3 + (Asymm-1)**1/3)**2 ) // - 2.E0 * UMASS(Zsymm,Nsymm,0) // - 1.44E0 * (Zsymm)**2 / (r_null * 2.E0 * (Asymm)**1/3)**2; */ // /* N1width = std::sqrt(0.5E0 * std::sqrt(1.E0/A_levdens*(Eld+E_saddle_scission)) / CN); */ // /* 8. 9. 1998: KHS (see also consideration in the first comment block) // sigma_N(Z=const) = A/Z * sigma_Z(A=const) // sigma_Z(A=const) = 0.4 to 0.5 (from Lang paper Nucl Phys. A345 (1980) 34) // sigma_N(Z=const) = 0.45 * A/Z (= 1.16 for 238U) // therefore: SIGZMIN = 1.16 */ if ( (imode == 1) || (imode == 2) ) { cn=(umass(z1,n1mean+1.,beta1) + umass(z1,n1mean-1.,beta1) + umass(z2,n2mean+1.,beta2) + umass(z2,n2mean-1.,beta2) + ecoul(z1,n1mean+1.,beta1,z2,n2mean-1.,beta2,2.0) + ecoul(z1,n1mean-1.,beta1,z2,n2mean+1.,beta2,2.0) - 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0) - 2.0 * umass(z1, n1mean, beta1) - 2.0 * umass(z2, n2mean, beta2) ) * 0.5; // /* Coulomb energy neglected for the moment! */ // IF (E_scission.lt.0.) Then // write(6,*)' E_scission < 0, MODE 1,2' // ENDIF // IF (CN.lt.0.) Then // write(6,*)'CN < 0, MODE 1,2' // ENDIF n1width=std::sqrt( (0.5 * (std::sqrt(1.0/a_levdens*(e_scission)))/cn) ); n1width=max(n1width, sigzmin); // /* random decision: N1R and N2R at scission, before evaporation: */ n1r = 1.0; n2r = 1.0; while ( (n1r<5.0) || (n2r<5.0) ) { // n1r = dble(gausshaz(k,sngl(n1mean),sngl(n1width))); // n1r = rnd.gaus(n1mean,n1width); n1r = gausshaz(k, n1mean, n1width); n2r = n - n1r; } // N1R = GAUSSHAZ(K,N1mean,N1width) if (itest == 1) { G4cout << "after neutron sample " << n1r << G4endl; } n1r = (float)( (int)((n1r+0.5)) ); n2r = n - n1r; even_odd(z1,r_e_o,i_help); // /* proton number with even-odd effect */ z1 = (float)(i_help); z2 = z - z1; a1r = z1 + n1r; a2r = z2 + n2r; } if (imode == 3) { //! /* When(3) */ if (nzpol > 0.0) { // /* We treat a simultaneous split in Z and N to determine polarisation */ cz = ( umass(z1-1., n1mean+1.,beta1) + umass(z2+1., n2mean-1.,beta1) + umass(z1+1., n1mean-1.,beta2) + umass(z2 - 1., n2mean + 1.,beta2) + ecoul(z1-1.,n1mean+1.,beta1,z2+1.,n2mean-1.,beta2,2.0) + ecoul(z1+1.,n1mean-1.,beta1,z2-1.,n2mean+1.,beta2,2.0) - 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0) - 2.0 * umass(z1, n1mean,beta1) - 2.0 * umass(z2, n2mean,beta2) ) * 0.5; // IF (E_scission.lt.0.) Then // write(6,*) ' E_scission < 0, MODE 1,2' // ENDIF // IF (CZ.lt.0.) Then // write(6,*) 'CZ < 0, MODE 1,2' // ENDIF za1width=std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*(e_scission)) / cz) ); za1width=std::sqrt( (max((za1width*za1width-(1.0/12.0)),0.1)) ); // /* Check the value of 0.1 ! */ // /* Shephard correction */ a1r = z1 + n1mean; a1r = (float)((int)((a1r+0.5))); a2r = a - a1r; // /* A1R and A2R are integer numbers now */ // /* $LIST(A1R,A2R,ZA1WIDTH); */ n1ucd = n/a * a1r; n2ucd = n/a * a2r; z1ucd = z/a * a1r; z2ucd = z/a * a2r; re1 = umass(z1ucd-1.,n1ucd+1.,beta1) + umass(z2ucd+1.,n2ucd-1.,beta2) + ecoul(z1ucd-1.,n1ucd+1.,beta1,z2ucd+1.,n2ucd-1.,beta2,d); re2 = umass(z1ucd,n1ucd,beta1) + umass(z2ucd,n2ucd,beta2) + ecoul(z1ucd,n1ucd,beta1,z2ucd,n2ucd,beta2,d); re3 = umass(z1ucd+1.,n1ucd-1.,beta1) + umass(z2ucd-1.,n2ucd+1.,beta2) + + ecoul(z1ucd+1.,n1ucd-1.,beta1,z2ucd-1.,n2ucd+1.,beta2,d); eps2 = (re1-2.0*re2+re3) / 2.0; eps1 = (re3 - re1)/2.0; dn1_pol = - eps1 / (2.0 * eps2); z1 = z1ucd + dn1_pol; if (itest == 1) { G4cout << "before proton sample " << z1 << G4endl; } // Z1 = dble(GAUSSHAZ(k,sngl(Z1),sngl(ZA1width))); // z1 = rnd.gaus(z1,za1width); z1 = gausshaz(k, z1, za1width); if (itest == 1) { G4cout << "after proton sample " << z1 << G4endl; } even_odd(z1,r_e_o,i_help); // /* proton number with even-odd effect */ z1 = (float)(i_help); z2 = (float)((int)( (z - z1 + 0.5)) ); n1r = a1r - z1; n2r = n - n1r; } else { // /* First division of protons, then adjustment of neutrons */ cn = ( umass(z1, n1mean+1.,beta1) + umass(z1, n1mean-1., beta1) + umass(z2, n2mean+1.,beta2) + umass(z2, n2mean-1., beta2) + ecoul(z1,n1mean+1.,beta1,z2,n2mean-1.,beta2,2.0) + ecoul(z1,n1mean-1.,beta1,z2,n2mean+1.,beta2,2.0) - 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0) - 2.0 * umass(z1, n1mean, 0.6) - 2.0 * umass(z2, n2mean, 0.6) ) * 0.5; // /* Coulomb energy neglected for the moment! */ // IF (E_scission.lt.0.) Then // write(6,*) ' E_scission < 0, MODE 1,2' // Endif // IF (CN.lt.0.) Then // write(6,*) 'CN < 0, MODE 1,2' // Endif n1width=std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*(e_scission)) / cn) ); n1width=max(n1width, sigzmin); // /* random decision: N1R and N2R at scission, before evaporation: */ // N1R = dble(GAUSSHAZ(k,sngl(N1mean),sngl(N1width))); // n1r = rnd.gaus(n1mean,n1width); n1r = gausshaz(k, n1mean, n1width); n1r = (float)( (int)((n1r+0.5)) ); n2r = n - n1r; even_odd(z1,r_e_o,i_help); // /* Integer proton number with even-odd effect */ z1 = (float)(i_help); z2 = z - z1; a1r = z1 + n1r; a2r = z2 + n2r; } } if (itest == 1) { G4cout << "remid imode = " << imode << G4endl; G4cout << "n1width = " << n1width << G4endl; G4cout << "n1r = " << n1r << G4endl; G4cout << "a1r = " << a1r << G4endl; G4cout << "n2r = " << n2r << G4endl; G4cout << "a2r = " << a2r << G4endl; } // Up to here: checked CS 11/10/05 // /* Extracted from Lang et al. Nucl. Phys. A 345 (1980) 34 */ e1exc_sigma = 5.5; e2exc_sigma = 5.5; neufcentquatrevingtsept: // E1final = dble(Gausshaz(k,sngl(E1exc),sngl(E1exc_sigma))); // E2final = dble(Gausshaz(k,sngl(E2exc),sngl(E2exc_sigma))); // e1final = rnd.gaus(e1exc,e1exc_sigma); // e2final = rnd.gaus(e2exc,e2exc_sigma); e1final = gausshaz(k, e1exc, e1exc_sigma); e2final = gausshaz(k, e2exc, e2exc_sigma); if ( (e1final < 0.0) || (e2final < 0.0) ) goto neufcentquatrevingtsept; if (itest == 1) { G4cout << "sampled exc 1 " << e1final << G4endl; G4cout << "sampled exc 2 " << e2final << G4endl; } // /* OUTPUT QUANTITIES OF THE EVENT GENERATOR: */ // /* Quantities before neutron evaporation */ // /* Neutron number of prefragments: N1R and N2R */ // /* Atomic number of fragments: Z1 and Z2 */ // /* Kinetic energy of fragments: EkinR1, EkinR2 *7 // /* Quantities after neutron evaporation: */ // /* Neutron number of fragments: N1 and N2 */ // /* Mass number of fragments: A1 and A2 */ // /* Atomic number of fragments: Z1 and Z2 */ // /* Number of evaporated neutrons: N1R-N1 and N2R-N2 */ // /* Kinetic energy of fragments: EkinR1*A1/A1R and // EkinR2*A2/A2R */ n1 = n1r; n2 = n2r; a1 = n1 + z1; a2 = n2 + z2; e1 = e1final; e2 = e2final; // /* Pre-neutron-emission total kinetic energy: */ tker = (z1 * z2 * 1.44) / ( r0 * std::pow(a1,0.33333) * (1.0 + 2.0/3.0 * beta1) + r0 * std::pow(a2,0.33333) * (1.0 + 2.0/3.0 * beta2) + 2.0 ); // /* Pre-neutron-emission kinetic energy of 1. fragment: */ ekinr1 = tker * a2 / a; // /* Pre-neutron-emission kinetic energy of 2. fragment: */ ekinr2 = tker * a1 / a; v1 = std::sqrt( (ekinr1/a1) ) * 1.3887; v2 = std::sqrt( (ekinr2/a2) ) * 1.3887; if (itest == 1) { G4cout << "ekinr1 " << ekinr1 << G4endl; G4cout << "ekinr2 " << ekinr2 << G4endl; } milledeux: //************************** //*** only symmetric fission //************************** // Symmetric fission: Ok! Checked CS 10/10/05 if ( (icz == -1) || (a1 < 0.0) || (a2 < 0.0) ) { // IF (z.eq.92) THEN // write(6,*)'symmetric fission' // write(6,*)'Z,A,E,A1,A2,icz,Atot',Z,A,E,A1,A2,icz,Atot // END IF if (itest == 1) { G4cout << "milledeux: liquid-drop option " << G4endl; } n = a-z; // proton number in symmetric fission (centre) * zsymm = z / 2.0; nsymm = n / 2.0; asymm = nsymm + zsymm; a_levdens = a / xlevdens; masscurv = 2.0; cz_symm = 8.0 / std::pow(z,2) * masscurv; wzsymm = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cz_symm) ) ; if (itest == 1) { G4cout << " symmetric high energy fission " << G4endl; G4cout << "wzsymm " << wzsymm << G4endl; } z1mean = zsymm; z1width = wzsymm; // random decision: Z1 and Z2 at scission: */ z1 = 1.0; z2 = 1.0; while ( (z1 < 5.0) || (z2 < 5.0) ) { // z1 = dble(gausshaz(kkk,sngl(z1mean),sngl(z1width))); // z1 = rnd.gaus(z1mean,z1width); z1 = gausshaz(kkk, z1mean, z1width); z2 = z - z1; } if (itest == 1) { G4cout << " z1 " << z1 << G4endl; G4cout << " z2 " << z2 << G4endl; } if (itest == 1) { G4cout << " zsymm " << zsymm << G4endl; G4cout << " nsymm " << nsymm << G4endl; G4cout << " asymm " << asymm << G4endl; } // CN = UMASS(Zsymm , Nsymm + 1.E0) + UMASS(Zsymm, Nsymm - 1.E0) // # + 1.44E0 * (Zsymm)**2 / // # (r_null**2 * ((Asymm+1)**(1./3.) + // # (Asymm-1)**(1./3.))**2 ) // # - 2.E0 * UMASS(Zsymm,Nsymm) // # - 1.44E0 * (Zsymm)**2 / // # (r_null * 2.E0 * (Asymm)**(1./3.))**2 n1ucd = z1 * n/z; n2ucd = z2 * n/z; re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) + ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,2.0); 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); 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); reps2 = (re1-2.0*re2+re3)/2.0; reps1 = re2 - re1 -reps2; rn1_pol = -reps1/(2.0*reps2); n1mean = n1ucd + rn1_pol; n2mean = n - n1mean; if (itest == 1) { G4cout << " n1mean " << n1mean << G4endl; G4cout << " n2mean " << n2mean << G4endl; } cn = (umass(z1,n1mean+1.,0.0) + umass(z1,n1mean-1.,0.0) + + umass(z2,n2mean+1.,0.0) + umass(z2,n2mean-1.,0.0) - 2.0 * umass(z1,n1mean,0.0) + - 2.0 * umass(z2,n2mean,0.0) ) * 0.5; // This is an approximation! Coulomb energy is neglected. n1width = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cn) ); if (itest == 1) { G4cout << " cn " << cn << G4endl; G4cout << " n1width " << n1width << G4endl; } // random decision: N1R and N2R at scission, before evaporation: */ // N1R = dfloat(NINT(GAUSSHAZ(KKK,sngl(N1mean),sngl(N1width)))); // n1r = (float)( (int)(rnd.gaus(n1mean,n1width)) ); n1r = (float)( (int)(gausshaz(k, n1mean,n1width)) ); n2r = n - n1r; // Mass of first and second fragment */ a1 = z1 + n1r; a2 = z2 + n2r; e1 = e*a1/(a1+a2); e2 = e - e*a1/(a1+a2); if (itest == 1) { G4cout << " n1r " << n1r << G4endl; G4cout << " n2r " << n2r << G4endl; } } if (itest == 1) { G4cout << " a1 " << a1 << G4endl; G4cout << " z1 " << z1 << G4endl; G4cout << " a2 " << a2 << G4endl; G4cout << " z2 " << z2 << G4endl; G4cout << " e1 " << e1 << G4endl; G4cout << " e2 " << e << G4endl; } // /* Pre-neutron-emission total kinetic energy: */ tker = (z1 * z2 * 1.44) / ( r0 * std::pow(a1,0.33333) * (1.0 + 2.0/3.0 * beta1) + r0 * std::pow(a2,0.33333) * (1.0 + 2.0/3.0 * beta2) + 2.0 ); // /* Pre-neutron-emission kinetic energy of 1. fragment: */ ekin1 = tker * a2 / a; // /* Pre-neutron-emission kinetic energy of 2. fragment: */ ekin2 = tker * a1 / a; v1 = std::sqrt( (ekin1/a1) ) * 1.3887; v2 = std::sqrt( (ekin2/a2) ) * 1.3887; if (itest == 1) { G4cout << " kinetic energies " << G4endl; G4cout << " ekin1 " << ekin1 << G4endl; G4cout << " ekin2 " << ekin2 << G4endl; } } void G4AblaFissionSimfis18::standardRandom(G4double *rndm, G4long *seed) { (*seed) = (*seed); // Avoid warning during compilation. // Use Geant4 G4UniformRand (*rndm) = randomGenerator->getRandom(); } G4double G4AblaFissionSimfis18::haz(G4int k) { const G4int pSize = 110; static G4double p[pSize]; static G4long ix = 0, i = 0; static G4double x = 0.0, y = 0.0, a = 0.0, haz = 0.0; // k =< -1 on initialise // k = -1 c'est reproductible // k < -1 || k > -1 ce n'est pas reproductible // Zero is invalid random seed. Set proper value from our random seed collection: if(ix == 0) { ix = hazard->ial; } if (k <= -1) { //then if(k == -1) { //then ix = 0; } else { x = 0.0; y = secnds(int(x)); ix = int(y * 100 + 43543000); if(mod(ix,2) == 0) { ix = ix + 1; } } // Here we are using random number generator copied from INCL code // instead of the CERNLIB one! This causes difficulties for // automatic testing since the random number generators, and thus // the behavior of the routines in C++ and FORTRAN versions is no // longer exactly the same! x = randomGenerator->getRandom(); // standardRandom(&x, &ix); for(G4int i = 0; i < pSize; i++) { //do i=1,110 p[i] = randomGenerator->getRandom(); // standardRandom(&(p[i]), &ix); } a = randomGenerator->getRandom(); standardRandom(&a, &ix); k = 0; } i = nint(100*a)+1; haz = p[i]; a = randomGenerator->getRandom(); // standardRandom(&a, &ix); p[i] = a; hazard->ial = ix; return haz; } G4double G4AblaFissionSimfis18::gausshaz(int k, double xmoy, double sig) { // Gaussian random numbers: // 1005 C*** TIRAGE ALEATOIRE DANS UNE GAUSSIENNE DE LARGEUR SIG ET MOYENNE XMOY static G4int iset = 0; static G4double v1,v2,r,fac,gset,gausshaz; if(iset == 0) { //then do { v1 = 2.0*haz(k) - 1.0; v2 = 2.0*haz(k) - 1.0; r = std::pow(v1,2) + std::pow(v2,2); } while(r >= 1); fac = std::sqrt(-2.*std::log(r)/r); gset = v1*fac; gausshaz = v2*fac*sig+xmoy; iset = 1; } else { gausshaz=gset*sig+xmoy; iset=0; } return gausshaz; } // Utilities G4double G4AblaFissionSimfis18::min(G4double a, G4double b) { if(a < b) { return a; } else { return b; } } G4int G4AblaFissionSimfis18::min(G4int a, G4int b) { if(a < b) { return a; } else { return b; } } G4double G4AblaFissionSimfis18::max(G4double a, G4double b) { if(a > b) { return a; } else { return b; } } G4int G4AblaFissionSimfis18::max(G4int a, G4int b) { if(a > b) { return a; } else { return b; } } G4int G4AblaFissionSimfis18::nint(G4double number) { G4double intpart = 0.0; G4double fractpart = 0.0; fractpart = std::modf(number, &intpart); if(number == 0) { return 0; } if(number > 0) { if(fractpart < 0.5) { return int(std::floor(number)); } else { return int(std::ceil(number)); } } if(number < 0) { if(fractpart < -0.5) { return int(std::floor(number)); } else { return int(std::ceil(number)); } } return int(std::floor(number)); } G4int G4AblaFissionSimfis18::secnds(G4int x) { time_t mytime; tm *mylocaltime; time(&mytime); mylocaltime = localtime(&mytime); if(x == 0) { return(mylocaltime->tm_hour*60*60 + mylocaltime->tm_min*60 + mylocaltime->tm_sec); } else { return(mytime - x); } } G4int G4AblaFissionSimfis18::mod(G4int a, G4int b) { if(b != 0) { return (a - (a/b)*b); } else { return 0; } } G4double G4AblaFissionSimfis18::dmod(G4double a, G4double b) { if(b != 0) { return (a - (a/b)*b); } else { return 0.0; } } G4double G4AblaFissionSimfis18::dint(G4double a) { G4double value = 0.0; if(a < 0.0) { value = double(std::ceil(a)); } else { value = double(std::floor(a)); } return value; } G4int G4AblaFissionSimfis18::idint(G4double a) { G4int value = 0; if(a < 0) { value = int(std::ceil(a)); } else { value = int(std::floor(a)); } return value; } G4int G4AblaFissionSimfis18::idnint(G4double value) { G4double valueCeil = int(std::ceil(value)); G4double valueFloor = int(std::floor(value)); if(std::fabs(value - valueCeil) < std::fabs(value - valueFloor)) { return int(valueCeil); } else { return int(valueFloor); } } G4double G4AblaFissionSimfis18::dmin1(G4double a, G4double b, G4double c) { if(a < b && a < c) { return a; } if(b < a && b < c) { return b; } if(c < a && c < b) { return c; } return a; } G4double G4AblaFissionSimfis18::utilabs(G4double a) { if(a > 0) { return a; } if(a < 0) { return (-1*a); } if(a == 0) { return a; } return a; }