#include <G4Abla.hh>
Definition at line 42 of file G4Abla.hh.
| G4Abla::G4Abla | ( | ) |
| G4Abla::G4Abla | ( | G4Hazard * | aHazard, | |
| G4Volant * | aVolant, | |||
| G4VarNtp * | aVarntp | |||
| ) |
This constructor is used by standalone test driver and the Geant4 interface.
| aHazard | random seeds | |
| aVolant | data structure for ABLA output | |
| aVarNtp | data structure for transfering ABLA output to Geant4 interface |
Definition at line 66 of file G4Abla.cc.
References G4Volant::iv.
00067 { 00068 verboseLevel = 0; 00069 ilast = 0; 00070 volant = aVolant; // ABLA internal particle data 00071 volant->iv = 0; 00072 hazard = aHazard; // Random seeds 00073 varntp = aVarntp; // Output data structure 00074 varntp->ntrack = 0; 00075 00076 pace = new G4Pace(); 00077 ald = new G4Ald(); 00078 ablamain = new G4Ablamain(); 00079 emdpar = new G4Emdpar(); 00080 eenuc = new G4Eenuc(); 00081 ec2sub = new G4Ec2sub(); 00082 ecld = new G4Ecld(); 00083 fb = new G4Fb(); 00084 fiss = new G4Fiss(); 00085 opt = new G4Opt(); 00086 }
| G4Abla::G4Abla | ( | G4Hazard * | hazard, | |
| G4Volant * | volant | |||
| ) |
Constructor that is to be used only for testing purposes.
| aHazard | random seeds | |
| aVolant | data structure for ABLA output |
Definition at line 45 of file G4Abla.cc.
References G4Volant::iv.
00046 { 00047 verboseLevel = 0; 00048 ilast = 0; 00049 volant = volant; // ABLA internal particle data 00050 volant->iv = 0; 00051 hazard = hazard; // Random seeds 00052 00053 varntp = new G4VarNtp(); 00054 pace = new G4Pace(); 00055 ald = new G4Ald(); 00056 ablamain = new G4Ablamain(); 00057 emdpar = new G4Emdpar(); 00058 eenuc = new G4Eenuc(); 00059 ec2sub = new G4Ec2sub(); 00060 ecld = new G4Ecld(); 00061 fb = new G4Fb(); 00062 fiss = new G4Fiss(); 00063 opt = new G4Opt(); 00064 }
| G4Abla::~G4Abla | ( | ) |
Procedure for calculating the pairing correction to the binding energy of a specific nucleus.
Definition at line 2707 of file G4Abla.cc.
References parite().
02708 { 02709 // CALCUL DE LA CORRECTION, DUE A L'APPARIEMENT, DE L'ENERGIE DE 02710 // LIAISON D'UN NOYAU 02711 // PROCEDURE FOR CALCULATING THE PAIRING CORRECTION TO THE BINDING 02712 // ENERGY OF A SPECIFIC NUCLEUS 02713 02714 double para = 0.0, parz = 0.0; 02715 // A MASS NUMBER 02716 // Z NUCLEAR CHARGE 02717 // PARA HELP VARIABLE FOR PARITY OF A 02718 // PARZ HELP VARIABLE FOR PARITY OF Z 02719 // DEL PAIRING CORRECTION 02720 02721 parite(a, ¶); 02722 02723 if (para < 0.0) { 02724 (*del) = 0.0; 02725 } 02726 else { 02727 parite(z, &parz); 02728 if (parz > 0.0) { 02729 // assert(isnan(std::sqrt(a)) == false); 02730 (*del) = -12.0/std::sqrt(a); 02731 } 02732 else { 02733 // assert(isnan(std::sqrt(a)) == false); 02734 (*del) = 12.0/std::sqrt(a); 02735 } 02736 } 02737 }
| void G4Abla::barfit | ( | G4int | iz, | |
| G4int | ia, | |||
| G4int | il, | |||
| G4double * | sbfis, | |||
| G4double * | segs, | |||
| G4double * | selmax | |||
| ) |
THIS SUBROUTINE RETURNS THE BARRIER HEIGHT BFIS, THE GROUND-STATE ENERGY SEGS, IN MEV, AND THE ANGULAR MOMENTUM AT WHICH THE FISSION BARRIER DISAPPEARS, LMAX, IN UNITS OF H-BAR, WHEN CALLED WITH INTEGER AGUMENTS IZ, THE ATOMIC NUMBER, IA, THE ATOMIC MASS NUMBER, AND IL, THE ANGULAR MOMENTUM IN UNITS OF H-BAR. (PLANCK'S CONSTANT DIVIDED BY 2*PI).
Definition at line 2870 of file G4Abla.cc.
References lpoly().
Referenced by direct(), and evapora().
02871 { 02872 // 2223 C VERSION FOR 32BIT COMPUTER 02873 // 2224 C THIS SUBROUTINE RETURNS THE BARRIER HEIGHT BFIS, THE 02874 // 2225 C GROUND-STATE ENERGY SEGS, IN MEV, AND THE ANGULAR MOMENTUM 02875 // 2226 C AT WHICH THE FISSION BARRIER DISAPPEARS, LMAX, IN UNITS OF 02876 // 2227 C H-BAR, WHEN CALLED WITH INTEGER AGUMENTS IZ, THE ATOMIC 02877 // 2228 C NUMBER, IA, THE ATOMIC MASS NUMBER, AND IL, THE ANGULAR 02878 // 2229 C MOMENTUM IN UNITS OF H-BAR. (PLANCK'S CONSTANT DIVIDED BY 02879 // 2230 C 2*PI). 02880 // 2231 C 02881 // 2232 C THE FISSION BARRIER FO IL = 0 IS CALCULATED FROM A 7TH 02882 // 2233 C ORDER FIT IN TWO VARIABLES TO 638 CALCULATED FISSION 02883 // 2234 C BARRIERS FOR Z VALUES FROM 20 TO 110. THESE 638 BARRIERS ARE 02884 // 2235 C FIT WITH AN RMS DEVIATION OF 0.10 MEV BY THIS 49-PARAMETER 02885 // 2236 C FUNCTION. 02886 // 2237 C IF BARFIT IS CALLED WITH (IZ,IA) VALUES OUTSIDE THE RANGE OF 02887 // 2238 C THE BARRIER HEIGHT IS SET TO 0.0, AND A MESSAGE IS PRINTED 02888 // 2239 C ON THE DEFAULT OUTPUT FILE. 02889 // 2240 C 02890 // 2241 C FOR IL VALUES NOT EQUAL TO ZERO, THE VALUES OF L AT WHICH 02891 // 2242 C THE BARRIER IS 80% AND 20% OF THE L=0 VALUE ARE RESPECTIVELY 02892 // 2243 C FIT TO 20-PARAMETER FUNCTIONS OF Z AND A, OVER A MORE 02893 // 2244 C RESTRICTED RANGE OF A VALUES, THAN IS THE CASE FOR L = 0. 02894 // 2245 C THE VALUE OF L WHERE THE BARRIER DISAPPEARS, LMAX IS FIT TO 02895 // 2246 C A 24-PARAMETER FUNCTION OF Z AND A, WITH THE SAME RANGE OF 02896 // 2247 C Z AND A VALUES AS L-80 AND L-20. 02897 // 2248 C ONCE AGAIN, IF AN (IZ,IA) PAIR IS OUTSIDE OF THE RANGE OF 02898 // 2249 C VALIDITY OF THE FIT, THE BARRIER VALUE IS SET TO 0.0 AND A 02899 // 2250 C MESSAGE IS PRINTED. THESE THREE VALUES (BFIS(L=0),L-80, AND 02900 // 2251 C L-20) AND THE CONSTRINTS OF BFIS = 0 AND D(BFIS)/DL = 0 AT 02901 // 2252 C L = LMAX AND L=0 LEAD TO A FIFTH-ORDER FIT TO BFIS(L) FOR 02902 // 2253 C L>L-20. THE FIRST THREE CONSTRAINTS LEAD TO A THIRD-ORDER FIT 02903 // 2254 C FOR THE REGION L < L-20. 02904 // 2255 C 02905 // 2256 C THE GROUND STATE ENERGIES ARE CALCULATED FROM A 02906 // 2257 C 120-PARAMETER FIT IN Z, A, AND L TO 214 GROUND-STATE ENERGIES 02907 // 2258 C FOR 36 DIFFERENT Z AND A VALUES. 02908 // 2259 C (THE RANGE OF Z AND A IS THE SAME AS FOR L-80, L-20, AND 02909 // 2260 C L-MAX) 02910 // 2261 C 02911 // 2262 C THE CALCULATED BARRIERS FROM WHICH THE FITS WERE MADE WERE 02912 // 2263 C CALCULATED IN 1983-1984 BY A. J. SIERK OF LOS ALAMOS 02913 // 2264 C NATIONAL LABORATORY GROUP T-9, USING YUKAWA-PLUS-EXPONENTIAL 02914 // 2265 C G4DOUBLE FOLDED NUCLEAR ENERGY, EXACT COULOMB DIFFUSENESS 02915 // 2266 C CORRECTIONS, AND DIFFUSE-MATTER MOMENTS OF INERTIA. 02916 // 2267 C THE PARAMETERS OF THE MODEL R-0 = 1.16 FM, AS 21.13 MEV, 02917 // 2268 C KAPPA-S = 2.3, A = 0.68 FM. 02918 // 2269 C THE DIFFUSENESS OF THE MATTER AND CHARGE DISTRIBUTIONS USED 02919 // 2270 C CORRESPONDS TO A SURFACE DIFFUSENESS PARAMETER (DEFINED BY 02920 // 2271 C MYERS) OF 0.99 FM. THE CALCULATED BARRIERS FOR L = 0 ARE 02921 // 2272 C ACCURATE TO A LITTLE LESS THAN 0.1 MEV; THE OUTPUT FROM 02922 // 2273 C THIS SUBROUTINE IS A LITTLE LESS ACCURATE. WORST ERRORS MAY BE 02923 // 2274 C AS LARGE AS 0.5 MEV; CHARACTERISTIC UNCERTAINY IS IN THE RANGE 02924 // 2275 C OF 0.1-0.2 MEV. THE RMS DEVIATION OF THE GROUND-STATE FIT 02925 // 2276 C FROM THE 214 INPUT VALUES IS 0.20 MEV. THE MAXIMUM ERROR 02926 // 2277 C OCCURS FOR LIGHT NUCLEI IN THE REGION WHERE THE GROUND STATE 02927 // 2278 C IS PROLATE, AND MAY BE GREATER THAN 1.0 MEV FOR VERY NEUTRON 02928 // 2279 C DEFICIENT NUCLEI, WITH L NEAR LMAX. FOR MOST NUCLEI LIKELY TO 02929 // 2280 C BE ENCOUNTERED IN REAL EXPERIMENTS, THE MAXIMUM ERROR IS 02930 // 2281 C CLOSER TO 0.5 MEV, AGAIN FOR LIGHT NUCLEI AND L NEAR LMAX. 02931 // 2282 C 02932 // 2283 C WRITTEN BY A. J. SIERK, LANL T-9 02933 // 2284 C VERSION 1.0 FEBRUARY, 1984 02934 // 2285 C 02935 // 2286 C THE FOLLOWING IS NECESSARY FOR 32-BIT MACHINES LIKE DEC VAX, 02936 // 2287 C IBM, ETC 02937 02938 G4double pa[7],pz[7],pl[10]; 02939 for(G4int init_i = 0; init_i < 7; init_i++) { 02940 pa[init_i] = 0.0; 02941 pz[init_i] = 0.0; 02942 } 02943 for(G4int init_i = 0; init_i < 10; init_i++) { 02944 pl[init_i] = 0.0; 02945 } 02946 02947 G4double a = 0.0, z = 0.0, amin = 0.0, amax = 0.0, amin2 = 0.0; 02948 G4double amax2 = 0.0, aa = 0.0, zz = 0.0, bfis = 0.0; 02949 G4double bfis0 = 0.0, ell = 0.0, el = 0.0, egs = 0.0, el80 = 0.0, el20 = 0.0; 02950 G4double elmax = 0.0, sel80 = 0.0, sel20 = 0.0, x = 0.0, y = 0.0, q = 0.0, qa = 0.0, qb = 0.0; 02951 G4double aj = 0.0, ak = 0.0, a1 = 0.0, a2 = 0.0; 02952 02953 G4int i = 0, j = 0, k = 0, m = 0; 02954 G4int l = 0; 02955 02956 G4double emncof[4][5] = {{-9.01100e+2,-1.40818e+3, 2.77000e+3,-7.06695e+2, 8.89867e+2}, 02957 {1.35355e+4,-2.03847e+4, 1.09384e+4,-4.86297e+3,-6.18603e+2}, 02958 {-3.26367e+3, 1.62447e+3, 1.36856e+3, 1.31731e+3, 1.53372e+2}, 02959 {7.48863e+3,-1.21581e+4, 5.50281e+3,-1.33630e+3, 5.05367e-2}}; 02960 02961 G4double elmcof[4][5] = {{1.84542e+3,-5.64002e+3, 5.66730e+3,-3.15150e+3, 9.54160e+2}, 02962 {-2.24577e+3, 8.56133e+3,-9.67348e+3, 5.81744e+3,-1.86997e+3}, 02963 {2.79772e+3,-8.73073e+3, 9.19706e+3,-4.91900e+3, 1.37283e+3}, 02964 {-3.01866e+1, 1.41161e+3,-2.85919e+3, 2.13016e+3,-6.49072e+2}}; 02965 02966 G4double emxcof[4][6] = {{9.43596e4,-2.241997e5,2.223237e5,-1.324408e5,4.68922e4,-8.83568e3}, 02967 {-1.655827e5,4.062365e5,-4.236128e5,2.66837e5,-9.93242e4,1.90644e4}, 02968 {1.705447e5,-4.032e5,3.970312e5,-2.313704e5,7.81147e4,-1.322775e4}, 02969 {-9.274555e4,2.278093e5,-2.422225e5,1.55431e5,-5.78742e4,9.97505e3}}; 02970 02971 G4double elzcof[7][7] = {{5.11819909e+5,-1.30303186e+6, 1.90119870e+6,-1.20628242e+6, 5.68208488e+5, 5.48346483e+4,-2.45883052e+4}, 02972 {-1.13269453e+6, 2.97764590e+6,-4.54326326e+6, 3.00464870e+6, -1.44989274e+6,-1.02026610e+5, 6.27959815e+4}, 02973 {1.37543304e+6,-3.65808988e+6, 5.47798999e+6,-3.78109283e+6, 1.84131765e+6, 1.53669695e+4,-6.96817834e+4}, 02974 {-8.56559835e+5, 2.48872266e+6,-4.07349128e+6, 3.12835899e+6, -1.62394090e+6, 1.19797378e+5, 4.25737058e+4}, 02975 {3.28723311e+5,-1.09892175e+6, 2.03997269e+6,-1.77185718e+6, 9.96051545e+5,-1.53305699e+5,-1.12982954e+4}, 02976 {4.15850238e+4, 7.29653408e+4,-4.93776346e+5, 6.01254680e+5, -4.01308292e+5, 9.65968391e+4,-3.49596027e+3}, 02977 {-1.82751044e+5, 3.91386300e+5,-3.03639248e+5, 1.15782417e+5, -4.24399280e+3,-6.11477247e+3, 3.66982647e+2}}; 02978 02979 const G4int sizex = 5; 02980 const G4int sizey = 6; 02981 const G4int sizez = 4; 02982 02983 G4double egscof[sizey][sizey][sizez]; 02984 02985 G4double egs1[sizey][sizex] = {{1.927813e5, 7.666859e5, 6.628436e5, 1.586504e5,-7.786476e3}, 02986 {-4.499687e5,-1.784644e6,-1.546968e6,-4.020658e5,-3.929522e3}, 02987 {4.667741e5, 1.849838e6, 1.641313e6, 5.229787e5, 5.928137e4}, 02988 {-3.017927e5,-1.206483e6,-1.124685e6,-4.478641e5,-8.682323e4}, 02989 {1.226517e5, 5.015667e5, 5.032605e5, 2.404477e5, 5.603301e4}, 02990 {-1.752824e4,-7.411621e4,-7.989019e4,-4.175486e4,-1.024194e4}}; 02991 02992 G4double egs2[sizey][sizex] = {{-6.459162e5,-2.903581e6,-3.048551e6,-1.004411e6,-6.558220e4}, 02993 {1.469853e6, 6.564615e6, 6.843078e6, 2.280839e6, 1.802023e5}, 02994 {-1.435116e6,-6.322470e6,-6.531834e6,-2.298744e6,-2.639612e5}, 02995 {8.665296e5, 3.769159e6, 3.899685e6, 1.520520e6, 2.498728e5}, 02996 {-3.302885e5,-1.429313e6,-1.512075e6,-6.744828e5,-1.398771e5}, 02997 {4.958167e4, 2.178202e5, 2.400617e5, 1.167815e5, 2.663901e4}}; 02998 02999 G4double egs3[sizey][sizex] = {{3.117030e5, 1.195474e6, 9.036289e5, 6.876190e4,-6.814556e4}, 03000 {-7.394913e5,-2.826468e6,-2.152757e6,-2.459553e5, 1.101414e5}, 03001 {7.918994e5, 3.030439e6, 2.412611e6, 5.228065e5, 8.542465e3}, 03002 {-5.421004e5,-2.102672e6,-1.813959e6,-6.251700e5,-1.184348e5}, 03003 {2.370771e5, 9.459043e5, 9.026235e5, 4.116799e5, 1.001348e5}, 03004 {-4.227664e4,-1.738756e5,-1.795906e5,-9.292141e4,-2.397528e4}}; 03005 03006 G4double egs4[sizey][sizex] = {{-1.072763e5,-5.973532e5,-6.151814e5, 7.371898e4, 1.255490e5}, 03007 {2.298769e5, 1.265001e6, 1.252798e6,-2.306276e5,-2.845824e5}, 03008 {-2.093664e5,-1.100874e6,-1.009313e6, 2.705945e5, 2.506562e5}, 03009 {1.274613e5, 6.190307e5, 5.262822e5,-1.336039e5,-1.115865e5}, 03010 {-5.715764e4,-2.560989e5,-2.228781e5,-3.222789e3, 1.575670e4}, 03011 {1.189447e4, 5.161815e4, 4.870290e4, 1.266808e4, 2.069603e3}}; 03012 03013 for(i = 0; i < sizey; i++) { 03014 for(j = 0; j < sizex; j++) { 03015 // egscof[i][j][0] = egs1[i][j]; 03016 // egscof[i][j][1] = egs2[i][j]; 03017 // egscof[i][j][2] = egs3[i][j]; 03018 // egscof[i][j][3] = egs4[i][j]; 03019 egscof[i][j][0] = egs1[i][j]; 03020 egscof[i][j][1] = egs2[i][j]; 03021 egscof[i][j][2] = egs3[i][j]; 03022 egscof[i][j][3] = egs4[i][j]; 03023 } 03024 } 03025 03026 // the program starts here 03027 if (iz < 19 || iz > 111) { 03028 goto barfit900; 03029 } 03030 03031 if(iz > 102 && il > 0) { 03032 goto barfit902; 03033 } 03034 03035 z=double(iz); 03036 a=double(ia); 03037 el=double(il); 03038 amin= 1.2e0*z + 0.01e0*z*z; 03039 amax= 5.8e0*z - 0.024e0*z*z; 03040 03041 if(a < amin || a > amax) { 03042 goto barfit910; 03043 } 03044 03045 // angul.mom.zero barrier 03046 aa=2.5e-3*a; 03047 zz=1.0e-2*z; 03048 ell=1.0e-2*el; 03049 bfis0 = 0.0; 03050 lpoly(zz,7,pz); 03051 lpoly(aa,7,pa); 03052 03053 for(i = 0; i < 7; i++) { //do 10 i=1,7 03054 for(j = 0; j < 7; j++) { //do 10 j=1,7 03055 bfis0=bfis0+elzcof[j][i]*pz[i]*pa[j]; 03056 //bfis0=bfis0+elzcof[i][j]*pz[j]*pa[i]; 03057 } 03058 } 03059 03060 bfis=bfis0; 03061 // assert(isnan(bfis) == false); 03062 03063 (*sbfis)=bfis; 03064 egs=0.0; 03065 (*segs)=egs; 03066 03067 // values of l at which the barrier 03068 // is 20%(el20) and 80%(el80) of l=0 value 03069 amin2 = 1.4e0*z + 0.009e0*z*z; 03070 amax2 = 20.e0 + 3.0e0*z; 03071 03072 if((a < amin2-5.e0 || a > amax2+10.e0) && il > 0) { 03073 goto barfit920; 03074 } 03075 03076 lpoly(zz,5,pz); 03077 lpoly(aa,4,pa); 03078 el80=0.0; 03079 el20=0.0; 03080 elmax=0.0; 03081 03082 for(i = 0; i < 4; i++) { 03083 for(j = 0; j < 5; j++) { 03084 // el80 = el80 + elmcof[j][i]*pz[j]*pa[i]; 03085 // el20 = el20 + emncof[j][i]*pz[j]*pa[i]; 03086 el80 = el80 + elmcof[i][j]*pz[j]*pa[i]; 03087 el20 = el20 + emncof[i][j]*pz[j]*pa[i]; 03088 } 03089 } 03090 03091 sel80 = el80; 03092 sel20 = el20; 03093 03094 // value of l (elmax) where barrier disapp. 03095 lpoly(zz,6,pz); 03096 lpoly(ell,9,pl); 03097 03098 for(i = 0; i < 4; i++) { //do 30 i= 1,4 03099 for(j = 0; j < 6; j++) { //do 30 j=1,6 03100 //elmax = elmax + emxcof[j][i]*pz[j]*pa[i]; 03101 // elmax = elmax + emxcof[j][i]*pz[i]*pa[j]; 03102 elmax = elmax + emxcof[i][j]*pz[j]*pa[i]; 03103 } 03104 } 03105 03106 // assert(isnan(elmax) == false); 03107 (*selmax)=elmax; 03108 03109 // value of barrier at ang.mom. l 03110 if(il < 1){ 03111 return; 03112 } 03113 03114 x = sel20/(*selmax); 03115 // assert(isnan(x) == false); 03116 y = sel80/(*selmax); 03117 // assert(isnan(y) == false); 03118 03119 if(el <= sel20) { 03120 // low l 03121 q = 0.2/(std::pow(sel20,2)*std::pow(sel80,2)*(sel20-sel80)); 03122 qa = q*(4.0*std::pow(sel80,3) - std::pow(sel20,3)); 03123 qb = -q*(4.0*std::pow(sel80,2) - std::pow(sel20,2)); 03124 bfis = bfis*(1.0 + qa*std::pow(el,2) + qb*std::pow(el,3)); 03125 } 03126 else { 03127 // high l 03128 aj = (-20.0*std::pow(x,5) + 25.e0*std::pow(x,4) - 4.0)*std::pow((y-1.0),2)*y*y; 03129 ak = (-20.0*std::pow(y,5) + 25.0*std::pow(y,4) - 1.0) * std::pow((x-1.0),2)*x*x; 03130 q = 0.2/(std::pow((y-x)*((1.0-x)*(1.0-y)*x*y),2)); 03131 qa = q*(aj*y - ak*x); 03132 qb = -q*(aj*(2.0*y + 1.0) - ak*(2.0*x + 1.0)); 03133 z = el/(*selmax); 03134 a1 = 4.0*std::pow(z,5) - 5.0*std::pow(z,4) + 1.0; 03135 a2 = qa*(2.e0*z + 1.e0); 03136 bfis=bfis*(a1 + (z - 1.e0)*(a2 + qb*z)*z*z*(z - 1.e0)); 03137 } 03138 03139 if(bfis <= 0.0) { 03140 bfis=0.0; 03141 } 03142 03143 if(el > (*selmax)) { 03144 bfis=0.0; 03145 } 03146 (*sbfis)=bfis; 03147 03148 // now calculate rotating ground state energy 03149 if(el > (*selmax)) { 03150 return; 03151 } 03152 03153 for(k = 0; k < 4; k++) { 03154 for(l = 0; l < 6; l++) { 03155 for(m = 0; m < 5; m++) { 03156 //egs = egs + egscof[l][m][k]*pz[l]*pa[k]*pl[2*m-1]; 03157 egs = egs + egscof[l][m][k]*pz[l]*pa[k]*pl[2*m]; 03158 // egs = egs + egscof[m][l][k]*pz[l]*pa[k]*pl[2*m-1]; 03159 } 03160 } 03161 } 03162 03163 (*segs)=egs; 03164 if((*segs) < 0.0) { 03165 (*segs)=0.0; 03166 } 03167 03168 return; 03169 03170 barfit900: //continue 03171 (*sbfis)=0.0; 03172 // for z<19 sbfis set to 1.0e3 03173 if (iz < 19) { 03174 (*sbfis) = 1.0e3; 03175 } 03176 (*segs)=0.0; 03177 (*selmax)=0.0; 03178 return; 03179 03180 barfit902: 03181 (*sbfis)=0.0; 03182 (*segs)=0.0; 03183 (*selmax)=0.0; 03184 return; 03185 03186 barfit910: 03187 (*sbfis)=0.0; 03188 (*segs)=0.0; 03189 (*selmax)=0.0; 03190 return; 03191 03192 barfit920: 03193 (*sbfis)=0.0; 03194 (*segs)=0.0; 03195 (*selmax)=0.0; 03196 return; 03197 }
This subroutine calculates the fission barriers of the liquid-drop model of Myers and Swiatecki (1967). Analytic parameterization of Dahlinger 1982 replaces tables. Barrier heights from Myers and Swiatecki
Definition at line 2507 of file G4Abla.cc.
02508 { 02509 // This subroutine calculates the fission barriers 02510 // of the liquid-drop model of Myers and Swiatecki (1967). 02511 // Analytic parameterization of Dahlinger 1982 02512 // replaces tables. Barrier heights from Myers and Swiatecki !!! 02513 02514 G4double nms = 0.0, ims = 0.0, ksims = 0.0, xms = 0.0, ums = 0.0; 02515 02516 nms = ams - zms; 02517 ims = (nms-zms)/ams; 02518 ksims= 50.15e0 * (1.- 1.78 * std::pow(ims,2)); 02519 xms = std::pow(zms,2) / (ams * ksims); 02520 ums = 0.368e0-5.057e0*xms+8.93e0*std::pow(xms,2)-8.71*std::pow(xms,3); 02521 return(0.7322e0*std::pow(zms,2)/std::pow(ams,(0.333333e0))*std::pow(10.e0,ums)); 02522 }
CALCULATION OF THE SURFACE BS OR CURVATURE BK OF A NUCLEUS RELATIVE TO THE SPHERICAL CONFIGURATION BASED ON MYERS, DROPLET MODEL FOR ARBITRARY SHAPES
Definition at line 2813 of file G4Abla.cc.
References G4cout, G4endl, and idint().
Referenced by direct().
02814 { 02815 // CALCULATION OF THE SURFACE BS OR CURVATURE BK OF A NUCLEUS 02816 // RELATIVE TO THE SPHERICAL CONFIGURATION 02817 // BASED ON MYERS, DROPLET MODEL FOR ARBITRARY SHAPES 02818 02819 // INPUT: IFLAG - 0/1 BK/BS CALCULATION 02820 // Y - (1 - X) COMPLEMENT OF THE FISSILITY 02821 02822 // LINEAR INTERPOLATION OF BS BK TABLE 02823 02824 int i = 0; 02825 02826 G4double bipolResult = 0.0; 02827 02828 const int bsbkSize = 54; 02829 02830 G4double bk[bsbkSize] = {0.0, 1.00000,1.00087,1.00352,1.00799,1.01433,1.02265,1.03306, 02831 1.04576,1.06099,1.07910,1.10056,1.12603,1.15651,1.19348, 02832 1.23915,1.29590,1.35951,1.41013,1.44103,1.46026,1.47339, 02833 1.48308,1.49068,1.49692,1.50226,1.50694,1.51114,1.51502, 02834 1.51864,1.52208,1.52539,1.52861,1.53177,1.53490,1.53803, 02835 1.54117,1.54473,1.54762,1.55096,1.55440,1.55798,1.56173, 02836 1.56567,1.56980,1.57413,1.57860,1.58301,1.58688,1.58688, 02837 1.58688,1.58740,1.58740, 0.0}; //Zeroes at bk[0], and at the end added by PK 02838 02839 G4double bs[bsbkSize] = {0.0, 1.00000,1.00086,1.00338,1.00750,1.01319, 02840 1.02044,1.02927,1.03974, 02841 1.05195,1.06604,1.08224,1.10085,1.12229,1.14717,1.17623,1.20963, 02842 1.24296,1.26532,1.27619,1.28126,1.28362,1.28458,1.28477,1.28450, 02843 1.28394,1.28320,1.28235,1.28141,1.28042,1.27941,1.27837,1.27732, 02844 1.27627,1.27522,1.27418,1.27314,1.27210,1.27108,1.27006,1.26906, 02845 1.26806,1.26707,1.26610,1.26514,1.26418,1.26325,1.26233,1.26147, 02846 1.26147,1.26147,1.25992,1.25992, 0.0}; 02847 02848 i = idint(y/(2.0e-02)) + 1; 02849 assert(i >= 1); 02850 02851 if(i >= bsbkSize) { 02852 if(verboseLevel > 2) { 02853 G4cout <<"G4Abla error: index i = " << i << " is greater than array size permits." << G4endl; 02854 } 02855 bipolResult = 0.0; 02856 } 02857 else { 02858 if (iflag == 1) { 02859 bipolResult = bs[i] + (bs[i+1] - bs[i])/2.0e-02 * (y - 2.0e-02*(i - 1)); 02860 } 02861 else { 02862 bipolResult = bk[i] + (bk[i+1] - bk[i])/2.0e-02 * (y - 2.0e-02*(i - 1)); 02863 } 02864 } 02865 02866 // assert(isnan(bipolResult) == false); 02867 return bipolResult; 02868 }
| void G4Abla::breakItUp | ( | G4double | nucleusA, | |
| G4double | nucleusZ, | |||
| G4double | nucleusMass, | |||
| G4double | excitationEnergy, | |||
| G4double | angularMomentum, | |||
| G4double | recoilEnergy, | |||
| G4double | momX, | |||
| G4double | momY, | |||
| G4double | momZ, | |||
| G4int | eventnumber | |||
| ) |
Main interface to the de-excitation code.
| nucleusA | mass number of the nucleus | |
| nucleusZ | charge number of the nucleus | |
| nucleusMass | mass of the nucleus | |
| excitationEnergy | excitation energy of the nucleus | |
| angularMomentum | angular momentum of the nucleus (produced as output by INCL4) | |
| recoilEnergy | recoil energy of the nucleus | |
| momX | momentum x-component | |
| momY | momentum y-component | |
| momZ | momentum z-component | |
| eventnumber | number of the event |
Definition at line 109 of file G4Abla.cc.
References G4Volant::acv, evapora(), G4cout, G4endl, G4Volant::iv, mglms(), pace2(), G4Volant::pcv, translab(), G4Volant::xcv, G4Volant::ycv, G4Volant::zcv, and G4Volant::zpcv.
00112 { 00113 const G4double uma = 931.4942; 00114 const G4double melec = 0.511; 00115 const G4double fmp = 938.27231; 00116 const G4double fmn = 939.56563; 00117 00118 G4double alrem = 0.0, berem = 0.0, garem = 0.0; 00119 G4double R[4][4]; // Rotation matrix 00120 G4double csdir1[4]; 00121 G4double csdir2[4]; 00122 G4double csrem[4]; 00123 G4double pfis_rem[4]; 00124 G4double pf1_rem[4]; 00125 for(G4int init_i = 0; init_i < 4; init_i++) { 00126 csdir1[init_i] = 0.0; 00127 csdir2[init_i] = 0.0; 00128 csrem[init_i] = 0.0; 00129 pfis_rem[init_i] = 0.0; 00130 pf1_rem[init_i] = 0.0; 00131 for(G4int init_j = 0; init_j < 4; init_j++) { 00132 R[init_i][init_j] = 0.0; 00133 } 00134 } 00135 00136 G4double plab1 = 0.0, gam1 = 0.0, eta1 = 0.0; 00137 G4double plab2 = 0.0, gam2 = 0.0, eta2 = 0.0; 00138 00139 G4double sitet = 0.0; 00140 G4double stet1 = 0.0; 00141 G4double stet2 = 0.0; 00142 00143 G4int nbpevap = 0; 00144 G4int mempaw = 0, memiv = 0; 00145 00146 G4double e_evapo = 0.0; 00147 G4double el = 0.0; 00148 G4double fmcv = 0.0; 00149 00150 G4double aff1 = 0.0; 00151 G4double zff1 = 0.0; 00152 G4double eff1 = 0.0; 00153 G4double aff2 = 0.0; 00154 G4double zff2 = 0.0; 00155 G4double eff2 = 0.0; 00156 00157 G4double v1 = 0.0, v2 = 0.0; 00158 00159 G4double t2 = 0.0; 00160 G4double ctet1 = 0.0; 00161 G4double ctet2 = 0.0; 00162 G4double phi1 = 0.0; 00163 G4double phi2 = 0.0; 00164 G4double p2 = 0.0; 00165 G4double epf2_out = 0.0 ; 00166 G4int lma_pf1 = 0, lmi_pf1 = 0; 00167 G4int lma_pf2 = 0, lmi_pf2 = 0; 00168 G4int nopart = 0; 00169 00170 G4double cst = 0.0, sst = 0.0, csf = 0.0, ssf = 0.0; 00171 00172 G4double zf = 0.0, af = 0.0, mtota = 0.0, pleva = 0.0, pxeva = 0.0, pyeva = 0.0; 00173 G4int ff = 0; 00174 G4int inum = eventnumber; 00175 G4int inttype = 0; 00176 G4double esrem = excitationEnergy; 00177 00178 G4double aprf = nucleusA; 00179 G4double zprf = nucleusZ; 00180 G4double mcorem = nucleusMass; 00181 G4double ee = excitationEnergy; 00182 G4double jprf = angularMomentum; // actually root-mean-squared 00183 00184 G4double erecrem = recoilEnergy; 00185 G4double trem = 0.0; 00186 G4double pxrem = momX; 00187 G4double pyrem = momY; 00188 G4double pzrem = momZ; 00189 00190 G4double remmass = 0.0; 00191 00192 varntp->ntrack = 0; 00193 // volant->iv = 0; 00194 volant->iv = 1; 00195 00196 G4double pcorem = std::sqrt(erecrem*(erecrem +2.*938.2796*nucleusA)); 00197 // G4double pcorem = std::sqrt(std::pow(momX,2) + std::pow(momY,2) + std::pow(momZ,2)); 00198 // assert(isnan(pcorem) == false); 00199 if(esrem >= 1.0e-3) { 00200 evapora(zprf,aprf,ee,jprf, &zf, &af, &mtota, &pleva, &pxeva, &pyeva, &ff, &inttype, &inum); 00201 } 00202 else { 00203 ff = 0; 00204 zf = zprf; 00205 af = aprf; 00206 pxeva = pxrem; 00207 pyeva = pyrem; 00208 pleva = pzrem; 00209 } 00210 // assert(isnan(zf) == false); 00211 // assert(isnan(af) == false); 00212 // assert(isnan(ee) == false); 00213 00214 if (ff == 1) { 00215 // Fission: 00216 // variable ee: Energy of fissioning nucleus above the fission barrier. 00217 // Calcul des impulsions des particules evaporees (avant fission) 00218 // dans le systeme labo. 00219 00220 trem = double(erecrem); 00221 remmass = pace2(aprf,zprf) + aprf*uma - zprf*melec; // canonic 00222 // remmass = mcorem + double(esrem); // ok 00223 // remmass = mcorem; //cugnon 00224 varntp->kfis = 1; 00225 G4double gamrem = (remmass + trem)/remmass; 00226 G4double etrem = std::sqrt(trem*(trem + 2.0*remmass))/remmass; 00227 // assert(isnan(etrem) == false); 00228 00229 // This is not treated as accurately as for the non fission case for which 00230 // the remnant mass is computed to satisfy the energy conservation 00231 // of evaporated particles. But it is not bad and more canonical! 00232 remmass = pace2(aprf,zprf) + aprf*uma - zprf*melec+double(esrem); // !canonic 00233 // Essais avec la masse de KHS (9/2002): 00234 el = 0.0; 00235 mglms(aprf,zprf,0,&el); 00236 remmass = zprf*fmp + (aprf-zprf)*fmn + el + double(esrem); 00237 00238 gamrem = std::sqrt(std::pow(pcorem,2) + std::pow(remmass,2))/remmass; 00239 // assert(isnan(gamrem) == false); 00240 etrem = pcorem/remmass; 00241 // assert(isnan(etrem) == false); 00242 00243 alrem = pxrem/pcorem; 00244 // assert(isnan(alrem) == false); 00245 berem = pyrem/pcorem; 00246 // assert(isnan(berem) == false); 00247 garem = pzrem/pcorem; 00248 // assert(isnan(garem) == false); 00249 00250 csrem[0] = 0.0; // Should not be used. 00251 csrem[1] = alrem; 00252 csrem[2] = berem; 00253 csrem[3] = garem; 00254 00255 // C Pour Vérif Remnant = evapo(Pre fission) + Noyau_fissionant (système Remnant) 00256 G4double bil_e = 0.0; 00257 G4double bil_px = 0.0; 00258 G4double bil_py = 0.0; 00259 G4double bil_pz = 0.0; 00260 G4double masse = 0.0; 00261 00262 for(G4int iloc = 1; iloc <= volant->iv; iloc++) { //DO iloc=1,iv 00263 // assert(isnan(volant->zpcv[iloc]) == false); 00264 // assert(volant->acv[iloc] != 0); 00265 // assert(volant->zpcv[iloc] != 0); 00266 mglms(double(volant->acv[iloc]),double(volant->zpcv[iloc]),0,&el); 00267 // assert(isnan(el) == false); 00268 masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el; 00269 // assert(isnan(masse) == false); 00270 bil_e = bil_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2)); 00271 // assert(isnan(bil_e) == false); 00272 bil_px = bil_px + volant->pcv[iloc]*(volant->xcv[iloc]); 00273 bil_py = bil_py + volant->pcv[iloc]*(volant->ycv[iloc]); 00274 bil_pz = bil_pz + volant->pcv[iloc]*(volant->zcv[iloc]); 00275 } // enddo 00276 // C Ce bilan (impulsion nulle) est parfait. (Bil_Px=Bil_Px+PXEVA....) 00277 00278 G4int ndec = 1; 00279 00280 if(volant->iv != 0) { //then 00281 if(verboseLevel > 2) { 00282 G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl; 00283 G4cout <<"1st Translab: Adding indices from " << ndec << " to " << volant->iv << G4endl; 00284 } 00285 nopart = varntp->ntrack - 1; 00286 translab(gamrem,etrem,csrem,nopart,ndec); 00287 if(verboseLevel > 2) { 00288 G4cout <<"Translab complete!" << G4endl; 00289 G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl; 00290 } 00291 } 00292 nbpevap = volant->iv; // nombre de particules d'evaporation traitees 00293 00294 // C 00295 // C Now calculation of the fission fragment distribution including 00296 // C evaporation from the fragments. 00297 // C 00298 00299 // C Distribution of the fission fragments: 00300 00301 // void fissionDistri(G4double a,G4double z,G4double e, 00302 // G4double &a1,G4double &z1,G4double &e1,G4double &v1, 00303 // G4double &a2,G4double &z2,G4double &e2,G4double &v2); 00304 00305 fissionDistri(af,zf,ee,aff1,zff1,eff1,v1,aff2,zff2,eff2,v2); 00306 00307 // C verif des A et Z decimaux: 00308 G4int na_f = int(std::floor(af + 0.5)); 00309 G4int nz_f = int(std::floor(zf + 0.5)); 00310 varntp->izfis = nz_f; // copie dans le ntuple 00311 varntp->iafis = na_f; 00312 G4int na_pf1 = int(std::floor(aff1 + 0.5)); 00313 G4int nz_pf1 = int(std::floor(zff1 + 0.5)); 00314 G4int na_pf2 = int(std::floor(aff2 + 0.5)); 00315 G4int nz_pf2 = int(std::floor(zff2 + 0.5)); 00316 00317 if((na_f != (na_pf1+na_pf2)) || (nz_f != (nz_pf1+nz_pf2))) { 00318 if(verboseLevel > 2) { 00319 G4cout <<"problemes arrondis dans la fission " << G4endl; 00320 G4cout << "af,zf,aff1,zff1,aff2,zff2" << G4endl; 00321 G4cout << af <<" , " << zf <<" , " << aff1 <<" , " << zff1 <<" , " << aff2 <<" , " << zff2 << G4endl; 00322 G4cout << "a,z,a1,z1,a2,z2 integer" << G4endl; 00323 G4cout << na_f <<" , " << nz_f <<" , " << na_pf1 <<" , " << nz_pf1 <<" , " << na_pf2 <<" , " << nz_pf2 << G4endl; 00324 } 00325 } 00326 00327 // Calcul de l'impulsion des PF dans le syteme noyau de fission: 00328 G4int kboud = idnint(zf); 00329 G4int jboud = idnint(af-zf); 00330 //G4double ef = fb->efa[kboud][jboud]; // barriere de fission 00331 G4double ef = fb->efa[jboud][kboud]; // barriere de fission 00332 // assert(isnan(ef) == false); 00333 varntp->estfis = ee + ef; // copie dans le ntuple 00334 00335 // C MASSEF = pace2(AF,ZF) 00336 // C MASSEF = MASSEF + AF*UMA - ZF*MELEC + EE + EF 00337 // C MASSE1 = pace2(DBLE(AFF1),DBLE(ZFF1)) 00338 // C MASSE1 = MASSE1 + AFF1*UMA - ZFF1*MELEC + EFF1 00339 // C MASSE2 = pace2(DBLE(AFF2),DBLE(ZFF2)) 00340 // C MASSE2 = MASSE2 + AFF2*UMA - ZFF2*MELEC + EFF2 00341 // C WRITE(6,*) 'MASSEF,MASSE1,MASSE2',MASSEF,MASSE1,MASSE2 00342 // C MGLMS est la fonction de masse cohérente avec KHS evapo-fis. 00343 // C Attention aux parametres, ici 0=OPTSHP, NO microscopic correct. 00344 mglms(af,zf,0,&el); 00345 // assert(isnan(el) == false); 00346 G4double massef = zf*fmp + (af - zf)*fmn + el + ee + ef; 00347 // assert(isnan(massef) == false); 00348 mglms(double(aff1),double(zff1),0,&el); 00349 // assert(isnan(el) == false); 00350 G4double masse1 = zff1*fmp + (aff1-zff1)*fmn + el + eff1; 00351 // assert(isnan(masse1) == false); 00352 mglms(aff2,zff2,0,&el); 00353 // assert(isnan(el) == false); 00354 G4double masse2 = zff2*fmp + (aff2 - zff2)*fmn + el + eff2; 00355 // assert(isnan(masse2) == false); 00356 // C WRITE(6,*) 'MASSEF,MASSE1,MASSE2',MASSEF,MASSE1,MASSE2 00357 G4double b = massef - masse1 - masse2; 00358 if(b < 0.0) { //then 00359 b=0.0; 00360 if(verboseLevel > 2) { 00361 G4cout <<"anomalie dans la fission: " << G4endl; 00362 G4cout << inum<< " , " << af<< " , " <<zf<< " , " <<massef<< " , " <<aff1<< " , " <<zff1<< " , " <<masse1<< " , " <<aff2<< " , " <<zff2<< " , " << masse2 << G4endl; 00363 } 00364 } //endif 00365 G4double t1 = b*(b + 2.0*masse2)/(2.0*massef); 00366 // assert(isnan(t1) == false); 00367 G4double p1 = std::sqrt(t1*(t1 + 2.0*masse1)); 00368 // assert(isnan(p1) == false); 00369 00370 G4double rndm; 00371 standardRandom(&rndm, &(hazard->igraine[13])); 00372 ctet1 = 2.0*rndm - 1.0; 00373 standardRandom(&rndm,&(hazard->igraine[9])); 00374 phi1 = rndm*2.0*3.141592654; 00375 00376 // C ----Coefs de la transformation de Lorentz (noyau de fission -> Remnant) 00377 G4double peva = std::pow(pxeva,2) + std::pow(pyeva,2) + std::pow(pleva,2); 00378 G4double gamfis = std::sqrt(std::pow(massef,2) + peva)/massef; 00379 // assert(isnan(gamfis) == false); 00380 peva = std::sqrt(peva); 00381 // assert(isnan(peva) == false); 00382 G4double etfis = peva/massef; 00383 00384 G4double epf1_in = 0.0; 00385 G4double epf1_out = 0.0; 00386 00387 // C ----Matrice de rotation (noyau de fission -> Remnant) 00388 if(peva >= 1.0e-4) { 00389 sitet = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2))/peva; 00390 // assert(isnan(sitet) == false); 00391 } 00392 if(sitet > 1.0e-5) { //then 00393 G4double cstet = pleva/peva; 00394 G4double siphi = pyeva/(sitet*peva); 00395 G4double csphi = pxeva/(sitet*peva); 00396 00397 R[1][1] = cstet*csphi; 00398 R[1][2] = -siphi; 00399 R[1][3] = sitet*csphi; 00400 R[2][1] = cstet*siphi; 00401 R[2][2] = csphi; 00402 R[2][3] = sitet*siphi; 00403 R[3][1] = -sitet; 00404 R[3][2] = 0.0; 00405 R[3][3] = cstet; 00406 } 00407 else { 00408 R[1][1] = 1.0; 00409 R[1][2] = 0.0; 00410 R[1][3] = 0.0; 00411 R[2][1] = 0.0; 00412 R[2][2] = 1.0; 00413 R[2][3] = 0.0; 00414 R[3][1] = 0.0; 00415 R[3][2] = 0.0; 00416 R[3][3] = 1.0; 00417 } // endif 00418 // c test de verif: 00419 00420 if((zff1 <= 0.0) || (aff1 <= 0.0) || (aff1 < zff1)) { //then 00421 if(verboseLevel > 2) { 00422 G4cout <<"zf = " << zf <<" af = " << af <<"ee = " << ee <<"zff1 = " << zff1 <<"aff1 = " << aff1 << G4endl; 00423 } 00424 } 00425 else { 00426 // C ---------------------- PF1 will evaporate 00427 epf1_in = double(eff1); 00428 epf1_out = epf1_in; 00429 // void evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf, 00430 // G4double *zf_par, G4double *af_par, G4double *mtota_par, 00431 // G4double *pleva_par, G4double *pxeva_par, G4double *pyeva_par, 00432 // G4double *ff_par, G4int *inttype_par, G4int *inum_par); 00433 G4double zf1, af1, malpha1, ffpleva1, ffpxeva1, ffpyeva1; 00434 G4int ff1, ftype1; 00435 evapora(zff1, aff1, epf1_out, 0.0, &zf1, &af1, &malpha1, &ffpleva1, 00436 &ffpxeva1, &ffpyeva1, &ff1, &ftype1, &inum); 00437 // C On ajoute le fragment: 00438 // assert(af1 > 0); 00439 volant->iv = volant->iv + 1; 00440 // assert(af1 != 0); 00441 // assert(zf1 != 0); 00442 volant->acv[volant->iv] = af1; 00443 volant->zpcv[volant->iv] = zf1; 00444 if(verboseLevel > 2) { 00445 G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl; 00446 } 00447 peva = std::sqrt(std::pow(ffpxeva1,2) + std::pow(ffpyeva1,2) + std::pow(ffpleva1,2)); 00448 // assert(isnan(peva) == false); 00449 volant->pcv[volant->iv] = peva; 00450 if(peva > 0.001) { // then 00451 volant->xcv[volant->iv] = ffpxeva1/peva; 00452 volant->ycv[volant->iv] = ffpyeva1/peva; 00453 volant->zcv[volant->iv] = ffpleva1/peva; 00454 } 00455 else { 00456 volant->xcv[volant->iv] = 1.0; 00457 volant->ycv[volant->iv] = 0.0; 00458 volant->zcv[volant->iv] = 0.0; 00459 } // end if 00460 00461 // C Pour Vérif evapo de PF1 dans le systeme du Noyau Fissionant 00462 G4double bil1_e = 0.0; 00463 G4double bil1_px = 0.0; 00464 G4double bil1_py=0.0; 00465 G4double bil1_pz=0.0; 00466 for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv 00467 // for(G4int iloc = nbpevap + 1; iloc <= volant->iv + 1; iloc++) { //do iloc=nbpevap+1,iv 00468 mglms(volant->acv[iloc], volant->zpcv[iloc],0,&el); 00469 masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el; 00470 // assert(isnan(masse) == false); 00471 bil1_e = bil1_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2)); 00472 // assert(isnan(bil1_e) == false); 00473 bil1_px = bil1_px + volant->pcv[iloc]*(volant->xcv[iloc]); 00474 bil1_py = bil1_py + volant->pcv[iloc]*(volant->ycv[iloc]); 00475 bil1_pz = bil1_pz + volant->pcv[iloc]*(volant->zcv[iloc]); 00476 } // enddo 00477 00478 // Calcul des cosinus directeurs de PF1 dans le Remnant et calcul 00479 // des coefs pour la transformation de Lorentz Systeme PF --> Systeme Remnant 00480 translabpf(masse1,t1,p1,ctet1,phi1,gamfis,etfis,R,&plab1,&gam1,&eta1,csdir1); 00481 00482 // calcul des impulsions des particules evaporees dans le systeme Remnant: 00483 if(verboseLevel > 2) { 00484 G4cout <<"2nd Translab (pf1 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl; 00485 } 00486 nopart = varntp->ntrack - 1; 00487 translab(gam1,eta1,csdir1,nopart,nbpevap+1); 00488 if(verboseLevel > 2) { 00489 G4cout <<"After translab call... varntp->ntrack = " << varntp->ntrack << G4endl; 00490 } 00491 memiv = nbpevap + 1; // memoires pour la future transformation 00492 mempaw = nopart; // remnant->labo pour pf1 et pf2. 00493 lmi_pf1 = nopart + nbpevap + 1; // indices min et max dans /var_ntp/ 00494 lma_pf1 = nopart + volant->iv; // des particules issues de pf1 00495 nbpevap = volant->iv; // nombre de particules d'evaporation traitees 00496 } // end if 00497 // C --------------------- End of PF1 calculation 00498 00499 // c test de verif: 00500 if((zff2 <= 0.0) || (aff2 <= 0.0) || (aff2 <= zff2)) { //then 00501 if(verboseLevel > 2) { 00502 G4cout << zf << " " << af << " " << ee << " " << zff2 << " " << aff2 << G4endl; 00503 } 00504 } 00505 else { 00506 // C ---------------------- PF2 will evaporate 00507 G4double epf2_in = double(eff2); 00508 G4double epf2_out = epf2_in; 00509 // void evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf, 00510 // G4double *zf_par, G4double *af_par, G4double *mtota_par, 00511 // G4double *pleva_par, G4double *pxeva_par, G4double *pyeva_par, 00512 // G4double *ff_par, G4int *inttype_par, G4int *inum_par); 00513 G4double zf2, af2, malpha2, ffpleva2, ffpxeva2, ffpyeva2; 00514 G4int ff2, ftype2; 00515 evapora(zff2,aff2,epf2_out,0.0,&zf2,&af2,&malpha2,&ffpleva2, 00516 &ffpxeva2,&ffpyeva2,&ff2,&ftype2,&inum); 00517 // C On ajoute le fragment: 00518 volant->iv = volant->iv + 1; 00519 volant->acv[volant->iv] = af2; 00520 volant->zpcv[volant->iv] = zf2; 00521 if(verboseLevel > 2) { 00522 G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl; 00523 } 00524 peva = std::sqrt(std::pow(ffpxeva2,2) + std::pow(ffpyeva2,2) + std::pow(ffpleva2,2)); 00525 // assert(isnan(peva) == false); 00526 volant->pcv[volant->iv] = peva; 00527 // exit(0); 00528 if(peva > 0.001) { //then 00529 volant->xcv[volant->iv] = ffpxeva2/peva; 00530 volant->ycv[volant->iv] = ffpyeva2/peva; 00531 volant->zcv[volant->iv] = ffpleva2/peva; 00532 } 00533 else { 00534 volant->xcv[volant->iv] = 1.0; 00535 volant->ycv[volant->iv] = 0.0; 00536 volant->zcv[volant->iv] = 0.0; 00537 } //end if 00538 // C Pour Vérif evapo de PF1 dans le systeme du Noyau Fissionant 00539 G4double bil2_e = 0.0; 00540 G4double bil2_px = 0.0; 00541 G4double bil2_py = 0.0; 00542 G4double bil2_pz = 0.0; 00543 // for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv 00544 for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv 00545 mglms(volant->acv[iloc],volant->zpcv[iloc],0,&el); 00546 masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el; 00547 bil2_e = bil2_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2)); 00548 // assert(isnan(bil2_e) == false); 00549 bil2_px = bil2_px + volant->pcv[iloc]*(volant->xcv[iloc]); 00550 bil2_py = bil2_py + volant->pcv[iloc]*(volant->ycv[iloc]); 00551 bil2_pz = bil2_pz + volant->pcv[iloc]*(volant->zcv[iloc]); 00552 } //enddo 00553 00554 // C ----Calcul des cosinus directeurs de PF2 dans le Remnant et calcul 00555 // c des coefs pour la transformation de Lorentz Systeme PF --> Systeme Remnant 00556 G4double t2 = b - t1; 00557 // G4double ctet2 = -ctet1; 00558 ctet2 = -1.0*ctet1; 00559 assert(std::fabs(ctet2) <= 1.0); 00560 // assert(isnan(ctet2) == false); 00561 phi2 = dmod(phi1+3.141592654,6.283185308); 00562 // assert(isnan(phi2) == false); 00563 G4double p2 = std::sqrt(t2*(t2+2.0*masse2)); 00564 // assert(isnan(p2) == false); 00565 00566 // void translabpf(G4double masse1, G4double t1, G4double p1, G4double ctet1, 00567 // G4double phi1, G4double gamrem, G4double etrem, G4double R[][4], 00568 // G4double *plab1, G4double *gam1, G4double *eta1, G4double csdir[]); 00569 translabpf(masse2,t2,p2,ctet2,phi2,gamfis,etfis,R,&plab2,&gam2,&eta2,csdir2); 00570 // C 00571 // C calcul des impulsions des particules evaporees dans le systeme Remnant: 00572 // c 00573 if(verboseLevel > 2) { 00574 G4cout <<"3rd Translab (pf2 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl; 00575 } 00576 nopart = varntp->ntrack - 1; 00577 translab(gam2,eta2,csdir2,nopart,nbpevap+1); 00578 lmi_pf2 = nopart + nbpevap + 1; // indices min et max dans /var_ntp/ 00579 lma_pf2 = nopart + volant->iv; // des particules issues de pf2 00580 } // end if 00581 // C --------------------- End of PF2 calculation 00582 00583 // C Pour vérifications: calculs du noyau fissionant et des PF dans 00584 // C le systeme du remnant. 00585 for(G4int iloc = 1; iloc <= 3; iloc++) { // do iloc=1,3 00586 pfis_rem[iloc] = 0.0; 00587 } // enddo 00588 G4double efis_rem, pfis_trav[4]; 00589 lorab(gamfis,etfis,massef,pfis_rem,&efis_rem,pfis_trav); 00590 rotab(R,pfis_trav,pfis_rem); 00591 00592 stet1 = std::sqrt(1.0 - std::pow(ctet1,2)); 00593 // assert(isnan(stet1) == false); 00594 pf1_rem[1] = p1*stet1*std::cos(phi1); 00595 pf1_rem[2] = p1*stet1*std::sin(phi1); 00596 pf1_rem[3] = p1*ctet1; 00597 G4double e1_rem; 00598 lorab(gamfis,etfis,masse1+t1,pf1_rem,&e1_rem,pfis_trav); 00599 rotab(R,pfis_trav,pf1_rem); 00600 00601 stet2 = std::sqrt(1.0 - std::pow(ctet2,2)); 00602 assert(std::pow(ctet2,2) >= 0.0); 00603 assert(std::pow(ctet2,2) <= 1.0); 00604 // assert(isnan(stet2) == false); 00605 00606 G4double pf2_rem[4]; 00607 G4double e2_rem; 00608 pf2_rem[1] = p2*stet2*std::cos(phi2); 00609 pf2_rem[2] = p2*stet2*std::sin(phi2); 00610 pf2_rem[3] = p2*ctet2; 00611 lorab(gamfis,etfis,masse2+t2,pf2_rem,&e2_rem,pfis_trav); 00612 rotab(R,pfis_trav,pf2_rem); 00613 // C Verif 0: Remnant = evapo_pre_fission + Noyau Fissionant 00614 bil_e = remmass - efis_rem - bil_e; 00615 bil_px = bil_px + pfis_rem[1]; 00616 bil_py = bil_py + pfis_rem[2]; 00617 bil_pz = bil_pz + pfis_rem[3]; 00618 // C Verif 1: noyau fissionant = PF1 + PF2 dans le systeme remnant 00619 // G4double bilan_e = efis_rem - e1_rem - e2_rem; 00620 // G4double bilan_px = pfis_rem[1] - pf1_rem[1] - pf2_rem[1]; 00621 // G4double bilan_py = pfis_rem[2] - pf1_rem[2] - pf2_rem[2]; 00622 // G4double bilan_pz = pfis_rem[3] - pf1_rem[3] - pf2_rem[3]; 00623 // C Verif 2: PF1 et PF2 egaux a toutes leurs particules evaporees 00624 // C (Systeme remnant) 00625 if((lma_pf1-lmi_pf1) != 0) { //then 00626 G4double bil_e_pf1 = e1_rem - epf1_out; 00627 G4double bil_px_pf1 = pf1_rem[1]; 00628 G4double bil_py_pf1 = pf1_rem[2]; 00629 G4double bil_pz_pf1 = pf1_rem[3]; 00630 for(G4int ipf1 = lmi_pf1; ipf1 <= lma_pf1; ipf1++) { //do ipf1=lmi_pf1,lma_pf1 00631 bil_e_pf1 = bil_e_pf1 - (std::pow(varntp->plab[ipf1],2) + std::pow(varntp->enerj[ipf1],2))/(2.0*(varntp->enerj[ipf1])); 00632 cst = std::cos(varntp->tetlab[ipf1]/57.2957795); 00633 sst = std::sin(varntp->tetlab[ipf1]/57.2957795); 00634 csf = std::cos(varntp->philab[ipf1]/57.2957795); 00635 ssf = std::sin(varntp->philab[ipf1]/57.2957795); 00636 bil_px_pf1 = bil_px_pf1 - varntp->plab[ipf1]*sst*csf; 00637 bil_py_pf1 = bil_py_pf1 - varntp->plab[ipf1]*sst*ssf; 00638 bil_pz_pf1 = bil_pz_pf1 - varntp->plab[ipf1]*cst; 00639 } // enddo 00640 } //endif 00641 00642 if((lma_pf2-lmi_pf2) != 0) { //then 00643 G4double bil_e_pf2 = e2_rem - epf2_out; 00644 G4double bil_px_pf2 = pf2_rem[1]; 00645 G4double bil_py_pf2 = pf2_rem[2]; 00646 G4double bil_pz_pf2 = pf2_rem[3]; 00647 for(G4int ipf2 = lmi_pf2; ipf2 <= lma_pf2; ipf2++) { //do ipf2=lmi_pf2,lma_pf2 00648 bil_e_pf2 = bil_e_pf2 - (std::pow(varntp->plab[ipf2],2) + std::pow(varntp->enerj[ipf2],2))/(2.0*(varntp->enerj[ipf2])); 00649 G4double cst = std::cos(varntp->tetlab[ipf2]/57.2957795); 00650 G4double sst = std::sin(varntp->tetlab[ipf2]/57.2957795); 00651 G4double csf = std::cos(varntp->philab[ipf2]/57.2957795); 00652 G4double ssf = std::sin(varntp->philab[ipf2]/57.2957795); 00653 bil_px_pf2 = bil_px_pf2 - varntp->plab[ipf2]*sst*csf; 00654 bil_py_pf2 = bil_py_pf2 - varntp->plab[ipf2]*sst*ssf; 00655 bil_pz_pf2 = bil_pz_pf2 - varntp->plab[ipf2]*cst; 00656 } // enddo 00657 } //endif 00658 // C 00659 // C ---- Transformation systeme Remnant -> systeme labo. (evapo de PF1 ET PF2) 00660 // C 00661 // G4double mempaw, memiv; 00662 if(verboseLevel > 2) { 00663 G4cout <<"4th Translab: Adding indices from " << memiv << " to " << volant->iv << G4endl; 00664 } 00665 translab(gamrem,etrem,csrem,mempaw,memiv); 00666 // C ******************* END of fission calculations ************************ 00667 } 00668 else { 00669 // C ************************ Evapo sans fission ***************************** 00670 // C Here, FF=0, --> Evapo sans fission, on ajoute le fragment: 00671 // C ************************************************************************* 00672 varntp->kfis = 0; 00673 if(verboseLevel > 2) { 00674 G4cout <<"Evaporation without fission" << G4endl; 00675 } 00676 volant->iv = volant->iv + 1; 00677 volant->acv[volant->iv] = af; 00678 volant->zpcv[volant->iv] = zf; 00679 G4double peva = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2)+std::pow(pleva,2)); 00680 // assert(isnan(peva) == false); 00681 volant->pcv[volant->iv] = peva; 00682 if(peva > 0.001) { //then 00683 volant->xcv[volant->iv] = pxeva/peva; 00684 volant->ycv[volant->iv] = pyeva/peva; 00685 volant->zcv[volant->iv] = pleva/peva; 00686 } 00687 else { 00688 volant->xcv[volant->iv] = 1.0; 00689 volant->ycv[volant->iv] = 0.0; 00690 volant->zcv[volant->iv] = 0.0; 00691 } // end if 00692 00693 // C 00694 // C calcul des impulsions des particules evaporees dans le systeme labo: 00695 // c 00696 trem = double(erecrem); 00697 // C REMMASS = pace2(APRF,ZPRF) + APRF*UMA - ZPRF*MELEC !Canonic 00698 // C REMMASS = MCOREM + DBLE(ESREM) !OK 00699 remmass = mcorem; //Cugnon 00700 // C GAMREM = (REMMASS + TREM)/REMMASS !OK 00701 // C ETREM = DSQRT(TREM*(TREM + 2.*REMMASS))/REMMASS !OK 00702 csrem[0] = 0.0; // Should not be used. 00703 csrem[1] = alrem; 00704 csrem[2] = berem; 00705 csrem[3] = garem; 00706 00707 // for(G4int j = 1; j <= volant->iv; j++) { //do j=1,iv 00708 for(G4int j = 1; j <= volant->iv; j++) { //do j=1,iv 00709 if(volant->acv[j] == 0) { 00710 if(verboseLevel > 2) { 00711 G4cout <<"volant->acv[" << j << "] = 0" << G4endl; 00712 G4cout <<"volant->iv = " << volant->iv << G4endl; 00713 } 00714 } 00715 if(volant->acv[j] > 0) { 00716 assert(volant->acv[j] != 0); 00717 // assert(volant->zpcv[j] != 0); 00718 mglms(volant->acv[j],volant->zpcv[j],0,&el); 00719 fmcv = volant->zpcv[j]*fmp + (volant->acv[j] - volant->zpcv[j])*fmn + el; 00720 e_evapo = e_evapo + std::sqrt(std::pow(volant->pcv[j],2) + std::pow(fmcv,2)); 00721 // assert(isnan(e_evapo) == false); 00722 } 00723 } // enddo 00724 00725 // C Redefinition pour conservation d'impulsion!!! 00726 // C this mass obtained by energy balance is very close to the 00727 // C mass of the remnant computed by pace2 + excitation energy (EE). (OK) 00728 remmass = e_evapo; 00729 00730 G4double gamrem = std::sqrt(std::pow(pcorem,2)+std::pow(remmass,2))/remmass; 00731 // assert(isnan(gamrem) == false); 00732 G4double etrem = pcorem/remmass; 00733 00734 if(verboseLevel > 2) { 00735 G4cout <<"5th Translab (no fission): Adding indices from " << 1 << " to " << volant->iv << G4endl; 00736 } 00737 nopart = varntp->ntrack - 1; 00738 translab(gamrem,etrem,csrem,nopart,1); 00739 00740 // C End of the (FISSION - NO FISSION) condition (FF=1 or 0) 00741 } //end if 00742 // C *********************** FIN de l'EVAPO KHS ******************** 00743 }
KRAMERS FAKTOR - REDUCTION OF THE FISSION PROBABILITY INDEPENDENT OF EXCITATION ENERGY
Definition at line 2795 of file G4Abla.cc.
Referenced by direct().
02796 { 02797 // INPUT : BET, HOMEGA NUCLEAR VISCOSITY + CURVATURE OF POTENTIAL 02798 // OUTPUT: KRAMERS FAKTOR - REDUCTION OF THE FISSION PROBABILITY 02799 // INDEPENDENT OF EXCITATION ENERGY 02800 02801 G4double rel = bet/(20.0*homega/6.582122); 02802 G4double cramResult = std::sqrt(1.0 + std::pow(rel,2)) - rel; 02803 // limitation introduced 6.1.2000 by khs 02804 02805 if (cramResult > 1.0) { 02806 cramResult = 1.0; 02807 } 02808 02809 // assert(isnan(cramResult) == false); 02810 return cramResult; 02811 }
| void G4Abla::densniv | ( | G4double | a, | |
| G4double | z, | |||
| G4double | ee, | |||
| G4double | esous, | |||
| G4double * | dens, | |||
| G4double | bshell, | |||
| G4double | bs, | |||
| G4double | bk, | |||
| G4double * | temp, | |||
| G4int | optshp, | |||
| G4int | optcol, | |||
| G4double | defbet | |||
| ) |
Level density parameters.
Definition at line 2265 of file G4Abla.cc.
References G4Ald::ak, G4Ald::as, G4Ald::av, fe, idnint(), G4Ald::optafan, parite(), and qrot().
Referenced by direct().
02267 { 02268 // 1498 C 02269 // 1499 C INPUT: 02270 // 1500 C A,EE,ESOUS,OPTSHP,BS,BK,BSHELL,DEFBET 02271 // 1501 C 02272 // 1502 C LEVEL DENSITY PARAMETERS 02273 // 1503 C COMMON /ALD/ AV,AS,AK,OPTAFAN 02274 // 1504 C AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE 02275 // 1505 C LEVEL DENSITY PARAMETER 02276 // 1506 C OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1 02277 // 1507 C RECOMMENDED IS OPTAFAN = 0 02278 // 1508 C--------------------------------------------------------------------- 02279 // 1509 C OUTPUT: DENS,TEMP 02280 // 1510 C 02281 // 1511 C ____________________________________________________________________ 02282 // 1512 C / 02283 // 1513 C / PROCEDURE FOR CALCULATING THE STATE DENSITY OF A COMPOUND NUCLEUS 02284 // 1514 C /____________________________________________________________________ 02285 // 1515 C 02286 // 1516 INTEGER AFP,IZ,OPTSHP,OPTCOL,J,OPTAFAN 02287 // 1517 REAL*8 A,EE,ESOUS,DENS,E,Y0,Y1,Y2,Y01,Y11,Y21,PA,BS,BK,TEMP 02288 // 1518 C=====INSERTED BY KUDYAEV=============================================== 02289 // 1519 COMMON /ALD/ AV,AS,AK,OPTAFAN 02290 // 1520 REAL*8 ECR,ER,DELTAU,Z,DELTPP,PARA,PARZ,FE,HE,ECOR,ECOR1,Pi6 02291 // 1521 REAL*8 BSHELL,DELTA0,AV,AK,AS,PONNIV,PONFE,DEFBET,QR,SIG,FP 02292 // 1522 C======================================================================= 02293 // 1523 C 02294 // 1524 C 02295 // 1525 C----------------------------------------------------------------------- 02296 // 1526 C A MASS NUMBER OF THE DAUGHTER NUCLEUS 02297 // 1527 C EE EXCITATION ENERGY OF THE MOTHER NUCLEUS 02298 // 1528 C ESOUS SEPARATION ENERGY PLUS EFFECTIVE COULOMB BARRIER 02299 // 1529 C DENS STATE DENSITY OF DAUGHTER NUCLEUS AT EE-ESOUS-EC 02300 // 1530 C BSHELL SHELL CORRECTION 02301 // 1531 C TEMP NUCLEAR TEMPERATURE 02302 // 1532 C E LOCAL EXCITATION ENERGY OF THE DAUGHTER NUCLEUS 02303 // 1533 C E1 LOCAL HELP VARIABLE 02304 // 1534 C Y0,Y1,Y2,Y01,Y11,Y21 02305 // 1535 C LOCAL HELP VARIABLES 02306 // 1536 C PA LOCAL STATE-DENSITY PARAMETER 02307 // 1537 C EC KINETIC ENERGY OF EMITTED PARTICLE WITHOUT 02308 // 1538 C COULOMB REPULSION 02309 // 1539 C IDEN FAKTOR FOR SUBSTRACTING KINETIC ENERGY IDEN*TEMP 02310 // 1540 C DELTA0 PAIRING GAP 12 FOR GROUND STATE 02311 // 1541 C 14 FOR SADDLE POINT 02312 // 1542 C EITERA HELP VARIABLE FOR TEMPERATURE ITERATION 02313 // 1543 C----------------------------------------------------------------------- 02314 // 1544 C 02315 // 1545 C 02316 G4double afp = 0.0; 02317 G4double delta0 = 0.0; 02318 G4double deltau = 0.0; 02319 G4double deltpp = 0.0; 02320 G4double e = 0.0; 02321 G4double ecor = 0.0; 02322 G4double ecor1 = 0.0; 02323 G4double ecr = 0.0; 02324 G4double er = 0.0; 02325 G4double fe = 0.0; 02326 G4double fp = 0.0; 02327 G4double he = 0.0; 02328 G4double iz = 0.0; 02329 G4double pa = 0.0; 02330 G4double para = 0.0; 02331 G4double parz = 0.0; 02332 G4double ponfe = 0.0; 02333 G4double ponniv = 0.0; 02334 G4double qr = 0.0; 02335 G4double sig = 0.0; 02336 G4double y01 = 0.0; 02337 G4double y11 = 0.0; 02338 G4double y2 = 0.0; 02339 G4double y21 = 0.0; 02340 G4double y1 = 0.0; 02341 G4double y0 = 0.0; 02342 02343 G4double pi6 = std::pow(3.1415926535,2) / 6.0; 02344 ecr=10.0; 02345 er=28.0; 02346 afp=idnint(a); 02347 iz=idnint(z); 02348 02349 // level density parameter 02350 if((ald->optafan == 1)) { 02351 pa = (ald->av)*a + (ald->as)*std::pow(a,(2.e0/3.e0)) + (ald->ak)*std::pow(a,(1.e0/3.e0)); 02352 } 02353 else { 02354 pa = (ald->av)*a + (ald->as)*bs*std::pow(a,(2.e0/3.e0)) + (ald->ak)*bk*std::pow(a,(1.e0/3.e0)); 02355 } 02356 02357 fp = 0.01377937231e0 * std::pow(a,(5.e0/3.e0)) * (1.e0 + defbet/3.e0); 02358 02359 // pairing corrections 02360 if (bs > 1.0) { 02361 delta0 = 14; 02362 } 02363 else { 02364 delta0 = 12; 02365 } 02366 02367 if (esous > 1.0e30) { 02368 (*dens) = 0.0; 02369 (*temp) = 0.0; 02370 goto densniv100; 02371 } 02372 02373 e = ee - esous; 02374 02375 if (e < 0.0) { 02376 (*dens) = 0.0; 02377 (*temp) = 0.0; 02378 goto densniv100; 02379 } 02380 02381 // shell corrections 02382 if (optshp > 0) { 02383 deltau = bshell; 02384 if (optshp == 2) { 02385 deltau = 0.0; 02386 } 02387 if (optshp >= 2) { 02388 // pairing energy shift with condensation energy a.r.j. 10.03.97 02389 // deltpp = -0.25e0* (delta0/std::pow(std::sqrt(a),2)) * pa /pi6 + 2.e0*delta0/std::sqrt(a); 02390 deltpp = -0.25e0* std::pow((delta0/std::sqrt(a)),2) * pa /pi6 + 2.e0*delta0/std::sqrt(a); 02391 // assert(isnan(deltpp) == false); 02392 02393 parite(a,¶); 02394 if (para < 0.0) { 02395 e = e - delta0/std::sqrt(a); 02396 // assert(isnan(e) == false); 02397 } else { 02398 parite(z,&parz); 02399 if (parz > 0.e0) { 02400 e = e - 2.0*delta0/std::sqrt(a); 02401 // assert(isnan(e) == false); 02402 } else { 02403 e = e; 02404 // assert(isnan(e) == false); 02405 } 02406 } 02407 } else { 02408 deltpp = 0.0; 02409 } 02410 } else { 02411 deltau = 0.0; 02412 deltpp = 0.0; 02413 } 02414 if (e < 0.0) { 02415 e = 0.0; 02416 (*temp) = 0.0; 02417 } 02418 02419 // washing out is made stronger ! g.k. 3.7.96 02420 ponfe = -2.5*pa*e*std::pow(a,(-4.0/3.0)); 02421 02422 if (ponfe < -700.0) { 02423 ponfe = -700.0; 02424 } 02425 fe = 1.0 - std::exp(ponfe); 02426 if (e < ecr) { 02427 // priv. comm. k.-h. schmidt 02428 he = 1.0 - std::pow((1.0 - e/ecr),2); 02429 } 02430 else { 02431 he = 1.0; 02432 } 02433 02434 // Excitation energy corrected for pairing and shell effects 02435 // washing out with excitation energy is included. 02436 ecor = e + deltau*fe + deltpp*he; 02437 02438 if (ecor <= 0.1) { 02439 ecor = 0.1; 02440 } 02441 02442 // statt 170.d0 a.r.j. 8.11.97 02443 02444 // iterative procedure according to grossjean and feldmeier 02445 // to avoid the singularity e = 0 02446 if (ee < 5.0) { 02447 y1 = std::sqrt(pa*ecor); 02448 // assert(isnan(y1) == false); 02449 for(int j = 0; j < 5; j++) { 02450 y2 = pa*ecor*(1.e0-std::exp(-y1)); 02451 // assert(isnan(y2) == false); 02452 y1 = std::sqrt(y2); 02453 // assert(isnan(y1) == false); 02454 } 02455 02456 y0 = pa/y1; 02457 // assert(isnan(y0) == false); 02458 assert(y0 != 0.0); 02459 (*temp)=1.0/y0; 02460 (*dens) = std::exp(y0*ecor)/ (std::pow((std::pow(ecor,3)*y0),0.5)*std::pow((1.0-0.5*y0*ecor*std::exp(-y1)),0.5))* std::exp(y1)*(1.0-std::exp(-y1))*0.1477045; 02461 if (ecor < 1.0) { 02462 ecor1=1.0; 02463 y11 = std::sqrt(pa*ecor1); 02464 // assert(isnan(y11) == false); 02465 for(int j = 0; j < 7; j++) { 02466 y21 = pa*ecor1*(1.0-std::exp(-y11)); 02467 // assert(isnan(21) == false); 02468 y11 = std::sqrt(y21); 02469 // assert(isnan(y11) == false); 02470 } 02471 02472 y01 = pa/y11; 02473 // assert(isnan(y01) == false); 02474 (*dens) = (*dens)*std::pow((y01/y0),1.5); 02475 (*temp) = (*temp)*std::pow((y01/y0),1.5); 02476 } 02477 } 02478 else { 02479 ponniv = 2.0*std::sqrt(pa*ecor); 02480 // assert(isnan(ponniv) == false); 02481 if (ponniv > 700.0) { 02482 ponniv = 700.0; 02483 } 02484 02485 // fermi gas state density 02486 (*dens) = std::pow(pa,(-0.25e0))*std::pow(ecor,(-1.25e0))*std::exp(ponniv) * 0.1477045e0; 02487 // assert(isnan(std::sqrt(ecor/pa)) == false); 02488 (*temp) = std::sqrt(ecor/pa); 02489 } 02490 densniv100: 02491 02492 // spin cutoff parameter 02493 sig = fp * (*temp); 02494 02495 // collective enhancement 02496 if (optcol == 1) { 02497 qrot(z,a,defbet,sig,ecor,&qr); 02498 } 02499 else { 02500 qr = 1.0; 02501 } 02502 02503 (*dens) = (*dens) * qr; 02504 }
| void G4Abla::direct | ( | G4double | zprf, | |
| G4double | a, | |||
| G4double | ee, | |||
| G4double | jprf, | |||
| G4double * | probp_par, | |||
| G4double * | probn_par, | |||
| G4double * | proba_par, | |||
| G4double * | probf_par, | |||
| G4double * | ptotl_par, | |||
| G4double * | sn_par, | |||
| G4double * | sbp_par, | |||
| G4double * | sba_par, | |||
| G4double * | ecn_par, | |||
| G4double * | ecp_par, | |||
| G4double * | eca_par, | |||
| G4double * | bp_par, | |||
| G4double * | ba_par, | |||
| G4int | inttype, | |||
| G4int | inum, | |||
| G4int | itest | |||
| ) |
Calculation of particle emission probabilities.
Definition at line 1533 of file G4Abla.cc.
References G4Fiss::akap, G4Ecld::alpha, barfit(), G4Fiss::bet, bipol(), cram(), densniv(), G4Ecld::ecfnz, G4Ecld::ecgnz, G4Fb::efa, fissility(), fmaxhaz(), G4cout, G4endl, G4Fiss::homega, idnint(), G4Fiss::ifis, mglms(), mglw(), G4Fiss::optcol, G4Fiss::optles, G4Fiss::optshp, G4Fiss::optxfis, parite(), G4INCL::Math::pi, G4InuclParticleNames::sp, spdef(), standardRandom(), tau(), and G4Ecld::vgsld.
Referenced by evapora().
01539 { 01540 G4int dummy0 = 0; 01541 01542 G4double probp = (*probp_par); 01543 G4double probn = (*probn_par); 01544 G4double proba = (*proba_par); 01545 G4double probf = (*probf_par); 01546 G4double ptotl = (*ptotl_par); 01547 G4double sn = (*sn_par); 01548 G4double sbp = (*sbp_par); 01549 G4double sba = (*sba_par); 01550 G4double ecn = (*ecn_par); 01551 G4double ecp = (*ecp_par); 01552 G4double eca = (*eca_par); 01553 G4double bp = (*bp_par); 01554 G4double ba = (*ba_par); 01555 01556 // CALCULATION OF PARTICLE-EMISSION PROBABILITIES & FISSION / 01557 // BASED ON THE SIMPLIFIED FORMULAS FOR THE DECAY WIDTH BY / 01558 // MORETTO, ROCHESTER MEETING TO AVOID COMPUTING TIME / 01559 // INTENSIVE INTEGRATION OF THE LEVEL DENSITIES / 01560 // USES EFFECTIVE COULOMB BARRIERS AND AN AVERAGE KINETIC ENERGY/ 01561 // OF THE EVAPORATED PARTICLES / 01562 // COLLECTIVE ENHANCMENT OF THE LEVEL DENSITY IS INCLUDED / 01563 // DYNAMICAL HINDRANCE OF FISSION IS INCLUDED BY A STEP FUNCTION/ 01564 // APPROXIMATION. SEE A.R. JUNGHANS DIPLOMA THESIS / 01565 // SHELL AND PAIRING STRUCTURES IN THE LEVEL DENSITY IS INCLUDED/ 01566 01567 // INPUT: 01568 // ZPRF,A,EE CHARGE, MASS, EXCITATION ENERGY OF COMPOUND 01569 // NUCLEUS 01570 // JPRF ROOT-MEAN-SQUARED ANGULAR MOMENTUM 01571 01572 // DEFORMATIONS AND G.S. SHELL EFFECTS 01573 // COMMON /ECLD/ ECGNZ,ECFNZ,VGSLD,ALPHA 01574 01575 // ECGNZ - GROUND STATE SHELL CORR. FRLDM FOR A SPHERICAL G.S. 01576 // ECFNZ - SHELL CORRECTION FOR THE SADDLE POINT (NOW: == 0) 01577 // VGSLD - DIFFERENCE BETWEEN DEFORMED G.S. AND LDM VALUE 01578 // ALPHA - ALPHA GROUND STATE DEFORMATION (THIS IS NOT BETA2!) 01579 // BETA2 = SQRT((4PI)/5) * ALPHA 01580 01581 //OPTIONS AND PARAMETERS FOR FISSION CHANNEL 01582 //COMMON /FISS/ AKAP,BET,HOMEGA,KOEFF,IFIS, 01583 // OPTSHP,OPTXFIS,OPTLES,OPTCOL 01584 // 01585 // AKAP - HBAR**2/(2* MN * R_0**2) = 10 MEV, R_0 = 1.4 FM 01586 // BET - REDUCED NUCLEAR FRICTION COEFFICIENT IN (10**21 S**-1) 01587 // HOMEGA - CURVATURE OF THE FISSION BARRIER = 1 MEV 01588 // KOEFF - COEFFICIENT FOR THE LD FISSION BARRIER == 1.0 01589 // IFIS - 0/1 FISSION CHANNEL OFF/ON 01590 // OPTSHP - INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY 01591 // = 0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY 01592 // = 1 SHELL , NO PAIRING 01593 // = 2 PAIRING, NO SHELL 01594 // = 3 SHELL AND PAIRING 01595 // OPTCOL - 0/1 COLLECTIVE ENHANCEMENT SWITCHED ON/OFF 01596 // OPTXFIS- 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV 01597 // FISSILITY PARAMETER. 01598 // OPTLES - CONSTANT TEMPERATURE LEVEL DENSITY FOR A,Z > TH-224 01599 // OPTCOL - 0/1 COLLECTIVE ENHANCEMENT OFF/ON 01600 01601 // LEVEL DENSITY PARAMETERS 01602 // COMMON /ALD/ AV,AS,AK,OPTAFAN 01603 // AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE 01604 // LEVEL DENSITY PARAMETER 01605 // OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1 01606 // RECOMMENDED IS OPTAFAN = 0 01607 01608 // FISSION BARRIERS 01609 // COMMON /FB/ EFA 01610 // EFA - ARRAY OF FISSION BARRIERS 01611 01612 01613 // OUTPUT: PROBN,PROBP,PROBA,PROBF,PTOTL: 01614 // - EMISSION PROBABILITIES FOR N EUTRON, P ROTON, A LPHA 01615 // PARTICLES, F ISSION AND NORMALISATION 01616 // SN,SBP,SBA: SEPARATION ENERGIES N P A 01617 // INCLUDING EFFECTIVE BARRIERS 01618 // ECN,ECP,ECA,BP,BA 01619 // - AVERAGE KINETIC ENERGIES (2*T) AND EFFECTIVE BARRIERS 01620 01621 static G4double bk = 0.0; 01622 static G4int afp = 0; 01623 static G4double at = 0.0; 01624 static G4double bs = 0.0; 01625 static G4double bshell = 0.0; 01626 static G4double cf = 0.0; 01627 static G4double dconst = 0.0; 01628 static G4double defbet = 0.0; 01629 static G4double denomi = 0.0; 01630 static G4double densa = 0.0; 01631 static G4double densf = 0.0; 01632 static G4double densg = 0.0; 01633 static G4double densn = 0.0; 01634 static G4double densp = 0.0; 01635 static G4double edyn = 0.0; 01636 static G4double eer = 0.0; 01637 static G4double ef = 0.0; 01638 static G4double ft = 0.0; 01639 static G4double ga = 0.0; 01640 static G4double gf = 0.0; 01641 static G4double gn = 0.0; 01642 static G4double gngf = 0.0; 01643 static G4double gp = 0.0; 01644 static G4double gsum = 0.0; 01645 static G4double hbar = 6.582122e-22; // = 0.0; 01646 static G4double iflag = 0.0; 01647 static G4int il = 0; 01648 static G4int imaxwell = 0; 01649 static G4int in = 0; 01650 static G4int iz = 0; 01651 static G4int j = 0; 01652 static G4int k = 0; 01653 static G4double ma1z = 0.0; 01654 static G4double ma1z1 = 0.0; 01655 static G4double ma4z2 = 0.0; 01656 static G4double maz = 0.0; 01657 static G4double nprf = 0.0; 01658 static G4double nt = 0.0; 01659 static G4double parc = 0.0; 01660 static G4double pi = 3.14159265; 01661 static G4double pt = 0.0; 01662 static G4double ra = 0.0; 01663 static G4double rat = 0.0; 01664 static G4double refmod = 0.0; 01665 static G4double rf = 0.0; 01666 static G4double rn = 0.0; 01667 static G4double rnd = 0.0; 01668 static G4double rnt = 0.0; 01669 static G4double rp = 0.0; 01670 static G4double rpt = 0.0; 01671 static G4double sa = 0.0; 01672 static G4double sbf = 0.0; 01673 static G4double sbfis = 0.0; 01674 static G4double segs = 0.0; 01675 static G4double selmax = 0.0; 01676 static G4double sp = 0.0; 01677 static G4double tauc = 0.0; 01678 static G4double tconst = 0.0; 01679 static G4double temp = 0.0; 01680 static G4double ts1 = 0.0; 01681 static G4double tsum = 0.0; 01682 static G4double wf = 0.0; 01683 static G4double wfex = 0.0; 01684 static G4double xx = 0.0; 01685 static G4double y = 0.0; 01686 01687 imaxwell = 1; 01688 inttype = 0; 01689 01690 // limiting of excitation energy where fission occurs 01691 // Note, this is not the dynamical hindrance (see end of routine) 01692 edyn = 1000.0; 01693 01694 // no limit if statistical model is calculated. 01695 if (fiss->bet <= 1.0e-16) { 01696 edyn = 10000.0; 01697 } 01698 01699 // just a change of name until the end of this subroutine 01700 eer = ee; 01701 if (inum == 1) { 01702 ilast = 1; 01703 } 01704 01705 // calculation of masses 01706 // refmod = 1 ==> myers,swiatecki model 01707 // refmod = 0 ==> weizsaecker model 01708 refmod = 1; // Default = 1 01709 01710 if (refmod == 1) { 01711 mglms(a,zprf,fiss->optshp,&maz); 01712 mglms(a-1.0,zprf,fiss->optshp,&ma1z); 01713 mglms(a-1.0,zprf-1.0,fiss->optshp,&ma1z1); 01714 mglms(a-4.0,zprf-2.0,fiss->optshp,&ma4z2); 01715 } 01716 else { 01717 mglw(a,zprf,&maz); 01718 mglw(a-1.0,zprf,&ma1z); 01719 mglw(a-1.0,zprf-1.0,&ma1z1); 01720 mglw(a-4.0,zprf-2.0,&ma4z2); 01721 } 01722 // assert(isnan(maz) == false); 01723 // assert(isnan(ma1z) == false); 01724 // assert(isnan(ma1z1) == false); 01725 // assert(isnan(ma4z2) == false); 01726 01727 // separation energies and effective barriers 01728 sn = ma1z - maz; 01729 sp = ma1z1 - maz; 01730 sa = ma4z2 - maz - 28.29688; 01731 if (zprf < 1.0e0) { 01732 sbp = 1.0e75; 01733 goto direct30; 01734 } 01735 01736 // parameterisation gaimard: 01737 // bp = 1.44*(zprf-1.d0)/(1.22*std::pow((a - 1.0),(1.0/3.0))+5.6) 01738 // parameterisation khs (12-99) 01739 bp = 1.44*(zprf - 1.0)/(2.1*std::pow((a - 1.0),(1.0/3.0)) + 0.0); 01740 01741 sbp = sp + bp; 01742 // assert(isnan(sbp) == false); 01743 // assert(isinf(sbp) == false); 01744 if (a-4.0 <= 0.0) { 01745 sba = 1.0e+75; 01746 goto direct30; 01747 } 01748 01749 // new effective barrier for alpha evaporation d=6.1: khs 01750 // ba = 2.88d0*(zprf-2.d0)/(1.22d0*(a-4.d0)**(1.d0/3.d0)+6.1d0) 01751 // parametrisation khs (12-99) 01752 ba = 2.88*(zprf - 2.0)/(2.2*std::pow((a - 4.0),(1.0/3.0)) + 0.0); 01753 01754 sba = sa + ba; 01755 // assert(isnan(sba) == false); 01756 // assert(isinf(sba) == false); 01757 direct30: 01758 01759 // calculation of surface and curvature integrals needed to 01760 // to calculate the level density parameter (in densniv) 01761 if (fiss->ifis > 0) { 01762 k = idnint(zprf); 01763 j = idnint(a - zprf); 01764 01765 // now ef is calculated from efa that depends on the subroutine 01766 // barfit which takes into account the modification on the ang. mom. 01767 // jb mvr 6-aug-1999 01768 // note *** shell correction! (ecgnz) jb mvr 20-7-1999 01769 il = idnint(jprf); 01770 barfit(k,k+j,il,&sbfis,&segs,&selmax); 01771 // assert(isnan(sbfis) == false); 01772 if ((fiss->optshp == 1) || (fiss->optshp == 3)) { 01773 // fb->efa[k][j] = G4double(sbfis) - ecld->ecgnz[j][k]; 01774 // fb->efa[j][k] = G4double(sbfis) - ecld->ecgnz[j][k]; 01775 fb->efa[j][k] = double(sbfis) - ecld->ecgnz[j][k]; 01776 } 01777 else { 01778 // fb->efa[k][j] = G4double(sbfis); 01779 fb->efa[j][k] = double(sbfis); 01780 } 01781 // ef = fb->efa[k][j]; 01782 ef = fb->efa[j][k]; 01783 01784 // to avoid negative values for impossible nuclei 01785 // the fission barrier is set to zero if smaller than zero. 01786 // if (fb->efa[k][j] < 0.0) { 01787 // fb->efa[k][j] = 0.0; 01788 // } 01789 if (fb->efa[j][k] < 0.0) { 01790 if(verboseLevel > 2) { 01791 G4cout <<"Setting fission barrier to 0" << G4endl; 01792 } 01793 fb->efa[j][k] = 0.0; 01794 } 01795 // assert(isnan(fb->efa[j][k]) == false); 01796 01797 // factor with jprf should be 0.0025d0 - 0.01d0 for 01798 // approximate influence of ang. momentum on bfis a.j. 22.07.96 01799 // 0.0 means no angular momentum 01800 01801 if (ef < 0.0) { 01802 ef = 0.0; 01803 } 01804 xx = fissility((k+j),k,fiss->optxfis); 01805 // assert(isnan(xx) == false); 01806 // assert(isinf(xx) == false); 01807 01808 y = 1.00 - xx; 01809 if (y < 0.0) { 01810 y = 0.0; 01811 } 01812 if (y > 1.0) { 01813 y = 1.0; 01814 } 01815 bs = bipol(1,y); 01816 // assert(isnan(bs) == false); 01817 // assert(isinf(bs) == false); 01818 bk = bipol(2,y); 01819 // assert(isnan(bk) == false); 01820 // assert(isinf(bk) == false); 01821 } 01822 else { 01823 ef = 1.0e40; 01824 bs = 1.0; 01825 bk = 1.0; 01826 } 01827 sbf = ee - ef; 01828 01829 afp = idnint(a); 01830 iz = idnint(zprf); 01831 in = afp - iz; 01832 bshell = ecld->ecfnz[in][iz]; 01833 // assert(isnan(bshell) == false); 01834 01835 // ld saddle point deformation 01836 // here: beta2 = std::sqrt(5/(4pi)) * alpha2 01837 01838 // for the ground state def. 1.5d0 should be used 01839 // because this was just the factor to produce the 01840 // alpha-deformation table 1.5d0 should be used 01841 // a.r.j. 6.8.97 01842 defbet = 1.58533e0 * spdef(idnint(a),idnint(zprf),fiss->optxfis); 01843 // assert(isnan(defbet) == false); 01844 01845 // level density and temperature at the saddle point 01846 // G4cout <<"a = " << a << G4endl; 01847 // G4cout <<"zprf = " << zprf << G4endl; 01848 // G4cout <<"ee = " << ee << G4endl; 01849 // G4cout <<"ef = " << ef << G4endl; 01850 // G4cout <<"bshell = " << bshell << G4endl; 01851 // G4cout <<"bs = " << bs << G4endl; 01852 // G4cout <<"bk = " << bk << G4endl; 01853 // G4cout <<"defbet = " << defbet << G4endl; 01854 densniv(a,zprf,ee,ef,&densf,bshell,bs,bk,&temp,int(fiss->optshp),int(fiss->optcol),defbet); 01855 // G4cout <<"densf = " << densf << G4endl; 01856 // G4cout <<"temp = " << temp << G4endl; 01857 // assert(isnan(densf) == false); 01858 // assert(isnan(temp) == false); 01859 // assert(temp != 0); 01860 ft = temp; 01861 if (iz >= 2) { 01862 bshell = ecld->ecgnz[in][iz-1] - ecld->vgsld[in][iz-1]; 01863 defbet = 1.5 * (ecld->alpha[in][iz-1]); 01864 01865 // level density and temperature in the proton daughter 01866 densniv(a-1.0,zprf-1.0e0,ee,sbp,&densp, bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet); 01867 assert(temp >= 0); 01868 // assert(isnan(temp) == false); 01869 pt = temp; 01870 if (imaxwell == 1) { 01871 // valentina - random kinetic energy in a maxwelliam distribution 01872 // modif juin/2002 a.b. c.v. for light targets; limit on the energy 01873 // from the maxwell distribution. 01874 rpt = pt; 01875 ecp = 2.0 * pt; 01876 if(rpt <= 1.0e-3) { 01877 goto direct2914; 01878 } 01879 iflag = 0; 01880 direct1914: 01881 ecp = fmaxhaz(rpt); 01882 iflag = iflag + 1; 01883 if(iflag >= 10) { 01884 standardRandom(&rnd,&(hazard->igraine[5])); 01885 ecp=std::sqrt(rnd)*(eer-sbp); 01886 // assert(isnan(ecp) == false); 01887 goto direct2914; 01888 } 01889 if((ecp+sbp) > eer) { 01890 goto direct1914; 01891 } 01892 } 01893 else { 01894 ecp = 2.0 * pt; 01895 } 01896 01897 direct2914: 01898 dummy0 = 0; 01899 // G4cout <<""<<G4endl; 01900 } 01901 else { 01902 densp = 0.0; 01903 ecp = 0.0; 01904 pt = 0.0; 01905 } 01906 01907 if (in >= 2) { 01908 bshell = ecld->ecgnz[in-1][iz] - ecld->vgsld[in-1][iz]; 01909 defbet = 1.5e0 * (ecld->alpha[in-1][iz]); 01910 01911 // level density and temperature in the neutron daughter 01912 densniv(a-1.0,zprf,ee,sn,&densn,bshell, 1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet); 01913 nt = temp; 01914 01915 if (imaxwell == 1) { 01916 // valentina - random kinetic energy in a maxwelliam distribution 01917 // modif juin/2002 a.b. c.v. for light targets; limit on the energy 01918 // from the maxwell distribution. 01919 rnt = nt; 01920 ecn = 2.0 * nt; 01921 if(rnt <= 1.e-3) { 01922 goto direct2915; 01923 } 01924 01925 iflag=0; 01926 01927 ecn = fmaxhaz(rnt); 01928 iflag=iflag+1; 01929 if(iflag >= 10) { 01930 standardRandom(&rnd,&(hazard->igraine[6])); 01931 ecn = std::sqrt(rnd)*(eer-sn); 01932 // assert(isnan(ecn) == false); 01933 goto direct2915; 01934 } 01935 // if((ecn+sn) > eer) { 01936 // goto direct1915; 01937 // } 01938 // else { 01939 // ecn = 2.e0 * nt; 01940 // } 01941 if((ecn + sn) <= eer) { 01942 ecn = 2.0 * nt; 01943 } 01944 direct2915: 01945 dummy0 = 0; 01946 // G4cout <<"" <<G4endl; 01947 } 01948 } 01949 else { 01950 densn = 0.0; 01951 ecn = 0.0; 01952 nt = 0.0; 01953 } 01954 01955 if ((in >= 3) && (iz >= 3)) { 01956 bshell = ecld->ecgnz[in-2][iz-2] - ecld->vgsld[in-2][iz-2]; 01957 defbet = 1.5 * (ecld->alpha[in-2][iz-2]); 01958 01959 // level density and temperature in the alpha daughter 01960 densniv(a-4.0,zprf-2.0e0,ee,sba,&densa,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet); 01961 01962 // valentina - random kinetic energy in a maxwelliam distribution 01963 at = temp; 01964 if (imaxwell == 1) { 01965 // modif juin/2002 a.b. c.v. for light targets; limit on the energy 01966 // from the maxwell distribution. 01967 rat = at; 01968 eca= 2.e0 * at; 01969 if(rat <= 1.e-3) { 01970 goto direct2916; 01971 } 01972 iflag=0; 01973 direct1916: 01974 eca = fmaxhaz(rat); 01975 iflag=iflag+1; 01976 if(iflag >= 10) { 01977 standardRandom(&rnd,&(hazard->igraine[7])); 01978 eca=std::sqrt(rnd)*(eer-sba); 01979 // assert(isnan(eca) == false); 01980 goto direct2916; 01981 } 01982 if((eca+sba) > eer) { 01983 goto direct1916; 01984 } 01985 else { 01986 eca = 2.0 * at; 01987 } 01988 direct2916: 01989 dummy0 = 0; 01990 // G4cout <<"" << G4endl; 01991 } 01992 else { 01993 densa = 0.0; 01994 eca = 0.0; 01995 at = 0.0; 01996 } 01997 } // PK 01998 01999 // special treatment for unbound nuclei 02000 if (sn < 0.0) { 02001 probn = 1.0; 02002 probp = 0.0; 02003 proba = 0.0; 02004 probf = 0.0; 02005 goto direct70; 02006 } 02007 if (sbp < 0.0) { 02008 probp = 1.0; 02009 probn = 0.0; 02010 proba = 0.0; 02011 probf = 0.0; 02012 goto direct70; 02013 } 02014 02015 if ((a < 50.e0) || (ee > edyn)) { // no fission if e*> edyn or mass < 50 02016 // G4cout <<"densf = 0.0" << G4endl; 02017 densf = 0.e0; 02018 } 02019 02020 bshell = ecld->ecgnz[in][iz] - ecld->vgsld[in][iz]; 02021 defbet = 1.5e0 * (ecld->alpha[in][iz]); 02022 02023 // compound nucleus level density 02024 densniv(a,zprf,ee,0.0e0,&densg,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet); 02025 // assert(isnan(densg) == false); 02026 // assert(isnan(temp) == false); 02027 02028 if ( densg > 0.e0) { 02029 // calculation of the partial decay width 02030 // used for both the time scale and the evaporation decay width 02031 gp = (std::pow(a,(2.0/3.0))/fiss->akap)*densp/densg/pi*std::pow(pt,2); 02032 gn = (std::pow(a,(2.0/3.0))/fiss->akap)*densn/densg/pi*std::pow(nt,2); 02033 ga = (std::pow(a,(2.0/3.0))/fiss->akap)*densa/densg/pi*2.0*std::pow(at,2); 02034 gf = densf/densg/pi/2.0*ft; 02035 // assert(isnan(gf) == false); 02036 02037 // assert(isnan(gp) == false); 02038 // assert(isnan(gn) == false); 02039 // assert(isnan(ga) == false); 02040 // assert(isnan(ft) == false); 02041 // assert(ft != 0); 02042 // assert(isnan(gf) == false); 02043 02044 if(itest == 1) { 02045 G4cout <<"gn,gp,ga,gf " << gn <<","<< gp <<","<< ga <<","<< gf << G4endl; 02046 } 02047 } 02048 else { 02049 if(verboseLevel > 2) { 02050 G4cout <<"direct: densg <= 0.e0 " << a <<","<< zprf <<","<< ee << G4endl; 02051 } 02052 } 02053 02054 gsum = ga + gp + gn; 02055 // assert(isinf(gsum) == false); 02056 // assert(isnan(gsum) == false); 02057 if (gsum > 0.0) { 02058 ts1 = hbar / gsum; 02059 } 02060 else { 02061 ts1 = 1.0e99; 02062 } 02063 02064 if (inum > ilast) { // new event means reset the time scale 02065 tsum = 0; 02066 } 02067 02068 // calculate the relative probabilities for all decay channels 02069 if (densf == 0.0) { 02070 if (densp == 0.0) { 02071 if (densn == 0.0) { 02072 if (densa == 0.0) { 02073 // no reaction is possible 02074 probf = 0.0; 02075 probp = 0.0; 02076 probn = 0.0; 02077 proba = 0.0; 02078 goto direct70; 02079 } 02080 02081 // alpha evaporation is the only open channel 02082 rf = 0.0; 02083 rp = 0.0; 02084 rn = 0.0; 02085 ra = 1.0; 02086 goto direct50; 02087 } 02088 02089 // alpha emission and neutron emission 02090 rf = 0.0; 02091 rp = 0.0; 02092 rn = 1.0; 02093 ra = densa*2.0/densn*std::pow((at/nt),2); 02094 goto direct50; 02095 } 02096 // alpha, proton and neutron emission 02097 rf = 0.0; 02098 rp = 1.0; 02099 rn = densn/densp*std::pow((nt/pt),2); 02100 // assert(isnan(rn) == false); 02101 ra = densa*2.0/densp*std::pow((at/pt),2); 02102 // assert(isnan(ra) == false); 02103 goto direct50; 02104 } 02105 02106 // here fission has taken place 02107 rf = 1.0; 02108 02109 // cramers and weidenmueller factors for the dynamical hindrances of 02110 // fission 02111 if (fiss->bet <= 1.0e-16) { 02112 cf = 1.0; 02113 wf = 1.0; 02114 } 02115 else if (sbf > 0.0e0) { 02116 cf = cram(fiss->bet,fiss->homega); 02117 // if fission barrier ef=0.d0 then fission is the only possible 02118 // channel. to avoid std::log(0) in function tau 02119 // a.j. 7/28/93 02120 if (ef <= 0.0) { 02121 rp = 0.0; 02122 rn = 0.0; 02123 ra = 0.0; 02124 goto direct50; 02125 } 02126 else { 02127 // transient time tau() 02128 tauc = tau(fiss->bet,fiss->homega,ef,ft); 02129 // assert(isnan(tauc) == false); 02130 } 02131 wfex = (tauc - tsum)/ts1; 02132 02133 if (wfex < 0.0) { 02134 wf = 1.0; 02135 } 02136 else { 02137 wf = std::exp( -wfex); 02138 } 02139 } 02140 else { 02141 cf=1.0; 02142 wf=1.0; 02143 } 02144 02145 if(verboseLevel > 2) { 02146 G4cout <<"tsum,wf,cf " << tsum <<","<< wf <<","<< cf << G4endl; 02147 } 02148 02149 tsum = tsum + ts1; 02150 02151 // change by g.k. and a.h. 5.9.95 02152 tconst = 0.7; 02153 dconst = 12.0/std::sqrt(a); 02154 // assert(isnan(dconst) == false); 02155 nprf = a - zprf; 02156 02157 if (fiss->optshp >= 2) { //then 02158 parite(nprf,&parc); 02159 // assert(isnan(parc) == false); 02160 dconst = dconst*parc; 02161 } 02162 else { 02163 dconst= 0.0; 02164 } 02165 if ((ee <= 17.e0) && (fiss->optles == 1) && (iz >= 90) && (in >= 134)) { //then 02166 // constant changed to 5.0 accord to moretto & vandenbosch a.j. 19.3.96 02167 gngf = std::pow(a,(2.0/3.0))*tconst/10.0*std::exp((ef-sn+dconst)/tconst); 02168 02169 // if the excitation energy is so low that densn=0 ==> gn = 0 02170 // fission remains the only channel. 02171 // a. j. 10.1.94 02172 if (gn == 0.0) { 02173 rn = 0.0; 02174 rp = 0.0; 02175 ra = 0.0; 02176 } 02177 else { 02178 rn=gngf; 02179 // assert(isnan(rn) == false); 02180 rp=gngf*gp/gn; 02181 // assert(isnan(rp) == false); 02182 ra=gngf*ga/gn; 02183 // assert(isnan(ra) == false); 02184 } 02185 } else { 02186 // assert(isnan(cf) == false); 02187 // assert(isinf(gn) == false); 02188 // assert(isinf(gf) == false); 02189 // assert(isinf(cf) == false); 02190 assert(gn > 0 || (gf != 0 && cf != 0)); 02191 rn = gn/(gf*cf); 02192 // G4cout <<"rn = " << G4endl; 02193 // G4cout <<"gn = " << gn << " gf = " << gf << " cf = " << cf << G4endl; 02194 // assert(isnan(rn) == false); 02195 rp = gp/(gf*cf); 02196 // assert(isnan(rp) == false); 02197 ra = ga/(gf*cf); 02198 // assert(isnan(ra) == false); 02199 } 02200 direct50: 02201 // relative decay probabilities 02202 // assert(isnan(ra) == false); 02203 // assert(isnan(rp) == false); 02204 // assert(isnan(rn) == false); 02205 // assert(isnan(rf) == false); 02206 02207 denomi = rp+rn+ra+rf; 02208 // assert(isnan(denomi) == false); 02209 assert(denomi > 0); 02210 // decay probabilities after transient time 02211 probf = rf/denomi; 02212 // assert(isnan(probf) == false); 02213 probp = rp/denomi; 02214 // assert(isnan(probp) == false); 02215 probn = rn/denomi; 02216 // assert(isnan(probn) == false); 02217 proba = ra/denomi; 02218 // assert(isnan(proba) == false); 02219 // assert(isinf(proba) == false); 02220 02221 // new treatment of grange-weidenmueller factor, 5.1.2000, khs !!! 02222 02223 // decay probabilites with transient time included 02224 // assert(isnan(wf) == false); 02225 assert(std::fabs(probf) <= 1.0); 02226 probf = probf * wf; 02227 if(probf == 1.0) { 02228 probp = 0.0; 02229 probn = 0.0; 02230 proba = 0.0; 02231 } 02232 else { 02233 probp = probp * (wf + (1.e0-wf)/(1.e0-probf)); 02234 probn = probn * (wf + (1.e0-wf)/(1.e0-probf)); 02235 proba = proba * (wf + (1.e0-wf)/(1.e0-probf)); 02236 } 02237 direct70: 02238 // assert(isnan(probp) == false); 02239 // assert(isnan(probn) == false); 02240 // assert(isnan(probf) == false); 02241 // assert(isnan(proba) == false); 02242 ptotl = probp+probn+proba+probf; 02243 // assert(isnan(ptotl) == false); 02244 02245 ee = eer; 02246 ilast = inum; 02247 02248 // Return values: 02249 // assert(isnan(proba) == false); 02250 (*probp_par) = probp; 02251 (*probn_par) = probn; 02252 (*proba_par) = proba; 02253 (*probf_par) = probf; 02254 (*ptotl_par) = ptotl; 02255 (*sn_par) = sn; 02256 (*sbp_par) = sbp; 02257 (*sba_par) = sba; 02258 (*ecn_par) = ecn; 02259 (*ecp_par) = ecp; 02260 (*eca_par) = eca; 02261 (*bp_par) = bp; 02262 (*ba_par) = ba; 02263 }
| G4double G4Abla::ecoul | ( | G4double | z1, | |
| G4double | n1, | |||
| G4double | beta1, | |||
| G4double | z2, | |||
| G4double | n2, | |||
| G4double | beta2, | |||
| G4double | d | |||
| ) |
Definition at line 3484 of file G4Abla.cc.
Referenced by fissionDistri().
03485 { 03486 // Coulomb potential between two nuclei 03487 // surfaces are in a distance of d 03488 // in a tip to tip configuration 03489 03490 // approximate formulation 03491 // On input: Z1 nuclear charge of first nucleus 03492 // N1 number of neutrons in first nucleus 03493 // beta1 deformation of first nucleus 03494 // Z2 nuclear charge of second nucleus 03495 // N2 number of neutrons in second nucleus 03496 // beta2 deformation of second nucleus 03497 // d distance of surfaces of the nuclei 03498 03499 // G4double Z1,N1,beta1,Z2,N2,beta2,d,ecoul; 03500 G4double ecoul = 0; 03501 G4double dtot = 0; 03502 const G4double r0 = 1.16; 03503 03504 dtot = r0 * ( std::pow((z1+n1),0.33333) * (1.0+(2.0/3.0)*beta1) 03505 + std::pow((z2+n2),0.33333) * (1.0+(2.0/3.0)*beta2) ) + d; 03506 ecoul = z1 * z2 * 1.44 / dtot; 03507 03508 // assert(isnan(ecoul) == false); 03509 return ecoul; 03510 }
This function will calculate the liquid-drop nuclear mass for spheri configuration according to the preprint NUCLEAR GROUND-STATE MASSES and DEFORMATIONS by P. Mo"ller et al. from August 16, 1993 p. All constants are taken from this publication for consistency.
Definition at line 2542 of file G4Abla.cc.
References f(), mod(), CLHEP::detail::n, G4INCL::Math::pi, and utilabs().
Referenced by mglms().
02543 { 02544 // CHANGED TO CALCULATE TOTAL BINDING ENERGY INSTEAD OF MASS EXCESS. 02545 // SWITCH FOR PAIRING INCLUDED AS WELL. 02546 // BINDING = EFLMAC(IA,IZ,0,OPTSHP) 02547 // FORTRAN TRANSCRIPT OF /U/GREWE/LANG/EEX/FRLDM.C 02548 // A.J. 15.07.96 02549 02550 // this function will calculate the liquid-drop nuclear mass for spheri 02551 // configuration according to the preprint NUCLEAR GROUND-STATE 02552 // MASSES and DEFORMATIONS by P. M"oller et al. from August 16, 1993 p. 02553 // All constants are taken from this publication for consistency. 02554 02555 // Parameters: 02556 // a: nuclear mass number 02557 // z: nuclear charge 02558 // flag: 0 - return mass excess 02559 // otherwise - return pairing (= -1/2 dpn + 1/2 (Dp + Dn)) 02560 02561 G4double eflmacResult = 0.0; 02562 02563 G4int in = 0; 02564 G4double z = 0.0, n = 0.0, a = 0.0, av = 0.0, as = 0.0; 02565 G4double a0 = 0.0, c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0; 02566 G4double f = 0.0, ca = 0.0, w = 0.0, dp = 0.0, dn = 0.0, dpn = 0.0, efl = 0.0; 02567 G4double rmac = 0.0, bs = 0.0, h = 0.0, r0 = 0.0, kf = 0.0, ks = 0.0; 02568 G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0; 02569 G4double mh = 0.0, mn = 0.0, esq = 0.0, ael = 0.0, i = 0.0; 02570 G4double pi = 3.141592653589793238e0; 02571 02572 // fundamental constants 02573 // hydrogen-atom mass excess 02574 mh = 7.289034; 02575 02576 // neutron mass excess 02577 mn = 8.071431; 02578 02579 // electronic charge squared 02580 esq = 1.4399764; 02581 02582 // constants from considerations other than nucl. masses 02583 // electronic binding 02584 ael = 1.433e-5; 02585 02586 // proton rms radius 02587 rp = 0.8; 02588 02589 // nuclear radius constant 02590 r0 = 1.16; 02591 02592 // range of yukawa-plus-expon. potential 02593 ay = 0.68; 02594 02595 // range of yukawa function used to generate 02596 // nuclear charge distribution 02597 aden= 0.70; 02598 02599 // constants from considering odd-even mass differences 02600 // average pairing gap 02601 rmac= 4.80; 02602 02603 // neutron-proton interaction 02604 h = 6.6; 02605 02606 // wigner constant 02607 w = 30.0; 02608 02609 // adjusted parameters 02610 // volume energy 02611 av = 16.00126; 02612 02613 // volume asymmetry 02614 kv = 1.92240; 02615 02616 // surface energy 02617 as = 21.18466; 02618 02619 // surface asymmetry 02620 ks = 2.345; 02621 // a^0 constant 02622 a0 = 2.615; 02623 02624 // charge asymmetry 02625 ca = 0.10289; 02626 02627 // we will account for deformation by using the microscopic 02628 // corrections tabulated from p. 68ff */ 02629 bs = 1.0; 02630 02631 z = double(iz); 02632 a = double(ia); 02633 in = ia - iz; 02634 n = double(in); 02635 dn = rmac*bs/std::pow(n,(1.0/3.0)); 02636 dp = rmac*bs/std::pow(z,(1.0/3.0)); 02637 dpn = h/bs/std::pow(a,(2.0/3.0)); 02638 // assert(isnan(dpn) == false); 02639 02640 c1 = 3.0/5.0*esq/r0; 02641 // assert(isnan(c1) == false); 02642 // assert(isinf(c1) == false); 02643 02644 c4 = 5.0/4.0*std::pow((3.0/(2.0*pi)),(2.0/3.0)) * c1; 02645 // assert(isnan(c4) == false); 02646 // assert(isinf(c4) == false); 02647 02648 // assert(isnan(pi) == false); 02649 // assert(isnan(z) == false); 02650 // assert(isnan(a) == false); 02651 // assert(isnan(r0) == false); 02652 kf = std::pow((9.0*pi*z/(4.0*a)),(1.0/3.0))/r0; 02653 // assert(isnan(kf) == false); 02654 // assert(isinf(kf) == false); 02655 02656 f = -1.0/8.0*rp*rp*esq/std::pow(r0,3) * (145.0/48.0 - 327.0/2880.0*std::pow(kf,2) * std::pow(rp,2) + 1527.0/1209600.0*std::pow(kf,4) * std::pow(rp,4)); 02657 i = (n-z)/a; 02658 02659 x0 = r0 * std::pow(a,(1.0/3.0)) / ay; 02660 y0 = r0 * std::pow(a,(1.0/3.0)) / aden; 02661 02662 b1 = 1.0 - 3.0/(std::pow(x0,2)) + (1.0 + x0) * (2.0 + 3.0/x0 + 3.0/std::pow(x0,2)) * std::exp(-2.0*x0); 02663 02664 b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0*y0) + 21.0/(8.0 * std::pow(y0,3)) 02665 - 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2) 02666 + 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0)); 02667 02668 // now calulation of total binding energy a.j. 16.7.96 02669 02670 efl = -1.0 * av*(1.0 - kv*i*i)*a + as*(1.0 - ks*i*i)*b1 * std::pow(a,(2.0/3.0)) + a0 02671 + c1*z*z*b3/std::pow(a,(1.0/3.0)) - c4*std::pow(z,(4.0/3.0))/std::pow(a,(1.e0/3.e0)) 02672 + f*std::pow(z,2)/a -ca*(n-z) - ael * std::pow(z,(2.39e0)); 02673 02674 if ((in == iz) && (mod(in,2) == 1) && (mod(iz,2) == 1)) { 02675 // n and z odd and equal 02676 efl = efl + w*(utilabs(i)+1.e0/a); 02677 } 02678 else { 02679 efl= efl + w* utilabs(i); 02680 } 02681 02682 // pairing is made optional 02683 if (optshp >= 2) { 02684 // average pairing 02685 if ((mod(in,2) == 1) && (mod(iz,2) == 1)) { 02686 efl = efl - dpn; 02687 } 02688 if (mod(in,2) == 1) { 02689 efl = efl + dn; 02690 } 02691 if (mod(iz,2) == 1) { 02692 efl = efl + dp; 02693 } 02694 // end if for pairing term 02695 } 02696 02697 if (flag != 0) { 02698 eflmacResult = (0.5*(dn + dp) - 0.5*dpn); 02699 } 02700 else { 02701 eflmacResult = efl; 02702 } 02703 02704 return eflmacResult; 02705 }
| void G4Abla::evapora | ( | G4double | zprf, | |
| G4double | aprf, | |||
| G4double | ee, | |||
| G4double | jprf, | |||
| G4double * | zf_par, | |||
| G4double * | af_par, | |||
| G4double * | mtota_par, | |||
| G4double * | pleva_par, | |||
| G4double * | pxeva_par, | |||
| G4double * | pyeva_par, | |||
| G4int * | ff_par, | |||
| G4int * | inttype_par, | |||
| G4int * | inum_par | |||
| ) |
Main evaporation routine.
Definition at line 1181 of file G4Abla.cc.
References G4Volant::acv, barfit(), direct(), dmin1(), G4Ecld::ecgnz, G4Fb::efa, G4cout, G4endl, haz(), idnint(), G4Volant::iv, G4Fiss::optshp, G4Volant::pcv, standardRandom(), G4Volant::xcv, G4Volant::ycv, G4Volant::zcv, and G4Volant::zpcv.
Referenced by breakItUp().
01185 { 01186 G4double zf = (*zf_par); 01187 G4double af = (*af_par); 01188 G4double mtota = (*mtota_par); 01189 G4double pleva = (*pleva_par); 01190 G4double pxeva = (*pxeva_par); 01191 G4double pyeva = (*pyeva_par); 01192 G4int ff = (*ff_par); 01193 G4int inttype = (*inttype_par); 01194 G4int inum = (*inum_par); 01195 01196 // 533 C 01197 // 534 C INPUT: 01198 // 535 C 01199 // 536 C ZPRF, APRF, EE(EE IS MODIFIED!), JPRF 01200 // 537 C 01201 // 538 C PROJECTILE AND TARGET PARAMETERS + CROSS SECTIONS 01202 // 539 C COMMON /ABRAMAIN/ AP,ZP,AT,ZT,EAP,BETA,BMAXNUC,CRTOT,CRNUC, 01203 // 540 C R_0,R_P,R_T, IMAX,IRNDM,PI, 01204 // 541 C BFPRO,SNPRO,SPPRO,SHELL 01205 // 542 C 01206 // 543 C AP,ZP,AT,ZT - PROJECTILE AND TARGET MASSES 01207 // 544 C EAP,BETA - BEAM ENERGY PER NUCLEON, V/C 01208 // 545 C BMAXNUC - MAX. IMPACT PARAMETER FOR NUCL. REAC. 01209 // 546 C CRTOT,CRNUC - TOTAL AND NUCLEAR REACTION CROSS SECTION 01210 // 547 C R_0,R_P,R_T, - RADIUS PARAMETER, PROJECTILE+ TARGET RADII 01211 // 548 C IMAX,IRNDM,PI - MAXIMUM NUMBER OF EVENTS, DUMMY, 3.141... 01212 // 549 C BFPRO - FISSION BARRIER OF THE PROJECTILE 01213 // 550 C SNPRO - NEUTRON SEPARATION ENERGY OF THE PROJECTILE 01214 // 551 C SPPRO - PROTON " " " " " 01215 // 552 C SHELL - GROUND STATE SHELL CORRECTION 01216 // 553 C 01217 // 554 C--------------------------------------------------------------------- 01218 // 555 C FISSION BARRIERS 01219 // 556 C COMMON /FB/ EFA 01220 // 557 C EFA - ARRAY OF FISSION BARRIERS 01221 // 558 C--------------------------------------------------------------------- 01222 // 559 C OUTPUT: 01223 // 560 C ZF, AF, MTOTA, PLEVA, PTEVA, FF, INTTYPE, INUM 01224 // 561 C 01225 // 562 C ZF,AF - CHARGE AND MASS OF FINAL FRAGMENT AFTER EVAPORATION 01226 // 563 C MTOTA _ NUMBER OF EVAPORATED ALPHAS 01227 // 564 C PLEVA,PXEVA,PYEVA - MOMENTUM RECOIL BY EVAPORATION 01228 // 565 C INTTYPE - TYPE OF REACTION 0/1 NUCLEAR OR ELECTROMAGNETIC 01229 // 566 C FF - 0/1 NO FISSION / FISSION EVENT 01230 // 567 C INUM - EVENTNUMBER 01231 // 568 C ____________________________________________________________________ 01232 // 569 C / 01233 // 570 C / CALCUL DE LA MASSE ET CHARGE FINALES D'UNE CHAINE D'EVAPORATION 01234 // 571 C / 01235 // 572 C / PROCEDURE FOR CALCULATING THE FINAL MASS AND CHARGE VALUES OF A 01236 // 573 C / SPECIFIC EVAPORATION CHAIN, STARTING POINT DEFINED BY (APRF, ZPRF, 01237 // 574 C / EE) 01238 // 575 C / On ajoute les 3 composantes de l'impulsion (PXEVA,PYEVA,PLEVA) 01239 // 576 C / (actuellement PTEVA n'est pas correct; mauvaise norme...) 01240 // 577 C /____________________________________________________________________ 01241 // 578 C 01242 // 612 C 01243 // 613 C----------------------------------------------------------------------- 01244 // 614 C IRNDM DUMMY ARGUMENT FOR RANDOM-NUMBER FUNCTION 01245 // 615 C SORTIE LOCAL HELP VARIABLE TO END THE EVAPORATION CHAIN 01246 // 616 C ZF NUCLEAR CHARGE OF THE FRAGMENT 01247 // 617 C ZPRF NUCLEAR CHARGE OF THE PREFRAGMENT 01248 // 618 C AF MASS NUMBER OF THE FRAGMENT 01249 // 619 C APRF MASS NUMBER OF THE PREFRAGMENT 01250 // 620 C EPSILN ENERGY BURNED IN EACH EVAPORATION STEP 01251 // 621 C MALPHA LOCAL MASS CONTRIBUTION TO MTOTA IN EACH EVAPORATION 01252 // 622 C STEP 01253 // 623 C EE EXCITATION ENERGY (VARIABLE) 01254 // 624 C PROBP PROTON EMISSION PROBABILITY 01255 // 625 C PROBN NEUTRON EMISSION PROBABILITY 01256 // 626 C PROBA ALPHA-PARTICLE EMISSION PROBABILITY 01257 // 627 C PTOTL TOTAL EMISSION PROBABILITY 01258 // 628 C E LOWEST PARTICLE-THRESHOLD ENERGY 01259 // 629 C SN NEUTRON SEPARATION ENERGY 01260 // 630 C SBP PROTON SEPARATION ENERGY PLUS EFFECTIVE COULOMB 01261 // 631 C BARRIER 01262 // 632 C SBA ALPHA-PARTICLE SEPARATION ENERGY PLUS EFFECTIVE 01263 // 633 C COULOMB BARRIER 01264 // 634 C BP EFFECTIVE PROTON COULOMB BARRIER 01265 // 635 C BA EFFECTIVE ALPHA COULOMB BARRIER 01266 // 636 C MTOTA TOTAL MASS OF THE EVAPORATED ALPHA PARTICLES 01267 // 637 C X UNIFORM RANDOM NUMBER FOR NUCLEAR CHARGE 01268 // 638 C AMOINS LOCAL MASS NUMBER OF EVAPORATED PARTICLE 01269 // 639 C ZMOINS LOCAL NUCLEAR CHARGE OF EVAPORATED PARTICLE 01270 // 640 C ECP KINETIC ENERGY OF PROTON WITHOUT COULOMB 01271 // 641 C REPULSION 01272 // 642 C ECN KINETIC ENERGY OF NEUTRON 01273 // 643 C ECA KINETIC ENERGY OF ALPHA PARTICLE WITHOUT COULOMB 01274 // 644 C REPULSION 01275 // 645 C PLEVA TRANSVERSAL RECOIL MOMENTUM OF EVAPORATION 01276 // 646 C PTEVA LONGITUDINAL RECOIL MOMENTUM OF EVAPORATION 01277 // 647 C FF FISSION FLAG 01278 // 648 C INTTYPE INTERACTION TYPE FLAG 01279 // 649 C RNDX RECOIL MOMENTUM IN X-DIRECTION IN A SINGLE STEP 01280 // 650 C RNDY RECOIL MOMENTUM IN Y-DIRECTION IN A SINGLE STEP 01281 // 651 C RNDZ RECOIL MOMENTUM IN Z-DIRECTION IN A SINGLE STEP 01282 // 652 C RNDN NORMALIZATION OF RECOIL MOMENTUM FOR EACH STEP 01283 // 653 C----------------------------------------------------------------------- 01284 // 654 C 01285 // 655 SAVE 01286 // SAVE -> static 01287 01288 static G4int sortie = 0; 01289 static G4double epsiln = 0.0, probp = 0.0, probn = 0.0, proba = 0.0, ptotl = 0.0, e = 0.0; 01290 static G4double sn = 0.0, sbp = 0.0, sba = 0.0, x = 0.0, amoins = 0.0, zmoins = 0.0; 01291 G4double ecn = 0.0, ecp = 0.0,eca = 0.0, bp = 0.0, ba = 0.0; 01292 static G4double pteva = 0.0; 01293 01294 static G4int itest = 0; 01295 static G4double probf = 0.0; 01296 01297 static G4int k = 0, j = 0, il = 0; 01298 01299 static G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0; 01300 static G4double sbfis = 0.0, rnd = 0.0; 01301 static G4double selmax = 0.0; 01302 static G4double segs = 0.0; 01303 static G4double ef = 0.0; 01304 static G4int irndm = 0; 01305 01306 static G4double pc = 0.0, malpha = 0.0; 01307 01308 zf = zprf; 01309 af = aprf; 01310 pleva = 0.0; 01311 pteva = 0.0; 01312 pxeva = 0.0; 01313 pyeva = 0.0; 01314 01315 sortie = 0; 01316 ff = 0; 01317 01318 itest = 0; 01319 if (itest == 1) { 01320 G4cout << "***************************" << G4endl; 01321 } 01322 01323 evapora10: 01324 01325 if (itest == 1) { 01326 G4cout <<"------zf,af,ee------" << idnint(zf) << "," << idnint(af) << "," << ee << G4endl; 01327 } 01328 01329 // calculation of the probabilities for the different decay channels 01330 // plus separation energies and kinetic energies of the particles 01331 direct(zf,af,ee,jprf,&probp,&probn,&proba,&probf,&ptotl, 01332 &sn,&sbp,&sba,&ecn,&ecp,&eca,&bp,&ba,inttype,inum,itest); //:::FIXME::: Call 01333 // assert(isnan(proba) == false); 01334 // assert(isnan(probp) == false); 01335 // assert(isnan(probn) == false); 01336 // assert(isnan(probf) == false); 01337 assert((eca+ba) >= 0); 01338 assert((ecp+bp) >= 0); 01339 // assert(isnan(ecp) == false); 01340 // assert(isnan(ecn) == false); 01341 // assert(isnan(bp) == false); 01342 // assert(isnan(ba) == false); 01343 k = idnint(zf); 01344 j = idnint(af-zf); 01345 01346 // now ef is calculated from efa that depends on the subroutine 01347 // barfit which takes into account the modification on the ang. mom. 01348 // jb mvr 6-aug-1999 01349 // note *** shell correction! (ecgnz) jb mvr 20-7-1999 01350 il = idnint(jprf); 01351 barfit(k,k+j,il,&sbfis,&segs,&selmax); 01352 // assert(isnan(sbfis) == false); 01353 01354 if ((fiss->optshp == 1) || (fiss->optshp == 3)) { //then 01355 // fb->efa[k][j] = G4double(sbfis) - ecld->ecgnz[j][k]; 01356 fb->efa[j][k] = G4double(sbfis) - ecld->ecgnz[j][k]; 01357 } 01358 else { 01359 fb->efa[j][k] = G4double(sbfis); 01360 // fb->efa[j][k] = G4double(sbfis); 01361 } //end if 01362 ef = fb->efa[j][k]; 01363 // ef = fb->efa[j][k]; 01364 // assert(isnan(fb->efa[j][k]) == false); 01365 // here the final steps of the evaporation are calculated 01366 if ((sortie == 1) || (ptotl == 0.e0)) { 01367 e = dmin1(sn,sbp,sba); 01368 if (e > 1.0e30) { 01369 if(verboseLevel > 2) { 01370 G4cout <<"erreur a la sortie evapora,e>1.e30,af=" << af <<" zf=" << zf << G4endl; 01371 } 01372 } 01373 if (zf <= 6.0) { 01374 goto evapora100; 01375 } 01376 if (e < 0.0) { 01377 if (sn == e) { 01378 af = af - 1.e0; 01379 } 01380 else if (sbp == e) { 01381 af = af - 1.0; 01382 zf = zf - 1.0; 01383 } 01384 else if (sba == e) { 01385 af = af - 4.0; 01386 zf = zf - 2.0; 01387 } 01388 if (af < 2.5) { 01389 goto evapora100; 01390 } 01391 goto evapora10; 01392 } 01393 goto evapora100; 01394 } 01395 irndm = irndm + 1; 01396 01397 // here the normal evaporation cascade starts 01398 01399 // random number for the evaporation 01400 // x = double(Rndm(irndm))*ptotl; 01401 x = double(haz(1))*ptotl; 01402 01403 // G4cout <<"proba = " << proba << G4endl; 01404 // G4cout <<"probp = " << probp << G4endl; 01405 // G4cout <<"probn = " << probn << G4endl; 01406 // G4cout <<"probf = " << probf << G4endl; 01407 01408 itest = 0; 01409 if (x < proba) { 01410 // alpha evaporation 01411 if (itest == 1) { 01412 G4cout <<"< alpha evaporation >" << G4endl; 01413 } 01414 amoins = 4.0; 01415 zmoins = 2.0; 01416 epsiln = sba + eca; 01417 assert((std::pow((1.0 + (eca+ba)/3.72834e3),2) - 1.0) >= 0); 01418 pc = std::sqrt(std::pow((1.0 + (eca+ba)/3.72834e3),2) - 1.0) * 3.72834e3; 01419 // assert(isnan(pc) == false); 01420 malpha = 4.0; 01421 01422 // volant: 01423 volant->iv = volant->iv + 1; 01424 volant->acv[volant->iv] = 4.; 01425 volant->zpcv[volant->iv] = 2.; 01426 volant->pcv[volant->iv] = pc; 01427 } 01428 else if (x < proba+probp) { 01429 // proton evaporation 01430 if (itest == 1) { 01431 G4cout <<"< proton evaporation >" << G4endl; 01432 } 01433 amoins = 1.0; 01434 zmoins = 1.0; 01435 epsiln = sbp + ecp; 01436 assert((std::pow((1.0 + (ecp + bp)/9.3827e2),2) - 1.0) >= 0); 01437 pc = std::sqrt(std::pow((1.0 + (ecp + bp)/9.3827e2),2) - 1.0) * 9.3827e2; 01438 // assert(isnan(pc) == false); 01439 malpha = 0.0; 01440 // volant: 01441 volant->iv = volant->iv + 1; 01442 volant->acv[volant->iv] = 1.0; 01443 volant->zpcv[volant->iv] = 1.; 01444 volant->pcv[volant->iv] = pc; 01445 } 01446 else if (x < proba+probp+probn) { 01447 // neutron evaporation 01448 if (itest == 1) { 01449 G4cout <<"< neutron evaporation >" << G4endl; 01450 } 01451 amoins = 1.0; 01452 zmoins = 0.0; 01453 epsiln = sn + ecn; 01454 assert((std::pow((1.0 + (ecn)/9.3956e2),2) - 1.0) >= 0); 01455 pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956e2),2) - 1.0) * 9.3956e2; 01456 // assert(isnan(pc) == false); 01457 malpha = 0.0; 01458 01459 // volant: 01460 volant->iv = volant->iv + 1; 01461 volant->acv[volant->iv] = 1.; 01462 volant->zpcv[volant->iv] = 0.; 01463 volant->pcv[volant->iv] = pc; 01464 } 01465 else { 01466 // fission 01467 // in case of fission-events the fragment nucleus is the mother nucleus 01468 // before fission occurs with excitation energy above the fis.- barrier. 01469 // fission fragment mass distribution is calulated in subroutine fisdis 01470 if (itest == 1) { 01471 G4cout <<"< fission >" << G4endl; 01472 } 01473 amoins = 0.0; 01474 zmoins = 0.0; 01475 epsiln = ef; 01476 01477 malpha = 0.0; 01478 pc = 0.0; 01479 ff = 1; 01480 // ff = 0; // For testing, allows to disable fission! 01481 } 01482 01483 if (itest == 1) { 01484 G4cout <<"sn,sbp,sba,ef" << sn << "," << sbp << "," << sba <<"," << ef << G4endl; 01485 G4cout <<"probn,probp,proba,probf,ptotl " <<","<< probn <<","<< probp <<","<< proba <<","<< probf <<","<< ptotl << G4endl; 01486 } 01487 01488 // calculation of the daughter nucleus 01489 af = af - amoins; 01490 zf = zf - zmoins; 01491 ee = ee - epsiln; 01492 if (ee <= 0.01) { 01493 ee = 0.01; 01494 } 01495 mtota = mtota + malpha; 01496 01497 if(ff == 0) { 01498 standardRandom(&rnd,&(hazard->igraine[8])); 01499 ctet1 = 2.0*rnd - 1.0; 01500 standardRandom(&rnd,&(hazard->igraine[4])); 01501 phi1 = rnd*2.0*3.141592654; 01502 stet1 = std::sqrt(1.0 - std::pow(ctet1,2)); 01503 // assert(isnan(stet1) == false); 01504 volant->xcv[volant->iv] = stet1*std::cos(phi1); 01505 volant->ycv[volant->iv] = stet1*std::sin(phi1); 01506 volant->zcv[volant->iv] = ctet1; 01507 pxeva = pxeva - pc * volant->xcv[volant->iv]; 01508 pyeva = pyeva - pc * volant->ycv[volant->iv]; 01509 pleva = pleva - pc * ctet1; 01510 // assert(isnan(pleva) == false); 01511 } 01512 01513 // condition for end of evaporation 01514 if ((af < 2.5) || (ff == 1)) { 01515 goto evapora100; 01516 } 01517 goto evapora10; 01518 01519 evapora100: 01520 (*zf_par) = zf; 01521 (*af_par) = af; 01522 (*mtota_par) = mtota; 01523 (*pleva_par) = pleva; 01524 (*pxeva_par) = pxeva; 01525 (*pyeva_par) = pyeva; 01526 (*ff_par) = ff; 01527 (*inttype_par) = inttype; 01528 (*inum_par) = inum; 01529 01530 return; 01531 }
Definition at line 3392 of file G4Abla.cc.
Referenced by fissionDistri().
03393 { 03394 // Procedure to calculate I_OUT from R_IN in a way that 03395 // on the average a flat distribution in R_IN results in a 03396 // fluctuating distribution in I_OUT with an even-odd effect as 03397 // given by R_EVEN_ODD 03398 03399 // /* ------------------------------------------------------------ */ 03400 // /* EXAMPLES : */ 03401 // /* ------------------------------------------------------------ */ 03402 // /* If R_EVEN_ODD = 0 : */ 03403 // /* CEIL(R_IN) ---- */ 03404 // /* */ 03405 // /* R_IN -> */ 03406 // /* (somewhere in between CEIL(R_IN) and FLOOR(R_IN)) */ */ 03407 // /* */ 03408 // /* FLOOR(R_IN) ---- --> I_OUT */ 03409 // /* ------------------------------------------------------------ */ 03410 // /* If R_EVEN_ODD > 0 : */ 03411 // /* The interval for the above treatment is */ 03412 // /* larger for FLOOR(R_IN) = even and */ 03413 // /* smaller for FLOOR(R_IN) = odd */ 03414 // /* For R_EVEN_ODD < 0 : just opposite treatment */ 03415 // /* ------------------------------------------------------------ */ 03416 03417 // /* ------------------------------------------------------------ */ 03418 // /* On input: R_ORIGIN nuclear charge (real number) */ 03419 // /* R_EVEN_ODD requested even-odd effect */ 03420 // /* Intermediate quantity: R_IN = R_ORIGIN + 0.5 */ 03421 // /* On output: I_OUT nuclear charge (integer) */ 03422 // /* ------------------------------------------------------------ */ 03423 03424 // G4double R_ORIGIN,R_IN,R_EVEN_ODD,R_REST,R_HELP; 03425 G4double r_in = 0.0, r_rest = 0.0, r_help = 0.0; 03426 G4double r_floor = 0.0; 03427 G4double r_middle = 0.0; 03428 // G4int I_OUT,N_FLOOR; 03429 G4int n_floor = 0; 03430 03431 r_in = r_origin + 0.5; 03432 r_floor = (float)((int)(r_in)); 03433 if (r_even_odd < 0.001) { 03434 i_out = (int)(r_floor); 03435 } 03436 else { 03437 r_rest = r_in - r_floor; 03438 r_middle = r_floor + 0.5; 03439 n_floor = (int)(r_floor); 03440 if (n_floor%2 == 0) { 03441 // even before modif. 03442 r_help = r_middle + (r_rest - 0.5) * (1.0 - r_even_odd); 03443 } 03444 else { 03445 // odd before modification 03446 r_help = r_middle + (r_rest - 0.5) * (1.0 + r_even_odd); 03447 } 03448 i_out = (int)(r_help); 03449 } 03450 }
Definition at line 1147 of file G4Abla.cc.
Referenced by direct(), and spdef().
01148 { 01149 // CALCULATION OF FISSILITY PARAMETER 01150 // 01151 // INPUT: A,Z INTEGER MASS & CHARGE OF NUCLEUS 01152 // OPTXFIS = 0 : MYERS, SWIATECKI 01153 // 1 : DAHLINGER 01154 // 2 : ANDREYEV 01155 01156 G4double aa = 0.0, zz = 0.0, i = 0.0; 01157 G4double fissilityResult = 0.0; 01158 01159 aa = double(a); 01160 zz = double(z); 01161 i = double(a-2*z) / aa; 01162 01163 // myers & swiatecki droplet modell 01164 if (optxfis == 0) { //then 01165 fissilityResult = std::pow(zz,2) / aa /50.8830e0 / (1.0e0 - 1.7826e0 * std::pow(i,2)); 01166 } 01167 01168 if (optxfis == 1) { 01169 // dahlinger fit: 01170 fissilityResult = std::pow(zz,2) / aa * std::pow((49.22e0*(1.e0 - 0.3803e0*std::pow(i,2) - 20.489e0*std::pow(i,4))),(-1)); 01171 } 01172 01173 if (optxfis == 2) { 01174 // dubna fit: 01175 fissilityResult = std::pow(zz,2) / aa /(48.e0*(1.e0 - 17.22e0*std::pow(i,4))); 01176 } 01177 01178 return fissilityResult; 01179 }
| void G4Abla::fissionDistri | ( | G4double & | a, | |
| G4double & | z, | |||
| G4double & | e, | |||
| G4double & | a1, | |||
| G4double & | z1, | |||
| G4double & | e1, | |||
| G4double & | v1, | |||
| G4double & | a2, | |||
| G4double & | z2, | |||
| G4double & | e2, | |||
| G4double & | v2 | |||
| ) |
Definition at line 3512 of file G4Abla.cc.
References ecoul(), even_odd(), G4cout, G4endl, gausshaz(), haz(), max(), min(), CLHEP::detail::n, and umass().
03515 { 03516 // On input: A, Z, E (mass, atomic number and exc. energy of compound nucleus 03517 // before fission) 03518 // On output: Ai, Zi, Ei (mass, atomic number and exc. energy of fragment 1 and 2 03519 // after fission) 03520 03521 // Additionally calculated but not put in the parameter list: 03522 // Kinetic energy of prefragments EkinR1, EkinR2 03523 03524 // Translation of SIMFIS18.PLI (KHS, 2.1.2001) 03525 03526 // This program calculates isotopic distributions of fission fragments 03527 // with a semiempirical model 03528 // Copy from SIMFIS3, KHS, 8. February 1995 03529 // Modifications made by Jose Benlliure and KHS in August 1996 03530 // Energy counted from lowest barrier (J. Benlliure, KHS 1997) 03531 // Some bugs corrected (J. Benlliure, KHS 1997) 03532 // Version used for thesis S. Steinhaueser (August 1997) 03533 // (Curvature of LD potential increased by factor of 2!) 03534 03535 // Weiter veraendert mit der Absicht, eine Version zu erhalten, die 03536 // derjenigen entspricht, die von J. Benlliure et al. 03537 // in Nucl. Phys. A 628 (1998) 458 verwendet wurde, 03538 // allerdings ohne volle Neutronenabdampfung. 03539 03540 // The excitation energy was calculate now for each fission channel 03541 // separately. The dissipation from saddle to scission was taken from 03542 // systematics, the deformation energy at scission considers the shell 03543 // effects in a simplified way, and the fluctuation is included. 03544 // KHS, April 1999 03545 03546 // The width in N/Z was carefully adapted to values given by Lang et al. 03547 03548 // The width and eventually a shift in N/Z (polarization) follows the 03549 // following rules: 03550 03551 // The line N/Z following UCD has an angle of std::atan(Zcn/Ncn) 03552 // to the horizontal axis on a chart of nuclides. 03553 // (For 238U the angle is 32.2 deg.) 03554 03555 // The following relations hold: (from Armbruster) 03556 // 03557 // sigma(N) (A=const) = sigma(Z) (A=const) 03558 // sigma(A) (N=const) = sigma(Z) (N=const) 03559 // sigma(A) (Z=const) = sigma(N) (Z=const) 03560 // 03561 // From this we get: 03562 // sigma(Z) (N=const) * N = sigma(N) (Z=const) * Z 03563 // sigma(A) (Z=const) = sigma(Z) (A=const) * A/Z 03564 // sigma(N) (Z=const) = sigma(Z) (A=const) * A/Z 03565 // Z*sigma(N) (Z=const) = N*sigma(Z) (N=const) = A*sigma(Z) (A=const) 03566 03567 // Excitation energy now calculated above the lowest potential point 03568 // Inclusion of a distribution of excitation energies 03569 03570 // Several modifications, starting from SIMFIS12: KHS November 2000 03571 // This version seems to work quite well for 238U. 03572 // The transition from symmetric to asymmetric fission around 226Th 03573 // is reasonably well reproduced, although St. I is too strong and St. II 03574 // is too weak. St. I and St. II are also weakly seen for 208Pb. 03575 03576 // Extensions for an event generator of fission events (21.11.2000,KHS) 03577 03578 // Defalt parameters (IPARS) rather carefully adjusted to 03579 // pre-neutron mass distributions of Vives et al. (238U + n) 03580 // Die Parameter Fgamma1 und Fgamma2 sind kleiner als die resultierenden 03581 // Breiten der Massenverteilungen!!! 03582 // Fgamma1 und Fgamma2 wurden angepa�, so da� 03583 // Sigma-A(ST-I) = 3.3, Sigma-A(St-II) = 5.8 (nach Vives) 03584 03585 // Parameters of the model carefully adjusted by KHS (2.2.2001) to 03586 // 238U + 208Pb, 1000 A MeV, Timo Enqvist et al. 03587 03588 03589 G4double n = 0.0; 03590 G4double nlight1 = 0.0, nlight2 = 0.0; 03591 G4double aheavy1 = 0.0,alight1 = 0.0, aheavy2 = 0.0, alight2 = 0.0; 03592 G4double eheavy1 = 0.0, elight1 = 0.0, eheavy2 = 0.0, elight2 = 0.0; 03593 G4double zheavy1_shell = 0.0, zheavy2_shell = 0.0; 03594 G4double zlight1 = 0.0, zlight2 = 0.0; 03595 G4double masscurv = 0.0; 03596 G4double sasymm1 = 0.0, sasymm2 = 0.0, ssymm = 0.0, ysum = 0.0, yasymm = 0.0; 03597 G4double ssymm_mode1 = 0.0, ssymm_mode2 = 0.0; 03598 G4double cz_asymm1_saddle = 0.0, cz_asymm2_saddle = 0.0; 03599 // Curvature at saddle, modified by ld-potential 03600 G4double wzasymm1_saddle, wzasymm2_saddle, wzsymm_saddle = 0.0; 03601 G4double wzasymm1_scission = 0.0, wzasymm2_scission = 0.0, wzsymm_scission = 0.0; 03602 G4double wzasymm1 = 0.0, wzasymm2 = 0.0, wzsymm = 0.0; 03603 G4double nlight1_eff = 0.0, nlight2_eff = 0.0; 03604 G4int imode = 0; 03605 G4double rmode = 0.0; 03606 G4double z1mean = 0.0, z2mean = 0.0, z1width = 0.0, za1width = 0.0; 03607 // G4double Z1,Z2,N1R,N2R,A1R,A2R,N1,N2,A1,A2; 03608 G4double n1r = 0.0, n2r = 0.0, a1r = 0.0, a2r = 0.0, n1 = 0.0, n2 = 0.0; 03609 03610 G4double zsymm = 0.0, nsymm = 0.0, asymm = 0.0; 03611 G4double n1mean = 0.0, n2mean, n1width; 03612 G4double dueff = 0.0; 03613 // effective shell effect at lowest barrier 03614 G4double eld = 0.0; 03615 // Excitation energy with respect to ld barrier 03616 G4double re1 = 0.0, re2 = 0.0, re3 = 0.0; 03617 G4double eps1 = 0.0, eps2 = 0.0; 03618 G4double n1ucd = 0.0, n2ucd = 0.0, z1ucd = 0.0, z2ucd = 0.0; 03619 G4double beta = 0.0, beta1 = 0.0, beta2 = 0.0; 03620 03621 G4double dn1_pol = 0.0; 03622 // shift of most probable neutron number for given Z, 03623 // according to polarization 03624 G4int i_help = 0; 03625 03626 // /* Parameters of the semiempirical fission model */ 03627 G4double a_levdens = 0.0; 03628 // /* level-density parameter */ 03629 G4double a_levdens_light1 = 0.0, a_levdens_light2 = 0.0; 03630 G4double a_levdens_heavy1 = 0.0, a_levdens_heavy2 = 0.0; 03631 const G4double r_null = 1.16; 03632 // /* radius parameter */ 03633 G4double epsilon_1_saddle = 0.0, epsilon0_1_saddle = 0.0; 03634 G4double epsilon_2_saddle = 0.0, epsilon0_2_saddle = 0.0, epsilon_symm_saddle = 0.0; 03635 G4double epsilon_1_scission = 0.0, epsilon0_1_scission = 0.0; 03636 G4double epsilon_2_scission = 0.0, epsilon0_2_scission = 0.0; 03637 G4double epsilon_symm_scission = 0.0; 03638 // /* modified energy */ 03639 G4double e_eff1_saddle = 0.0, e_eff2_saddle = 0.0; 03640 G4double epot0_mode1_saddle = 0.0, epot0_mode2_saddle = 0.0, epot0_symm_saddle = 0.0; 03641 G4double epot_mode1_saddle = 0.0, epot_mode2_saddle = 0.0, epot_symm_saddle = 0.0; 03642 G4double e_defo = 0.0, e_defo1 = 0.0, e_defo2 = 0.0, e_scission = 0.0, e_asym = 0.0; 03643 G4double e1exc = 0.0, e2exc = 0.0; 03644 G4double e1exc_sigma = 0.0, e2exc_sigma = 0.0; 03645 G4double e1final = 0.0, e2final = 0.0; 03646 03647 const G4double r0 = 1.16; 03648 G4double tker = 0.0; 03649 G4double ekin1 = 0.0, ekin2 = 0.0; 03650 // G4double EkinR1,EkinR2,E1,E2,V1,V2; 03651 G4double ekinr1 = 0.0, ekinr2 = 0.0; 03652 G4int icz = 0, k = 0; 03653 03654 // Input parameters: 03655 //OMMENT(Nuclear charge number); 03656 // G4double Z; 03657 //OMMENT(Nuclear mass number); 03658 // G4double A; 03659 //OMMENT(Excitation energy above fission barrier); 03660 // G4double E; 03661 03662 // Model parameters: 03663 //OMMENT(position of heavy peak valley 1); 03664 const G4double nheavy1 = 83.0; 03665 //OMMENT(position of heavy peak valley 2); 03666 const G4double nheavy2 = 90.0; 03667 //OMMENT(Shell effect for valley 1); 03668 const G4double delta_u1_shell = -2.65; 03669 // Parameter (Delta_U1_shell = -2) 03670 //OMMENT(Shell effect for valley 2); 03671 const G4double delta_u2_shell = -3.8; 03672 // Parameter (Delta_U2_shell = -3.2) 03673 //OMMENT(I: used shell effect); 03674 G4double delta_u1 = 0.0; 03675 //omment(I: used shell effect); 03676 G4double delta_u2 = 0.0; 03677 //OMMENT(Curvature of asymmetric valley 1); 03678 const G4double cz_asymm1_shell = 0.7; 03679 //OMMENT(Curvature of asymmetric valley 2); 03680 const G4double cz_asymm2_shell = 0.15; 03681 //OMMENT(Factor for width of distr. valley 1); 03682 const G4double fwidth_asymm1 = 0.63; 03683 //OMMENT(Factor for width of distr. valley 2); 03684 const G4double fwidth_asymm2 = 0.97; 03685 // Parameter (CZ_asymm2_scission = 0.12) 03686 //OMMENT(Parameter x: a = A/x); 03687 const G4double xlevdens = 12.0; 03688 //OMMENT(Factor to gamma_heavy1); 03689 const G4double fgamma1 = 2.0; 03690 //OMMENT(I: fading of shells (general)); 03691 G4double gamma = 0.0; 03692 //OMMENT(I: fading of shell 1); 03693 G4double gamma_heavy1 = 0.0; 03694 //OMMENT(I: fading of shell 2); 03695 G4double gamma_heavy2 = 0.0; 03696 //OMMENT(Zero-point energy at saddle); 03697 const G4double e_zero_point = 0.5; 03698 //OMMENT(I: friction from saddle to scission); 03699 G4double e_saddle_scission = 0.0; 03700 //OMMENT(Friction factor); 03701 const G4double friction_factor = 1.0; 03702 //OMMENT(I: Internal counter for different modes); INIT(0,0,0) 03703 // Integer*4 I_MODE(3) 03704 //OMMENT(I: Yield of symmetric mode); 03705 G4double ysymm = 0.0; 03706 //OMMENT(I: Yield of asymmetric mode 1); 03707 G4double yasymm1 = 0.0; 03708 //OMMENT(I: Yield of asymmetric mode 2); 03709 G4double yasymm2 = 0.0; 03710 //OMMENT(I: Effective position of valley 1); 03711 G4double nheavy1_eff = 0.0; 03712 //OMMENT(I: position of heavy peak valley 1); 03713 G4double zheavy1 = 0.0; 03714 //omment(I: Effective position of valley 2); 03715 G4double nheavy2_eff = 0.0; 03716 //OMMENT(I: position of heavy peak valley 2); 03717 G4double zheavy2 = 0.0; 03718 //omment(I: Excitation energy above saddle 1); 03719 G4double eexc1_saddle = 0.0; 03720 //omment(I: Excitation energy above saddle 2); 03721 G4double eexc2_saddle = 0.0; 03722 //omment(I: Excitation energy above lowest saddle); 03723 G4double eexc_max = 0.0; 03724 //omment(I: Effective mass mode 1); 03725 G4double aheavy1_mean = 0.0; 03726 //omment(I: Effective mass mode 2); 03727 G4double aheavy2_mean = 0.0; 03728 //omment(I: Width of symmetric mode); 03729 G4double wasymm_saddle = 0.0; 03730 //OMMENT(I: Width of asymmetric mode 1); 03731 G4double waheavy1_saddle = 0.0; 03732 //OMMENT(I: Width of asymmetric mode 2); 03733 G4double waheavy2_saddle = 0.0; 03734 //omment(I: Width of symmetric mode); 03735 G4double wasymm = 0.0; 03736 //OMMENT(I: Width of asymmetric mode 1); 03737 G4double waheavy1 = 0.0; 03738 //OMMENT(I: Width of asymmetric mode 2); 03739 G4double waheavy2 = 0.0; 03740 //OMMENT(I: Even-odd effect in Z); 03741 G4double r_e_o = 0.0, r_e_o_exp = 0.0; 03742 //OMMENT(I: Curveture of symmetric valley); 03743 G4double cz_symm = 0.0; 03744 //OMMENT(I: Curvature of mass distribution for fixed Z); 03745 G4double cn = 0.0; 03746 //OMMENT(I: Curvature of Z distribution for fixed A); 03747 G4double cz = 0.0; 03748 //OMMENT(Minimum neutron width for constant Z); 03749 const G4double sigzmin = 1.16; 03750 //OMMENT(Surface distance of scission configuration); 03751 const G4double d = 2.0; 03752 03753 // /* Charge polarisation from Wagemanns p. 397: */ 03754 //OMMENT(Charge polarisation standard I); 03755 const G4double cpol1 = 0.65; 03756 //OMMENT(Charge polarisation standard II); 03757 const G4double cpol2 = 0.55; 03758 //OMMENT(=1: Polarisation simult. in N and Z); 03759 const G4int nzpol = 1; 03760 //OMMENT(=1: test output, =0: no test output); 03761 const G4int itest = 0; 03762 03763 // G4double UMASS, ECOUL, reps1, reps2, rn1_pol; 03764 G4double reps1 = 0.0, reps2 = 0.0, rn1_pol = 0.0; 03765 // Float_t HAZ,GAUSSHAZ; 03766 G4int kkk = 0; 03767 // G4int kkk = 10; // PK 03768 03769 // I_MODE = 0; 03770 03771 if(itest == 1) { 03772 G4cout << " cn mass " << a << G4endl; 03773 G4cout << " cn charge " << z << G4endl; 03774 G4cout << " cn energy " << e << G4endl; 03775 } 03776 03777 // /* average Z of asymmetric and symmetric components: */ 03778 n = a - z; /* neutron number of the fissioning nucleus */ 03779 03780 k = 0; 03781 icz = 0; 03782 if ( (std::pow(z,2)/a < 25.0) || (n < nheavy2) || (e > 500.0) ) { 03783 icz = -1; 03784 // GOTO 1002; 03785 goto milledeux; 03786 } 03787 03788 nlight1 = n - nheavy1; 03789 nlight2 = n - nheavy2; 03790 03791 // /* Polarisation assumed for standard I and standard II: 03792 // Z - Zucd = cpol (for A = const); 03793 // from this we get (see Armbruster) 03794 // Z - Zucd = Acn/Ncn * cpol (for N = const) */ 03795 03796 zheavy1_shell = ((nheavy1/n) * z) - ((a/n) * cpol1); 03797 zheavy2_shell = ((nheavy2/n) * z) - ((a/n) * cpol2); 03798 03799 e_saddle_scission = 03800 (-24.0 + 0.02227 * (std::pow(z,2))/(std::pow(a,0.33333)) ) * friction_factor; 03801 03802 // /* Energy dissipated from saddle to scission */ 03803 // /* F. Rejmund et al., Nucl. Phys. A 678 (2000) 215, fig. 4 b */ 03804 // E_saddle_scission = DMAX1(0.,E_saddle_scission); 03805 if (e_saddle_scission > 0.) { 03806 e_saddle_scission = e_saddle_scission; 03807 } 03808 else { 03809 e_saddle_scission = 0.; 03810 } 03811 // /* Semiempirical fission model: */ 03812 03813 // /* Fit to experimental result on curvature of potential at saddle */ 03814 // /* reference: */ 03815 // /* IF Z**2/A < 33.15E0 THEN 03816 // MassCurv = 30.5438538E0 - 4.00212049E0 * Z**2/A 03817 // + 0.11983384E0 * Z**4 / (A**2) ; 03818 // ELSE 03819 // MassCurv = 10.E0 ** (7.16993332E0 - 0.26602401E0 * Z**2/A 03820 // + 0.00283802E0 * Z**4 / (A**2)) ; */ 03821 // /* New parametrization of T. Enqvist according to Mulgin et al. 1998 */ 03822 if ( (std::pow(z,2))/a < 34.0) { 03823 masscurv = std::pow( 10.0,(-1.093364 + 0.082933 * (std::pow(z,2)/a) 03824 - 0.0002602 * (std::pow(z,4)/std::pow(a,2))) ); 03825 } else { 03826 masscurv = std::pow( 10.0,(3.053536 - 0.056477 * (std::pow(z,2)/a) 03827 + 0.0002454 * (std::pow(z,4)/std::pow(a,2))) ); 03828 } 03829 03830 cz_symm = (8.0/std::pow(z,2)) * masscurv; 03831 03832 if(itest == 1) { 03833 G4cout << "cz_symmetry= " << cz_symm << G4endl; 03834 } 03835 03836 if (cz_symm < 0) { 03837 icz = -1; 03838 // GOTO 1002; 03839 goto milledeux; 03840 } 03841 03842 // /* proton number in symmetric fission (centre) */ 03843 zsymm = z/2.0; 03844 nsymm = n/2.0; 03845 asymm = nsymm + zsymm; 03846 03847 zheavy1 = (cz_symm*zsymm + cz_asymm1_shell*zheavy1_shell)/(cz_symm + cz_asymm1_shell); 03848 zheavy2 = (cz_symm*zsymm + cz_asymm2_shell*zheavy2_shell)/(cz_symm + cz_asymm2_shell); 03849 // /* position of valley due to influence of liquid-drop potential */ 03850 nheavy1_eff = (zheavy1 + (a/n * cpol1))*(n/z); 03851 nheavy2_eff = (zheavy2 + (a/n * cpol2))*(n/z); 03852 nlight1_eff = n - nheavy1_eff; 03853 nlight2_eff = n - nheavy2_eff; 03854 // /* proton number of light fragments (centre) */ 03855 zlight1 = z - zheavy1; 03856 // /* proton number of light fragments (centre) */ 03857 zlight2 = z - zheavy2; 03858 aheavy1 = nheavy1_eff + zheavy1; 03859 aheavy2 = nheavy2_eff + zheavy2; 03860 aheavy1_mean = aheavy1; 03861 aheavy2_mean = aheavy2; 03862 alight1 = nlight1_eff + zlight1; 03863 alight2 = nlight2_eff + zlight2; 03864 03865 a_levdens = a / xlevdens; 03866 a_levdens_heavy1 = aheavy1 / xlevdens; 03867 a_levdens_heavy2 = aheavy2 / xlevdens; 03868 a_levdens_light1 = alight1 / xlevdens; 03869 a_levdens_light2 = alight2 / xlevdens; 03870 gamma = a_levdens / (0.4 * (std::pow(a,1.3333)) ); 03871 gamma_heavy1 = ( a_levdens_heavy1 / (0.4 * (std::pow(aheavy1,1.3333)) ) ) * fgamma1; 03872 gamma_heavy2 = a_levdens_heavy2 / (0.4 * (std::pow(aheavy2,1.3333)) ); 03873 03874 cz_asymm1_saddle = cz_asymm1_shell + cz_symm; 03875 cz_asymm2_saddle = cz_asymm2_shell + cz_symm; 03876 03877 // Up to here: Ok! Checked CS 10/10/05 03878 03879 cn = umass(zsymm,(nsymm+1.),0.0) + umass(zsymm,(nsymm-1.),0.0) 03880 + 1.44 * (std::pow(zsymm,2))/ 03881 ( (std::pow(r_null,2)) * 03882 ( std::pow((asymm+1.0),0.33333) + std::pow((asymm-1.0),0.33333) ) * 03883 ( std::pow((asymm+1.0),0.33333) + std::pow((asymm-1.0),0.33333) ) ) 03884 - 2.0 * umass(zsymm,nsymm,0.0) 03885 - 1.44 * (std::pow(zsymm,2))/ 03886 ( ( 2.0 * r_null * (std::pow(asymm,0.33333)) ) * 03887 ( 2.0 * r_null * (std::pow(asymm,0.33333)) ) ); 03888 03889 // /* shell effect in valley of mode 1 */ 03890 delta_u1 = delta_u1_shell + (std::pow((zheavy1_shell-zheavy1),2))*cz_asymm1_shell; 03891 // /* shell effect in valley of mode 2 */ 03892 delta_u2 = delta_u2_shell + (std::pow((zheavy2_shell-zheavy2),2))*cz_asymm2_shell; 03893 03894 // /* liquid drop energies 03895 // at the centres of the different shell effects 03896 // with respect to liquid drop at symmetry: */ 03897 epot0_mode1_saddle = (std::pow((zheavy1-zsymm),2)) * cz_symm; 03898 epot0_mode2_saddle = (std::pow((zheavy2-zsymm),2)) * cz_symm; 03899 epot0_symm_saddle = 0.0; 03900 03901 if (itest == 1) { 03902 G4cout << "check zheavy1 = " << zheavy1 << G4endl; 03903 G4cout << "check zheavy2 = " << zheavy2 << G4endl; 03904 G4cout << "check zsymm = " << zsymm << G4endl; 03905 G4cout << "check czsymm = " << cz_symm << G4endl; 03906 G4cout << "check epot0_mode1_saddle = " << epot0_mode1_saddle << G4endl; 03907 G4cout << "check epot0_mode2_saddle = " << epot0_mode2_saddle << G4endl; 03908 G4cout << "check epot0_symm_saddle = " << epot0_symm_saddle << G4endl; 03909 G4cout << "delta_u1 = " << delta_u1 << G4endl; 03910 G4cout << "delta_u2 = " << delta_u2 << G4endl; 03911 } 03912 03913 // /* energies including shell effects 03914 // at the centres of the different shell effects 03915 // with respect to liquid drop at symmetry: */ 03916 epot_mode1_saddle = epot0_mode1_saddle + delta_u1; 03917 epot_mode2_saddle = epot0_mode2_saddle + delta_u2; 03918 epot_symm_saddle = epot0_symm_saddle; 03919 if (itest == 1) { 03920 G4cout << "check epot_mode1_saddle = " << epot_mode1_saddle << G4endl; 03921 G4cout << "check epot_mode2_saddle = " << epot_mode2_saddle << G4endl; 03922 G4cout << "check epot_symm_saddle = " << epot_symm_saddle << G4endl; 03923 } 03924 03925 // /* Minimum of potential with respect to ld potential at symmetry */ 03926 dueff = min(epot_mode1_saddle,epot_mode2_saddle); 03927 dueff = min(dueff,epot_symm_saddle); 03928 dueff = dueff - epot_symm_saddle; 03929 03930 eld = e + dueff + e_zero_point; 03931 03932 if (itest == 1) { 03933 G4cout << "check dueff = " << dueff << G4endl; 03934 G4cout << "check e = " << e << G4endl; 03935 G4cout << "check e_zero_point = " << e_zero_point << G4endl; 03936 G4cout << "check eld = " << eld << G4endl; 03937 } 03938 // Up to here: Ok! Checked CS 10/10/05 03939 03940 // /* E = energy above lowest effective barrier */ 03941 // /* Eld = energy above liquid-drop barrier */ 03942 03943 // /* Due to this treatment the energy E on input means the excitation */ 03944 // /* energy above the lowest saddle. */ 03945 03946 // /* These energies are not used */ 03947 eheavy1 = e * aheavy1 / a; 03948 eheavy2 = e * aheavy2 / a; 03949 elight1 = e * alight1 / a; 03950 elight2 = e * alight2 / a; 03951 03952 epsilon0_1_saddle = eld - e_zero_point - epot0_mode1_saddle; 03953 // /* excitation energy at saddle mode 1 without shell effect */ 03954 epsilon0_2_saddle = eld - e_zero_point - epot0_mode2_saddle; 03955 // /* excitation energy at saddle mode 2 without shell effect */ 03956 03957 epsilon_1_saddle = eld - e_zero_point - epot_mode1_saddle; 03958 // /* excitation energy at saddle mode 1 with shell effect */ 03959 epsilon_2_saddle = eld - e_zero_point - epot_mode2_saddle; 03960 // /* excitation energy at saddle mode 2 with shell effect */ 03961 epsilon_symm_saddle = eld - e_zero_point - epot_symm_saddle; 03962 03963 // /* global parameters */ 03964 eexc1_saddle = epsilon_1_saddle; 03965 eexc2_saddle = epsilon_2_saddle; 03966 eexc_max = max(eexc1_saddle,eexc2_saddle); 03967 eexc_max = max(eexc_max,eld); 03968 03969 // /* EEXC_MAX is energy above the lowest saddle */ 03970 03971 03972 epsilon0_1_scission = eld + e_saddle_scission - epot0_mode1_saddle; 03973 // /* excitation energy without shell effect */ 03974 epsilon0_2_scission = eld + e_saddle_scission - epot0_mode2_saddle; 03975 // /* excitation energy without shell effect */ 03976 03977 epsilon_1_scission = eld + e_saddle_scission - epot_mode1_saddle; 03978 // /* excitation energy at scission */ 03979 epsilon_2_scission = eld+ e_saddle_scission - epot_mode2_saddle; 03980 // /* excitation energy at scission */ 03981 epsilon_symm_scission = eld + e_saddle_scission - epot_symm_saddle; 03982 // /* excitation energy of symmetric fragment at scission */ 03983 03984 // /* Calculate widhts at the saddle: */ 03985 03986 e_eff1_saddle = epsilon0_1_saddle - delta_u1 * (std::exp((-epsilon_1_saddle*gamma))); 03987 03988 if (e_eff1_saddle > 0.0) { 03989 wzasymm1_saddle = std::sqrt( (0.5 * 03990 (std::sqrt(1.0/a_levdens*e_eff1_saddle)) / 03991 (cz_asymm1_shell * std::exp((-epsilon_1_saddle*gamma)) + cz_symm) ) ); 03992 } 03993 else { 03994 wzasymm1_saddle = 1.0; 03995 } 03996 03997 e_eff2_saddle = epsilon0_2_saddle - delta_u2 * (std::exp((-epsilon_2_saddle*gamma))); 03998 if (e_eff2_saddle > 0.0) { 03999 wzasymm2_saddle = std::sqrt( (0.5 * 04000 (std::sqrt(1.0/a_levdens*e_eff2_saddle)) / 04001 (cz_asymm2_shell * std::exp((-epsilon_2_saddle*gamma)) + cz_symm) ) ); 04002 } 04003 else { 04004 wzasymm2_saddle = 1.0; 04005 } 04006 04007 if (eld > e_zero_point) { 04008 if ( (eld + epsilon_symm_saddle) < 0.0) { 04009 G4cout << "<e> eld + epsilon_symm_saddle < 0" << G4endl; 04010 } 04011 wzsymm_saddle = std::sqrt( (0.5 * 04012 (std::sqrt(1.0/a_levdens*(eld+epsilon_symm_saddle))) / cz_symm ) ); 04013 } else { 04014 wzsymm_saddle = 1.0; 04015 } 04016 04017 if (itest == 1) { 04018 G4cout << "wz1(saddle) = " << wzasymm1_saddle << G4endl; 04019 G4cout << "wz2(saddle) = " << wzasymm2_saddle << G4endl; 04020 G4cout << "wzsymm(saddle) = " << wzsymm_saddle << G4endl; 04021 } 04022 04023 // /* Calculate widhts at the scission point: */ 04024 // /* fits of ref. Beizin 1991 (Plots brought to GSI by Sergei Zhdanov) */ 04025 04026 wzsymm_scission = wzsymm_saddle; 04027 04028 if (e_saddle_scission == 0.0) { 04029 04030 wzasymm1_scission = wzasymm1_saddle; 04031 wzasymm2_scission = wzasymm2_saddle; 04032 04033 } 04034 else { 04035 04036 if (nheavy1_eff > 75.0) { 04037 wzasymm1_scission = (std::sqrt(21.0)) * z/a; 04038 wzasymm2_scission = (std::sqrt (max( (70.0-28.0)/3.0*(z*z/a-35.0)+28.,0.0 )) ) * z/a; 04039 } 04040 else { 04041 wzasymm1_scission = wzasymm1_saddle; 04042 wzasymm2_scission = wzasymm2_saddle; 04043 } 04044 04045 } 04046 04047 wzasymm1_scission = max(wzasymm1_scission,wzasymm1_saddle); 04048 wzasymm2_scission = max(wzasymm2_scission,wzasymm2_saddle); 04049 04050 wzasymm1 = wzasymm1_scission * fwidth_asymm1; 04051 wzasymm2 = wzasymm2_scission * fwidth_asymm2; 04052 wzsymm = wzsymm_scission; 04053 04054 /* if (ITEST == 1) { 04055 G4cout << "WZ1(scission) = " << WZasymm1_scission << G4endl; 04056 G4cout << "WZ2(scission) = " << WZasymm2_scission << G4endl; 04057 G4cout << "WZsymm(scission) = " << WZsymm_scission << G4endl; 04058 } 04059 if (ITEST == 1) { 04060 G4cout << "WZ1(scission) final= " << WZasymm1 << G4endl; 04061 G4cout << "WZ2(scission) final= " << WZasymm2 << G4endl; 04062 G4cout << "WZsymm(scission) final= " << WZsymm << G4endl; 04063 } */ 04064 04065 wasymm = wzsymm * a/z; 04066 waheavy1 = wzasymm1 * a/z; 04067 waheavy2 = wzasymm2 * a/z; 04068 04069 wasymm_saddle = wzsymm_saddle * a/z; 04070 waheavy1_saddle = wzasymm1_saddle * a/z; 04071 waheavy2_saddle = wzasymm2_saddle * a/z; 04072 04073 if (itest == 1) { 04074 G4cout << "wasymm = " << wzsymm << G4endl; 04075 G4cout << "waheavy1 = " << waheavy1 << G4endl; 04076 G4cout << "waheavy2 = " << waheavy2 << G4endl; 04077 } 04078 // Up to here: Ok! Checked CS 11/10/05 04079 04080 if ( (epsilon0_1_saddle - delta_u1*std::exp((-epsilon_1_saddle*gamma_heavy1))) < 0.0) { 04081 sasymm1 = -10.0; 04082 } 04083 else { 04084 sasymm1 = 2.0 * std::sqrt( a_levdens * (epsilon0_1_saddle - 04085 delta_u1*(std::exp((-epsilon_1_saddle*gamma_heavy1))) ) ); 04086 } 04087 04088 if ( (epsilon0_2_saddle - delta_u2*std::exp((-epsilon_2_saddle*gamma_heavy2))) < 0.0) { 04089 sasymm2 = -10.0; 04090 } 04091 else { 04092 sasymm2 = 2.0 * std::sqrt( a_levdens * (epsilon0_2_saddle - 04093 delta_u2*(std::exp((-epsilon_2_saddle*gamma_heavy2))) ) ); 04094 } 04095 04096 if (epsilon_symm_saddle > 0.0) { 04097 ssymm = 2.0 * std::sqrt( a_levdens*(epsilon_symm_saddle) ); 04098 } 04099 else { 04100 ssymm = -10.0; 04101 } 04102 04103 if (ssymm > -10.0) { 04104 ysymm = 1.0; 04105 04106 if (epsilon0_1_saddle < 0.0) { 04107 // /* low energy */ 04108 yasymm1 = std::exp((sasymm1-ssymm)) * wzasymm1_saddle / wzsymm_saddle * 2.0; 04109 // /* factor of 2 for symmetry classes */ 04110 } 04111 else { 04112 // /* high energy */ 04113 ssymm_mode1 = 2.0 * std::sqrt( a_levdens*(epsilon0_1_saddle) ); 04114 yasymm1 = ( std::exp((sasymm1-ssymm)) - std::exp((ssymm_mode1 - ssymm)) ) 04115 * wzasymm1_saddle / wzsymm_saddle * 2.0; 04116 } 04117 04118 if (epsilon0_2_saddle < 0.0) { 04119 // /* low energy */ 04120 yasymm2 = std::exp((sasymm2-ssymm)) * wzasymm2_saddle / wzsymm_saddle * 2.0; 04121 // /* factor of 2 for symmetry classes */ 04122 } 04123 else { 04124 // /* high energy */ 04125 ssymm_mode2 = 2.0 * std::sqrt( a_levdens*(epsilon0_2_saddle) ); 04126 yasymm2 = ( std::exp((sasymm2-ssymm)) - std::exp((ssymm_mode2 - ssymm)) ) 04127 * wzasymm2_saddle / wzsymm_saddle * 2.0; 04128 } 04129 // /* difference in the exponent in order */ 04130 // /* to avoid numerical overflow */ 04131 04132 } 04133 else { 04134 if ( (sasymm1 > -10.0) && (sasymm2 > -10.0) ) { 04135 ysymm = 0.0; 04136 yasymm1 = std::exp(sasymm1) * wzasymm1_saddle * 2.0; 04137 yasymm2 = std::exp(sasymm2) * wzasymm2_saddle * 2.0; 04138 } 04139 } 04140 04141 // /* normalize */ 04142 ysum = ysymm + yasymm1 + yasymm2; 04143 if (ysum > 0.0) { 04144 ysymm = ysymm / ysum; 04145 yasymm1 = yasymm1 / ysum; 04146 yasymm2 = yasymm2 / ysum; 04147 yasymm = yasymm1 + yasymm2; 04148 } 04149 else { 04150 ysymm = 0.0; 04151 yasymm1 = 0.0; 04152 yasymm2 = 0.0; 04153 // /* search minimum threshold and attribute all events to this mode */ 04154 if ( (epsilon_symm_saddle < epsilon_1_saddle) && (epsilon_symm_saddle < epsilon_2_saddle) ) { 04155 ysymm = 1.0; 04156 } 04157 else { 04158 if (epsilon_1_saddle < epsilon_2_saddle) { 04159 yasymm1 = 1.0; 04160 } 04161 else { 04162 yasymm2 = 1.0; 04163 } 04164 } 04165 } 04166 04167 if (itest == 1) { 04168 G4cout << "ysymm normalized= " << ysymm << G4endl; 04169 G4cout << "yasymm1 normalized= " << yasymm1 << G4endl; 04170 G4cout << "yasymm2 normalized= " << yasymm2 << G4endl; 04171 } 04172 // Up to here: Ok! Ckecked CS 11/10/05 04173 04174 // /* even-odd effect */ 04175 // /* simple parametrization KHS, Nov. 2000. From Rejmund et al. */ 04176 if ((int)(z) % 2 == 0) { 04177 r_e_o_exp = -0.017 * (e_saddle_scission + eld) * (e_saddle_scission + eld); 04178 if ( r_e_o_exp < -307.0) { 04179 r_e_o_exp = -307.0; 04180 r_e_o = std::pow(10.0,r_e_o_exp); 04181 } 04182 else { 04183 r_e_o = std::pow(10.0,r_e_o_exp); 04184 } 04185 } 04186 else { 04187 r_e_o = 0.0; 04188 } 04189 04190 // $LOOP; /* event loop */ 04191 // I_COUNT = I_COUNT + 1; 04192 04193 // /* random decision: symmetric or asymmetric */ 04194 // /* IMODE = 1 means asymmetric fission, mode 1, 04195 // IMODE = 2 means asymmetric fission, mode 2, 04196 // IMODE = 3 means symmetric */ 04197 // RMODE = dble(HAZ(k)); 04198 // rmode = rnd.rndm(); 04199 04200 // Safety check added to make sure we always select well defined 04201 // fission mode. 04202 do { 04203 rmode = haz(k); 04204 // Cast for test CS 11/10/05 04205 // RMODE = 0.54; 04206 // rmode = 0.54; 04207 if (rmode < yasymm1) { 04208 imode = 1; 04209 } else if ( (rmode > yasymm1) && (rmode < (yasymm1+yasymm2)) ) { 04210 imode = 2; 04211 } else if ( (rmode > yasymm1) && (rmode > (yasymm1+yasymm2)) ) { 04212 imode = 3; 04213 } 04214 } while(imode == 0); 04215 04216 // /* determine parameters of the Z distribution */ 04217 // force imode (for testing, PK) 04218 // imode = 3; 04219 if (imode == 1) { 04220 z1mean = zheavy1; 04221 z1width = wzasymm1; 04222 } 04223 if (imode == 2) { 04224 z1mean = zheavy2; 04225 z1width = wzasymm2; 04226 } 04227 if (imode == 3) { 04228 z1mean = zsymm; 04229 z1width = wzsymm; 04230 } 04231 04232 if (itest == 1) { 04233 G4cout << "nbre aleatoire tire " << rmode << G4endl; 04234 G4cout << "fission mode " << imode << G4endl; 04235 G4cout << "z1mean= " << z1mean << G4endl; 04236 G4cout << "z1width= " << z1width << G4endl; 04237 } 04238 04239 // /* random decision: Z1 and Z2 at scission: */ 04240 z1 = 1.0; 04241 z2 = 1.0; 04242 while ( (z1<5.0) || (z2<5.0) ) { 04243 // Z1 = dble(GAUSSHAZ(K,sngl(Z1mean),sngl(Z1width))); 04244 // z1 = rnd.gaus(z1mean,z1width); 04245 z1 = gausshaz(k, z1mean, z1width); 04246 z2 = z - z1; 04247 } 04248 if (itest == 1) { 04249 G4cout << "ff charge sample " << G4endl; 04250 G4cout << "z1 = " << z1 << G4endl; 04251 G4cout << "z2 = " << z2 << G4endl; 04252 } 04253 04254 // CALL EVEN_ODD(Z1,R_E_O,I_HELP); 04255 // /* Integer proton number with even-odd effect */ 04256 // Z1 = REAL(I_HELP) 04257 // /* Z1 = INT(Z1+0.5E0); */ 04258 z2 = z - z1; 04259 04260 // /* average N of both fragments: */ 04261 if (imode == 1) { 04262 n1mean = (z1 + cpol1 * a/n) * n/z; 04263 } 04264 if (imode == 2) { 04265 n1mean = (z1 + cpol2 * a/n) * n/z; 04266 } 04267 /* CASE(99) ! only for testing; 04268 N1UCD = Z1 * N/Z; 04269 N2UCD = Z2 * N/Z; 04270 re1 = UMASS(Z1,N1UCD,0.6) +; 04271 & UMASS(Z2,N2UCD,0.6) +; 04272 & ECOUL(Z1,N1UCD,0.6,Z2,N2UCD,0.6,d); 04273 re2 = UMASS(Z1,N1UCD+1.,0.6) +; 04274 & UMASS(Z2,N2UCD-1.,0.6) +; 04275 & ECOUL(Z1,N1UCD+1.,0.6,Z2,N2UCD-1.,0.6,d); 04276 re3 = UMASS(Z1,N1UCD+2.,0.6) +; 04277 & UMASS(Z2,N2UCD-2.,0.6) +; 04278 & ECOUL(Z1,N1UCD+2.,0.6,Z2,N2UCD-2.,0.6,d); 04279 eps2 = (re1-2.0*re2+re3) / 2.0; 04280 eps1 = re2 - re1 - eps2; 04281 DN1_POL = - eps1 / (2.0 * eps2); 04282 N1mean = N1UCD + DN1_POL; */ 04283 if (imode == 3) { 04284 n1ucd = z1 * n/z; 04285 n2ucd = z2 * n/z; 04286 re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) + ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,d); 04287 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); 04288 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); 04289 eps2 = (re1-2.0*re2+re3) / 2.0; 04290 eps1 = re2 - re1 - eps2; 04291 dn1_pol = - eps1 / (2.0 * eps2); 04292 n1mean = n1ucd + dn1_pol; 04293 } 04294 // all fission modes features have been checked CS 11/10/05 04295 n2mean = n - n1mean; 04296 z2mean = z - z1mean; 04297 04298 // /* Excitation energies */ 04299 // /* formulated in energies in close consistency with the fission model */ 04300 04301 // /* E_defo = UMASS(Z*0.5E0,N*0.5E0,0.6E0) - 04302 // UMASS(Z*0.5E0,N*0.5E0,0); */ 04303 // /* calculates the deformation energy of the liquid drop for 04304 // deformation beta = 0.6 which is most probable at scission */ 04305 04306 // /* N1R and N2R provisionaly taken without fluctuations in 04307 // polarisation: */ 04308 n1r = n1mean; 04309 n2r = n2mean; 04310 a1r = n1r + z1; 04311 a2r = n2r + z2; 04312 04313 if (imode == 1) { /* N = 82 */; 04315 e_scission = max(epsilon_1_scission,1.0); 04316 if (n1mean > (n * 0.5) ) { 04318 beta1 = 0.0; 04319 beta2 = 0.6; 04320 e1exc = epsilon_1_scission * a1r / a; 04321 e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0); 04322 e2exc = epsilon_1_scission * a2r / a + e_defo; 04323 } 04324 else { 04326 beta1 = 0.6; 04327 beta2 = 0.0; 04328 e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0); 04329 e1exc = epsilon_1_scission * a1r / a + e_defo; 04330 e2exc = epsilon_1_scission * a2r / a; 04331 } 04332 } 04333 04334 if (imode == 2) { 04336 e_scission = max(epsilon_2_scission,1.0); 04337 if (n1mean > (n * 0.5) ) { 04339 beta1 = (n1r - nheavy2) * 0.034 + 0.3; 04340 e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0); 04341 e1exc = epsilon_2_scission * a1r / a + e_defo; 04342 beta2 = 0.6 - beta1; 04343 e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0); 04344 e2exc = epsilon_2_scission * a2r / a + e_defo; 04345 } 04346 else { 04348 beta2 = (n2r - nheavy2) * 0.034 + 0.3; 04349 e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0); 04350 e2exc = epsilon_2_scission * a2r / a + e_defo; 04351 beta1 = 0.6 - beta2; 04352 e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0); 04353 e1exc = epsilon_2_scission * a1r / a + e_defo; 04354 } 04355 } 04356 04357 if (imode == 3) { 04358 // ! /* Symmetric fission channel */ 04359 04360 // /* the fit function for beta is the deformation for 04361 // optimum energy at the scission point, d = 2 */ 04362 // /* beta : deformation of symmetric fragments */ 04363 // /* beta1 : deformation of first fragment */ 04364 // /* beta2 : deformation of second fragment */ 04365 beta = 0.177963 + 0.0153241 * zsymm - 0.000162037 * zsymm*zsymm; 04366 beta1 = 0.177963 + 0.0153241 * z1 - 0.000162037 * z1*z1; 04367 // beta1 = 0.6 04368 e_defo1 = umass(z1,n1r,beta1) - umass(z1,n1r,0.0); 04369 beta2 = 0.177963 + 0.0153241 * z2 - 0.000162037 * z2*z2; 04370 // beta2 = 0.6 04371 e_defo2 = umass(z2,n2r,beta2) - umass(z2,n2r,0.0); 04372 e_asym = umass(z1 , n1r, beta1) + umass(z2, n2r ,beta2) 04373 + ecoul(z1,n1r,beta1,z2,n2r,beta2,2.0) 04374 - 2.0 * umass(zsymm,nsymm,beta) 04375 - ecoul(zsymm,nsymm,beta,zsymm,nsymm,beta,2.0); 04376 // E_asym = CZ_symm * (Z1 - Zsymm)**2 04377 e_scission = max((epsilon_symm_scission - e_asym),1.0); 04378 // /* $LIST(Z1,N1R,Z2,N2R,E_asym,E_scission); */ 04379 e1exc = e_scission * a1r / a + e_defo1; 04380 e2exc = e_scission * a2r / a + e_defo2; 04381 } 04382 // Energies checked for all the modes CS 11/10/05 04383 04384 // /* random decision: N1R and N2R at scission, before evaporation: */ 04385 // /* CN = UMASS(Zsymm , Nsymm + 1.E0,0) + 04386 // UMASS(Zsymm, Nsymm - 1.E0,0) 04387 // + 1.44E0 * (Zsymm)**2 / 04388 // (r_null**2 * ((Asymm+1)**1/3 + (Asymm-1)**1/3)**2 ) 04389 // - 2.E0 * UMASS(Zsymm,Nsymm,0) 04390 // - 1.44E0 * (Zsymm)**2 / (r_null * 2.E0 * (Asymm)**1/3)**2; */ 04391 04392 04393 // /* N1width = std::sqrt(0.5E0 * std::sqrt(1.E0/A_levdens*(Eld+E_saddle_scission)) / CN); */ 04394 // /* 8. 9. 1998: KHS (see also consideration in the first comment block) 04395 // sigma_N(Z=const) = A/Z * sigma_Z(A=const) 04396 // sigma_Z(A=const) = 0.4 to 0.5 (from Lang paper Nucl Phys. A345 (1980) 34) 04397 // sigma_N(Z=const) = 0.45 * A/Z (= 1.16 for 238U) 04398 // therefore: SIGZMIN = 1.16 */ 04399 04400 if ( (imode == 1) || (imode == 2) ) { 04401 cn=(umass(z1,n1mean+1.,beta1) + umass(z1,n1mean-1.,beta1) 04402 + umass(z2,n2mean+1.,beta2) + umass(z2,n2mean-1.,beta2) 04403 + ecoul(z1,n1mean+1.,beta1,z2,n2mean-1.,beta2,2.0) 04404 + ecoul(z1,n1mean-1.,beta1,z2,n2mean+1.,beta2,2.0) 04405 - 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0) 04406 - 2.0 * umass(z1, n1mean, beta1) 04407 - 2.0 * umass(z2, n2mean, beta2) ) * 0.5; 04408 // /* Coulomb energy neglected for the moment! */ 04409 // IF (E_scission.lt.0.) Then 04410 // write(6,*)'<E> E_scission < 0, MODE 1,2' 04411 // ENDIF 04412 // IF (CN.lt.0.) Then 04413 // write(6,*)'CN < 0, MODE 1,2' 04414 // ENDIF 04415 n1width=std::sqrt( (0.5 * (std::sqrt(1.0/a_levdens*(e_scission)))/cn) ); 04416 n1width=max(n1width, sigzmin); 04417 04418 // /* random decision: N1R and N2R at scission, before evaporation: */ 04419 n1r = 1.0; 04420 n2r = 1.0; 04421 while ( (n1r<5.0) || (n2r<5.0) ) { 04422 // n1r = dble(gausshaz(k,sngl(n1mean),sngl(n1width))); 04423 // n1r = rnd.gaus(n1mean,n1width); 04424 n1r = gausshaz(k, n1mean, n1width); 04425 n2r = n - n1r; 04426 } 04427 // N1R = GAUSSHAZ(K,N1mean,N1width) 04428 if (itest == 1) { 04429 G4cout << "after neutron sample " << n1r << G4endl; 04430 } 04431 n1r = (float)( (int)((n1r+0.5)) ); 04432 n2r = n - n1r; 04433 04434 even_odd(z1,r_e_o,i_help); 04435 // /* proton number with even-odd effect */ 04436 z1 = (float)(i_help); 04437 z2 = z - z1; 04438 04439 a1r = z1 + n1r; 04440 a2r = z2 + n2r; 04441 } 04442 04443 if (imode == 3) { 04445 if (nzpol > 0.0) { 04446 // /* We treat a simultaneous split in Z and N to determine polarisation */ 04447 cz = ( umass(z1-1., n1mean+1.,beta1) 04448 + umass(z2+1., n2mean-1.,beta1) 04449 + umass(z1+1., n1mean-1.,beta2) 04450 + umass(z2 - 1., n2mean + 1.,beta2) 04451 + ecoul(z1-1.,n1mean+1.,beta1,z2+1.,n2mean-1.,beta2,2.0) 04452 + ecoul(z1+1.,n1mean-1.,beta1,z2-1.,n2mean+1.,beta2,2.0) 04453 - 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0) 04454 - 2.0 * umass(z1, n1mean,beta1) 04455 - 2.0 * umass(z2, n2mean,beta2) ) * 0.5; 04456 // IF (E_scission.lt.0.) Then 04457 // write(6,*) '<E> E_scission < 0, MODE 1,2' 04458 // ENDIF 04459 // IF (CZ.lt.0.) Then 04460 // write(6,*) 'CZ < 0, MODE 1,2' 04461 // ENDIF 04462 za1width=std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*(e_scission)) / cz) ); 04463 za1width=std::sqrt( (max((za1width*za1width-(1.0/12.0)),0.1)) ); 04464 // /* Check the value of 0.1 ! */ 04465 // /* Shephard correction */ 04466 a1r = z1 + n1mean; 04467 a1r = (float)((int)((a1r+0.5))); 04468 a2r = a - a1r; 04469 // /* A1R and A2R are integer numbers now */ 04470 // /* $LIST(A1R,A2R,ZA1WIDTH); */ 04471 04472 n1ucd = n/a * a1r; 04473 n2ucd = n/a * a2r; 04474 z1ucd = z/a * a1r; 04475 z2ucd = z/a * a2r; 04476 04477 re1 = umass(z1ucd-1.,n1ucd+1.,beta1) + umass(z2ucd+1.,n2ucd-1.,beta2) 04478 + ecoul(z1ucd-1.,n1ucd+1.,beta1,z2ucd+1.,n2ucd-1.,beta2,d); 04479 re2 = umass(z1ucd,n1ucd,beta1) + umass(z2ucd,n2ucd,beta2) 04480 + ecoul(z1ucd,n1ucd,beta1,z2ucd,n2ucd,beta2,d); 04481 re3 = umass(z1ucd+1.,n1ucd-1.,beta1) + umass(z2ucd-1.,n2ucd+1.,beta2) + 04482 + ecoul(z1ucd+1.,n1ucd-1.,beta1,z2ucd-1.,n2ucd+1.,beta2,d); 04483 04484 eps2 = (re1-2.0*re2+re3) / 2.0; 04485 eps1 = (re3 - re1)/2.0; 04486 dn1_pol = - eps1 / (2.0 * eps2); 04487 z1 = z1ucd + dn1_pol; 04488 if (itest == 1) { 04489 G4cout << "before proton sample " << z1 << G4endl; 04490 } 04491 // Z1 = dble(GAUSSHAZ(k,sngl(Z1),sngl(ZA1width))); 04492 // z1 = rnd.gaus(z1,za1width); 04493 z1 = gausshaz(k, z1, za1width); 04494 if (itest == 1) { 04495 G4cout << "after proton sample " << z1 << G4endl; 04496 } 04497 even_odd(z1,r_e_o,i_help); 04498 // /* proton number with even-odd effect */ 04499 z1 = (float)(i_help); 04500 z2 = (float)((int)( (z - z1 + 0.5)) ); 04501 04502 n1r = a1r - z1; 04503 n2r = n - n1r; 04504 } 04505 else { 04506 // /* First division of protons, then adjustment of neutrons */ 04507 cn = ( umass(z1, n1mean+1.,beta1) + umass(z1, n1mean-1., beta1) 04508 + umass(z2, n2mean+1.,beta2) + umass(z2, n2mean-1., beta2) 04509 + ecoul(z1,n1mean+1.,beta1,z2,n2mean-1.,beta2,2.0) 04510 + ecoul(z1,n1mean-1.,beta1,z2,n2mean+1.,beta2,2.0) 04511 - 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0) 04512 - 2.0 * umass(z1, n1mean, 0.6) 04513 - 2.0 * umass(z2, n2mean, 0.6) ) * 0.5; 04514 // /* Coulomb energy neglected for the moment! */ 04515 // IF (E_scission.lt.0.) Then 04516 // write(6,*) '<E> E_scission < 0, MODE 1,2' 04517 // Endif 04518 // IF (CN.lt.0.) Then 04519 // write(6,*) 'CN < 0, MODE 1,2' 04520 // Endif 04521 n1width=std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*(e_scission)) / cn) ); 04522 n1width=max(n1width, sigzmin); 04523 04524 // /* random decision: N1R and N2R at scission, before evaporation: */ 04525 // N1R = dble(GAUSSHAZ(k,sngl(N1mean),sngl(N1width))); 04526 // n1r = rnd.gaus(n1mean,n1width); 04527 n1r = gausshaz(k, n1mean, n1width); 04528 n1r = (float)( (int)((n1r+0.5)) ); 04529 n2r = n - n1r; 04530 04531 even_odd(z1,r_e_o,i_help); 04532 // /* Integer proton number with even-odd effect */ 04533 z1 = (float)(i_help); 04534 z2 = z - z1; 04535 04536 a1r = z1 + n1r; 04537 a2r = z2 + n2r; 04538 04539 } 04540 } 04541 04542 if (itest == 1) { 04543 G4cout << "remid imode = " << imode << G4endl; 04544 G4cout << "n1width = " << n1width << G4endl; 04545 G4cout << "n1r = " << n1r << G4endl; 04546 G4cout << "a1r = " << a1r << G4endl; 04547 G4cout << "n2r = " << n2r << G4endl; 04548 G4cout << "a2r = " << a2r << G4endl; 04549 } 04550 // Up to here: checked CS 11/10/05 04551 04552 // /* Extracted from Lang et al. Nucl. Phys. A 345 (1980) 34 */ 04553 e1exc_sigma = 5.5; 04554 e2exc_sigma = 5.5; 04555 04556 neufcentquatrevingtsept: 04557 // E1final = dble(Gausshaz(k,sngl(E1exc),sngl(E1exc_sigma))); 04558 // E2final = dble(Gausshaz(k,sngl(E2exc),sngl(E2exc_sigma))); 04559 // e1final = rnd.gaus(e1exc,e1exc_sigma); 04560 // e2final = rnd.gaus(e2exc,e2exc_sigma); 04561 e1final = gausshaz(k, e1exc, e1exc_sigma); 04562 e2final = gausshaz(k, e2exc, e2exc_sigma); 04563 if ( (e1final < 0.0) || (e2final < 0.0) ) goto neufcentquatrevingtsept; 04564 if (itest == 1) { 04565 G4cout << "sampled exc 1 " << e1final << G4endl; 04566 G4cout << "sampled exc 2 " << e2final << G4endl; 04567 } 04568 04569 // /* OUTPUT QUANTITIES OF THE EVENT GENERATOR: */ 04570 04571 // /* Quantities before neutron evaporation */ 04572 04573 // /* Neutron number of prefragments: N1R and N2R */ 04574 // /* Atomic number of fragments: Z1 and Z2 */ 04575 // /* Kinetic energy of fragments: EkinR1, EkinR2 *7 04576 04577 // /* Quantities after neutron evaporation: */ 04578 04579 // /* Neutron number of fragments: N1 and N2 */ 04580 // /* Mass number of fragments: A1 and A2 */ 04581 // /* Atomic number of fragments: Z1 and Z2 */ 04582 // /* Number of evaporated neutrons: N1R-N1 and N2R-N2 */ 04583 // /* Kinetic energy of fragments: EkinR1*A1/A1R and 04584 // EkinR2*A2/A2R */ 04585 04586 n1 = n1r; 04587 n2 = n2r; 04588 a1 = n1 + z1; 04589 a2 = n2 + z2; 04590 e1 = e1final; 04591 e2 = e2final; 04592 04593 // /* Pre-neutron-emission total kinetic energy: */ 04594 tker = (z1 * z2 * 1.44) / 04595 ( r0 * std::pow(a1,0.33333) * (1.0 + 2.0/3.0 * beta1) + 04596 r0 * std::pow(a2,0.33333) * (1.0 + 2.0/3.0 * beta2) + 2.0 ); 04597 // /* Pre-neutron-emission kinetic energy of 1. fragment: */ 04598 ekinr1 = tker * a2 / a; 04599 // /* Pre-neutron-emission kinetic energy of 2. fragment: */ 04600 ekinr2 = tker * a1 / a; 04601 04602 v1 = std::sqrt( (ekinr1/a1) ) * 1.3887; 04603 v2 = std::sqrt( (ekinr2/a2) ) * 1.3887; 04604 04605 if (itest == 1) { 04606 G4cout << "ekinr1 " << ekinr1 << G4endl; 04607 G4cout << "ekinr2 " << ekinr2 << G4endl; 04608 } 04609 04610 milledeux: 04611 //************************** 04612 //*** only symmetric fission 04613 //************************** 04614 // Symmetric fission: Ok! Checked CS 10/10/05 04615 if ( (icz == -1) || (a1 < 0.0) || (a2 < 0.0) ) { 04616 // IF (z.eq.92) THEN 04617 // write(6,*)'symmetric fission' 04618 // write(6,*)'Z,A,E,A1,A2,icz,Atot',Z,A,E,A1,A2,icz,Atot 04619 // END IF 04620 04621 if (itest == 1) { 04622 G4cout << "milledeux: liquid-drop option " << G4endl; 04623 } 04624 04625 n = a-z; 04626 // proton number in symmetric fission (centre) * 04627 zsymm = z / 2.0; 04628 nsymm = n / 2.0; 04629 asymm = nsymm + zsymm; 04630 04631 a_levdens = a / xlevdens; 04632 04633 masscurv = 2.0; 04634 cz_symm = 8.0 / std::pow(z,2) * masscurv; 04635 04636 wzsymm = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cz_symm) ) ; 04637 04638 if (itest == 1) { 04639 G4cout << " symmetric high energy fission " << G4endl; 04640 G4cout << "wzsymm " << wzsymm << G4endl; 04641 } 04642 04643 z1mean = zsymm; 04644 z1width = wzsymm; 04645 04646 // random decision: Z1 and Z2 at scission: */ 04647 z1 = 1.0; 04648 z2 = 1.0; 04649 while ( (z1 < 5.0) || (z2 < 5.0) ) { 04650 // z1 = dble(gausshaz(kkk,sngl(z1mean),sngl(z1width))); 04651 // z1 = rnd.gaus(z1mean,z1width); 04652 z1 = gausshaz(kkk, z1mean, z1width); 04653 z2 = z - z1; 04654 } 04655 04656 if (itest == 1) { 04657 G4cout << " z1 " << z1 << G4endl; 04658 G4cout << " z2 " << z2 << G4endl; 04659 } 04660 if (itest == 1) { 04661 G4cout << " zsymm " << zsymm << G4endl; 04662 G4cout << " nsymm " << nsymm << G4endl; 04663 G4cout << " asymm " << asymm << G4endl; 04664 } 04665 // CN = UMASS(Zsymm , Nsymm + 1.E0) + UMASS(Zsymm, Nsymm - 1.E0) 04666 // # + 1.44E0 * (Zsymm)**2 / 04667 // # (r_null**2 * ((Asymm+1)**(1./3.) + 04668 // # (Asymm-1)**(1./3.))**2 ) 04669 // # - 2.E0 * UMASS(Zsymm,Nsymm) 04670 // # - 1.44E0 * (Zsymm)**2 / 04671 // # (r_null * 2.E0 * (Asymm)**(1./3.))**2 04672 04673 n1ucd = z1 * n/z; 04674 n2ucd = z2 * n/z; 04675 re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) + 04676 ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,2.0); 04677 re2 = umass(z1,n1ucd+1.,0.6) + umass(z2,n2ucd-1.,0.6) + 04678 ecoul(z1,n1ucd+1.,0.6,z2,n2ucd-1.,0.6,2.0); 04679 re3 = umass(z1,n1ucd+2.,0.6) + umass(z2,n2ucd-2.,0.6) + 04680 ecoul(z1,n1ucd+2.,0.6,z2,n2ucd-2.,0.6,2.0); 04681 reps2 = (re1-2.0*re2+re3)/2.0; 04682 reps1 = re2 - re1 -reps2; 04683 rn1_pol = -reps1/(2.0*reps2); 04684 n1mean = n1ucd + rn1_pol; 04685 n2mean = n - n1mean; 04686 04687 if (itest == 1) { 04688 G4cout << " n1mean " << n1mean << G4endl; 04689 G4cout << " n2mean " << n2mean << G4endl; 04690 } 04691 04692 cn = (umass(z1,n1mean+1.,0.0) + umass(z1,n1mean-1.,0.0) + 04693 + umass(z2,n2mean+1.,0.0) + umass(z2,n2mean-1.,0.0) 04694 - 2.0 * umass(z1,n1mean,0.0) + 04695 - 2.0 * umass(z2,n2mean,0.0) ) * 0.5; 04696 // This is an approximation! Coulomb energy is neglected. 04697 04698 n1width = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cn) ); 04699 04700 if (itest == 1) { 04701 G4cout << " cn " << cn << G4endl; 04702 G4cout << " n1width " << n1width << G4endl; 04703 } 04704 04705 // random decision: N1R and N2R at scission, before evaporation: */ 04706 // N1R = dfloat(NINT(GAUSSHAZ(KKK,sngl(N1mean),sngl(N1width)))); 04707 // n1r = (float)( (int)(rnd.gaus(n1mean,n1width)) ); 04708 n1r = (float)( (int)(gausshaz(k, n1mean,n1width)) ); 04709 n2r = n - n1r; 04710 // Mass of first and second fragment */ 04711 a1 = z1 + n1r; 04712 a2 = z2 + n2r; 04713 04714 e1 = e*a1/(a1+a2); 04715 e2 = e - e*a1/(a1+a2); 04716 if (itest == 1) { 04717 G4cout << " n1r " << n1r << G4endl; 04718 G4cout << " n2r " << n2r << G4endl; 04719 } 04720 04721 } 04722 04723 if (itest == 1) { 04724 G4cout << " a1 " << a1 << G4endl; 04725 G4cout << " z1 " << z1 << G4endl; 04726 G4cout << " a2 " << a2 << G4endl; 04727 G4cout << " z2 " << z2 << G4endl; 04728 G4cout << " e1 " << e1 << G4endl; 04729 G4cout << " e2 " << e << G4endl; 04730 } 04731 04732 // /* Pre-neutron-emission total kinetic energy: */ 04733 tker = (z1 * z2 * 1.44) / 04734 ( r0 * std::pow(a1,0.33333) * (1.0 + 2.0/3.0 * beta1) + 04735 r0 * std::pow(a2,0.33333) * (1.0 + 2.0/3.0 * beta2) + 2.0 ); 04736 // /* Pre-neutron-emission kinetic energy of 1. fragment: */ 04737 ekin1 = tker * a2 / a; 04738 // /* Pre-neutron-emission kinetic energy of 2. fragment: */ 04739 ekin2 = tker * a1 / a; 04740 04741 v1 = std::sqrt( (ekin1/a1) ) * 1.3887; 04742 v2 = std::sqrt( (ekin2/a2) ) * 1.3887; 04743 04744 if (itest == 1) { 04745 G4cout << " kinetic energies " << G4endl; 04746 G4cout << " ekin1 " << ekin1 << G4endl; 04747 G4cout << " ekin2 " << ekin2 << G4endl; 04748 } 04749 }
tirage aleatoire dans une maxwellienne
Definition at line 3220 of file G4Abla.cc.
References f(), fd(), nint(), and standardRandom().
Referenced by direct().
03221 { 03222 // tirage aleatoire dans une maxwellienne 03223 // t : temperature 03224 // 03225 // declaration des variables 03226 // 03227 03228 const int pSize = 101; 03229 static G4double p[pSize]; 03230 03231 // ial generateur pour le cascade (et les iy pour eviter les correlations) 03232 static G4int i = 0; 03233 static G4int itest = 0; 03234 // programme principal 03235 03236 // calcul des p(i) par approximation de newton 03237 p[pSize-1] = 8.0; 03238 G4double x = 0.1; 03239 G4double x1 = 0.0; 03240 G4double y = 0.0; 03241 03242 if (itest == 1) { 03243 goto fmaxhaz120; 03244 } 03245 03246 for(i = 1; i <= 99; i++) { 03247 fmaxhaz20: 03248 x1 = x - (f(x) - double(i)/100.0)/fd(x); 03249 x = x1; 03250 if (std::fabs(f(x) - double(i)/100.0) < 1e-5) { 03251 goto fmaxhaz100; 03252 } 03253 goto fmaxhaz20; 03254 fmaxhaz100: 03255 p[i] = x; 03256 } //end do 03257 03258 // itest = 1; 03259 itest = 0; 03260 // tirage aleatoire et calcul du x correspondant 03261 // par regression lineaire 03262 fmaxhaz120: 03263 standardRandom(&y, &(hazard->igraine[17])); 03264 i = nint(y*100); 03265 03266 // 2590 c ici on evite froidement les depassements de tableaux....(a.b. 3/9/99) 03267 if(i == 0) { 03268 goto fmaxhaz120; 03269 } 03270 03271 if (i == 1) { 03272 x = p[i]*y*100; 03273 } 03274 else { 03275 x = (p[i] - p[i-1])*(y*100 - i) + p[i]; 03276 } 03277 03278 return(x*T); 03279 }
| G4double G4Abla::gausshaz | ( | int | k, | |
| double | xmoy, | |||
| double | sig | |||
| ) |
Definition at line 5071 of file G4Abla.cc.
References haz().
Referenced by fissionDistri().
05072 { 05073 // Gaussian random numbers: 05074 05075 // 1005 C*** TIRAGE ALEATOIRE DANS UNE GAUSSIENNE DE LARGEUR SIG ET MOYENNE XMOY 05076 static G4int iset = 0; 05077 static G4double v1,v2,r,fac,gset,gausshaz; 05078 05079 if(iset == 0) { //then 05080 do { 05081 v1 = 2.0*haz(k) - 1.0; 05082 v2 = 2.0*haz(k) - 1.0; 05083 r = std::pow(v1,2) + std::pow(v2,2); 05084 } while(r >= 1); 05085 05086 fac = std::sqrt(-2.*std::log(r)/r); 05087 // assert(isnan(fac) == false); 05088 gset = v1*fac; 05089 gausshaz = v2*fac*sig+xmoy; 05090 iset = 1; 05091 } 05092 else { 05093 gausshaz=gset*sig+xmoy; 05094 iset=0; 05095 } 05096 05097 return gausshaz; 05098 }
Definition at line 3319 of file G4Abla.cc.
Referenced by pace2().
03320 { 03321 // TABLE DE MASSES ET FORMULE DE MASSE TIRE DU PAPIER DE BRACK-GUET 03322 // Gives the theoritical value for mass excess... 03323 // Révisée pour x, z flottants 25/4/2002 03324 03325 //real*8 x,z 03326 // dimension q(0:50,0:70) 03327 G4double x = (*x_par); 03328 G4double z = (*z_par); 03329 G4double find = (*find_par); 03330 03331 const G4int qrows = 50; 03332 const G4int qcols = 70; 03333 G4double q[qrows][qcols]; 03334 for(G4int init_i = 0; init_i < qrows; init_i++) { 03335 for(G4int init_j = 0; init_j < qcols; init_j++) { 03336 q[init_i][init_j] = 0.0; 03337 } 03338 } 03339 03340 G4int ix=G4int(std::floor(x+0.5)); 03341 G4int iz=G4int(std::floor(z+0.5)); 03342 G4double zz = iz; 03343 G4double xx = ix; 03344 find = 0.0; 03345 G4double avol = 15.776; 03346 G4double asur = -17.22; 03347 G4double ac = -10.24; 03348 G4double azer = 8.0; 03349 G4double xjj = -30.03; 03350 G4double qq = -35.4; 03351 G4double c1 = -0.737; 03352 G4double c2 = 1.28; 03353 03354 if(ix <= 7) { 03355 q[0][1]=939.50; 03356 q[1][1]=938.21; 03357 q[1][2]=1876.1; 03358 q[1][3]=2809.39; 03359 q[2][4]=3728.34; 03360 q[2][3]=2809.4; 03361 q[2][5]=4668.8; 03362 q[2][6]=5606.5; 03363 q[3][5]=4669.1; 03364 q[3][6]=5602.9; 03365 q[3][7]=6535.27; 03366 q[4][6]=5607.3; 03367 q[4][7]=6536.1; 03368 q[5][7]=6548.3; 03369 find=q[iz][ix]; 03370 } 03371 else { 03372 G4double xneu=xx-zz; 03373 G4double si=(xneu-zz)/xx; 03374 G4double x13=std::pow(xx,.333); 03375 G4double ee1=c1*zz*zz/x13; 03376 G4double ee2=c2*zz*zz/xx; 03377 G4double aux=1.+(9.*xjj/4./qq/x13); 03378 G4double ee3=xjj*xx*si*si/aux; 03379 G4double ee4=avol*xx+asur*(std::pow(xx,.666))+ac*x13+azer; 03380 G4double tota = ee1 + ee2 + ee3 + ee4; 03381 find = 939.55*xneu+938.77*zz - tota; 03382 } 03383 03384 (*x_par) = x; 03385 (*z_par) = z; 03386 (*find_par) = find; 03387 }
Definition at line 5019 of file G4Abla.cc.
References mod(), nint(), secnds(), and standardRandom().
Referenced by evapora(), expohaz(), fissionDistri(), and gausshaz().
05020 { 05021 const G4int pSize = 110; 05022 static G4double p[pSize]; 05023 static G4long ix = 0, i = 0; 05024 static G4double x = 0.0, y = 0.0, a = 0.0, haz = 0.0; 05025 // k =< -1 on initialise 05026 // k = -1 c'est reproductible 05027 // k < -1 || k > -1 ce n'est pas reproductible 05028 05029 // Zero is invalid random seed. Set proper value from our random seed collection: 05030 if(ix == 0) { 05031 ix = hazard->ial; 05032 } 05033 05034 if (k <= -1) { //then 05035 if(k == -1) { //then 05036 ix = 0; 05037 } 05038 else { 05039 x = 0.0; 05040 y = secnds(int(x)); 05041 ix = int(y * 100 + 43543000); 05042 if(mod(ix,2) == 0) { 05043 ix = ix + 1; 05044 } 05045 } 05046 05047 // Here we are using random number generator copied from INCL code 05048 // instead of the CERNLIB one! This causes difficulties for 05049 // automatic testing since the random number generators, and thus 05050 // the behavior of the routines in C++ and FORTRAN versions is no 05051 // longer exactly the same! 05052 standardRandom(&x, &ix); 05053 for(G4int i = 0; i < pSize; i++) { //do i=1,110 05054 standardRandom(&(p[i]), &ix); 05055 } 05056 standardRandom(&a, &ix); 05057 k = 0; 05058 } 05059 05060 i = nint(100*a)+1; 05061 haz = p[i]; 05062 standardRandom(&a, &ix); 05063 p[i] = a; 05064 05065 hazard->ial = ix; 05066 05067 return haz; 05068 }
Definition at line 5235 of file G4Abla.cc.
Referenced by densniv(), direct(), evapora(), mglms(), mglw(), parite(), and spdef().
05236 { 05237 G4double valueCeil = int(std::ceil(value)); 05238 G4double valueFloor = int(std::floor(value)); 05239 05240 if(std::fabs(value - valueCeil) < std::fabs(value - valueFloor)) { 05241 return int(valueCeil); 05242 } 05243 else { 05244 return int(valueFloor); 05245 } 05246 }
| void G4Abla::initEvapora | ( | ) |
Initialize ABLA evaporation code.
Definition at line 746 of file G4Abla.cc.
References G4Ald::ak, G4Fiss::akap, G4Ecld::alpha, G4Ald::as, G4Ald::av, G4Fiss::bet, G4Pace::dm, G4Ecld::ecfnz, G4Ecld::ecgnz, G4Ec2sub::ecnz, G4Opt::eefac, G4cout, G4endl, G4Exception(), G4Fiss::homega, G4Fiss::ifis, G4Fiss::koeff, CLHEP::detail::n, G4Ald::optafan, G4Opt::optcha, G4Fiss::optcol, G4Opt::optemd, G4Fiss::optles, G4Fiss::optshp, G4Fiss::optxfis, and G4Ecld::vgsld.
00747 { 00748 // 37 C PROJECTILE AND TARGET PARAMETERS + CROSS SECTIONS 00749 // 38 C COMMON /ABLAMAIN/ AP,ZP,AT,ZT,EAP,BETA,BMAXNUC,CRTOT,CRNUC, 00750 // 39 C R_0,R_P,R_T, IMAX,IRNDM,PI, 00751 // 40 C BFPRO,SNPRO,SPPRO,SHELL 00752 // 41 C 00753 // 42 C AP,ZP,AT,ZT - PROJECTILE AND TARGET MASSES 00754 // 43 C EAP,BETA - BEAM ENERGY PER NUCLEON, V/C 00755 // 44 C BMAXNUC - MAX. IMPACT PARAMETER FOR NUCL. REAC. 00756 // 45 C CRTOT,CRNUC - TOTAL AND NUCLEAR REACTION CROSS SECTION 00757 // 46 C R_0,R_P,R_T, - RADIUS PARAMETER, PROJECTILE+ TARGET RADII 00758 // 47 C IMAX,IRNDM,PI - MAXIMUM NUMBER OF EVENTS, DUMMY, 3.141... 00759 // 48 C BFPRO - FISSION BARRIER OF THE PROJECTILE 00760 // 49 C SNPRO - NEUTRON SEPARATION ENERGY OF THE PROJECTILE 00761 // 50 C SPPRO - PROTON " " " " " 00762 // 51 C SHELL - GROUND STATE SHELL CORRECTION 00763 // 52 C--------------------------------------------------------------------- 00764 // 53 C 00765 // 54 C ENERGIES WIDTHS AND CROSS SECTIONS FOR EM EXCITATION 00766 // 55 C COMMON /EMDPAR/ EGDR,EGQR,FWHMGDR,FWHMGQR,CREMDE1,CREMDE2, 00767 // 56 C AE1,BE1,CE1,AE2,BE2,CE2,SR1,SR2,XR 00768 // 57 C 00769 // 58 C EGDR,EGQR - MEAN ENERGY OF GDR AND GQR 00770 // 59 C FWHMGDR,FWHMGQR - FWHM OF GDR, GQR 00771 // 60 C CREMDE1,CREMDE2 - EM CROSS SECTION FOR E1 AND E2 00772 // 61 C AE1,BE1,CE1 - ARRAYS TO CALCULATE 00773 // 62 C AE2,BE2,CE2 - THE EXCITATION ENERGY AFTER E.M. EXC. 00774 // 63 C SR1,SR2,XR - WITH MONTE CARLO 00775 // 64 C--------------------------------------------------------------------- 00776 // 65 C 00777 // 66 C DEFORMATIONS AND G.S. SHELL EFFECTS 00778 // 67 C COMMON /ECLD/ ECGNZ,ECFNZ,VGSLD,ALPHA 00779 // 68 C 00780 // 69 C ECGNZ - GROUND STATE SHELL CORR. FRLDM FOR A SPHERICAL G.S. 00781 // 70 C ECFNZ - SHELL CORRECTION FOR THE SADDLE POINT (NOW: == 0) 00782 // 71 C VGSLD - DIFFERENCE BETWEEN DEFORMED G.S. AND LDM VALUE 00783 // 72 C ALPHA - ALPHA GROUND STATE DEFORMATION (THIS IS NOT BETA2!) 00784 // 73 C BETA2 = SQRT(5/(4PI)) * ALPHA 00785 // 74 C--------------------------------------------------------------------- 00786 // 75 C 00787 // 76 C ARRAYS FOR EXCITATION ENERGY BY STATISTICAL HOLE ENERY MODEL 00788 // 77 C COMMON /EENUC/ SHE, XHE 00789 // 78 C 00790 // 79 C SHE, XHE - ARRAYS TO CALCULATE THE EXC. ENERGY AFTER 00791 // 80 C ABRASION BY THE STATISTICAL HOLE ENERGY MODEL 00792 // 81 C--------------------------------------------------------------------- 00793 // 82 C 00794 // 83 C G.S. SHELL EFFECT 00795 // 84 C COMMON /EC2SUB/ ECNZ 00796 // 85 C 00797 // 86 C ECNZ G.S. SHELL EFFECT FOR THE MASSES (IDENTICAL TO ECGNZ) 00798 // 87 C--------------------------------------------------------------------- 00799 // 88 C 00800 // 89 C OPTIONS AND PARAMETERS FOR FISSION CHANNEL 00801 // 90 C COMMON /FISS/ AKAP,BET,HOMEGA,KOEFF,IFIS, 00802 // 91 C OPTSHP,OPTXFIS,OPTLES,OPTCOL 00803 // 92 C 00804 // 93 C AKAP - HBAR**2/(2* MN * R_0**2) = 10 MEV 00805 // 94 C BET - REDUCED NUCLEAR FRICTION COEFFICIENT IN (10**21 S**-1) 00806 // 95 C HOMEGA - CURVATURE OF THE FISSION BARRIER = 1 MEV 00807 // 96 C KOEFF - COEFFICIENT FOR THE LD FISSION BARRIER == 1.0 00808 // 97 C IFIS - 0/1 FISSION CHANNEL OFF/ON 00809 // 98 C OPTSHP - INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY 00810 // 99 C = 0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY 00811 // 100 C = 1 SHELL , NO PAIRING 00812 // 101 C = 2 PAIRING, NO SHELL 00813 // 102 C = 3 SHELL AND PAIRING 00814 // 103 C OPTCOL - 0/1 COLLECTIVE ENHANCEMENT SWITCHED ON/OFF 00815 // 104 C OPTXFIS- 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV 00816 // 105 C FISSILITY PARAMETER. 00817 // 106 C OPTLES - CONSTANT TEMPERATURE LEVEL DENSITY FOR A,Z > TH-224 00818 // 107 C OPTCOL - 0/1 COLLECTIVE ENHANCEMENT OFF/ON 00819 // 108 C--------------------------------------------------------------------- 00820 // 109 C 00821 // 110 C OPTIONS 00822 // 111 C COMMON /OPT/ OPTEMD,OPTCHA,EEFAC 00823 // 112 C 00824 // 113 C OPTEMD - 0/1 NO EMD / INCL. EMD 00825 // 114 C OPTCHA - 0/1 0 GDR / 1 HYPERGEOMETRICAL PREFRAGMENT-CHARGE-DIST. 00826 // 115 C *** RECOMMENDED IS OPTCHA = 1 *** 00827 // 116 C EEFAC - EXCITATION ENERGY FACTOR, 2.0 RECOMMENDED 00828 // 117 C--------------------------------------------------------------------- 00829 // 118 C 00830 // 119 C FISSION BARRIERS 00831 // 120 C COMMON /FB/ EFA 00832 // 121 C EFA - ARRAY OF FISSION BARRIERS 00833 // 122 C--------------------------------------------------------------------- 00834 // 123 C 00835 // 124 C p LEVEL DENSITY PARAMETERS 00836 // 125 C COMMON /ALD/ AV,AS,AK,OPTAFAN 00837 // 126 C AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE 00838 // 127 C LEVEL DENSITY PARAMETER 00839 // 128 C OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1 00840 // 129 C RECOMMENDED IS OPTAFAN = 0 00841 // 130 C--------------------------------------------------------------------- 00842 // 131 C ____________________________________________________________________ 00843 // 132 C / 00844 // 133 C / INITIALIZES PARAMETERS IN COMMON /ABRAMAIN/, /EMDPAR/, /ECLD/ ... 00845 // 134 C / PROJECTILE PARAMETERS, EMD PARAMETERS, SHELL CORRECTION TABLES. 00846 // 135 C / CALCULATES MAXIMUM IMPACT PARAMETER FOR NUCLEAR COLLISIONS AND 00847 // 136 C / TOTAL GEOMETRICAL CROSS SECTION + EMD CROSS SECTIONS 00848 // 137 C ____________________________________________________________________ 00849 // 138 C 00850 // 139 C 00851 // 201 C 00852 // 202 C---------- SET INPUT VALUES 00853 // 203 C 00854 // 204 C *** INPUT FROM UNIT 10 IN THE FOLLOWING SEQUENCE ! 00855 // 205 C AP1 = INTEGER ! 00856 // 206 C ZP1 = INTEGER ! 00857 // 207 C AT1 = INTEGER ! 00858 // 208 C ZT1 = INTEGER ! 00859 // 209 C EAP = REAL ! 00860 // 210 C IMAX = INTEGER ! 00861 // 211 C IFIS = INTEGER SWITCH FOR FISSION 00862 // 212 C OPTSHP = INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY 00863 // 213 C =0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY 00864 // 214 C =1 SHELL , NO PAIRING CORRECTION 00865 // 215 C =2 PAIRING, NO SHELL CORRECTION 00866 // 216 C =3 SHELL AND PAIRING CORRECTION IN MASSES AND ENERGY 00867 // 217 C OPTEMD =0,1 0 NO EMD, 1 INCL. EMD 00868 // 218 C ELECTROMAGNETIC DISSOZIATION IS CALCULATED AS WELL. 00869 // 219 C OPTCHA =0,1 0 GDR- , 1 HYPERGEOMETRICAL PREFRAGMENT-CHARGE-DIST. 00870 // 220 C RECOMMENDED IS OPTCHA=1 00871 // 221 C OPTCOL =0,1 COLLECTIVE ENHANCEMENT SWITCHED ON 1 OR OFF 0 IN DENSN 00872 // 222 C OPTAFAN=0,1 SWITCH FOR AF/AN = 1 IN DENSNIV 0 AF/AN>1 1 AF/AN=1 00873 // 223 C AKAP = REAL ALWAYS EQUALS 10 00874 // 224 C BET = REAL REDUCED FRICTION COEFFICIENT / 10**(+21) S**(-1) 00875 // 225 C HOMEGA = REAL CURVATURE / MEV RECOMMENDED = 1. MEV 00876 // 226 C KOEFF = REAL COEFFICIENT FOR FISSION BARRIER 00877 // 227 C OPTXFIS= INTEGER 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV 00878 // 228 C FISSILITY PARAMETER. 00879 // 229 C EEFAC = REAL EMPIRICAL FACTOR FOR THE EXCITATION ENERGY 00880 // 230 C RECOMMENDED 2.D0, STATISTICAL ABRASION MODELL 1.D0 00881 // 231 C AV = REAL KOEFFICIENTS FOR CALCULATION OF A(TILDE) 00882 // 232 C AS = REAL LEVEL DENSITY PARAMETER 00883 // 233 C AK = REAL 00884 // 234 C 00885 // 235 C This following inputs will be initialized in the main through the 00886 // 236 C common /ABLAMAIN/ (A.B.) 00887 // 237 00888 00889 // switch-fission.1=on.0=off 00890 fiss->ifis = 1; 00891 00892 // shell+pairing.0-1-2-3 00893 fiss->optshp = 0; 00894 00895 // optemd =0,1 0 no emd, 1 incl. emd 00896 opt->optemd = 1; 00897 // read(10,*,iostat=io) dum(10),optcha 00898 opt->optcha = 1; 00899 00900 // not.to.be.changed.(akap) 00901 fiss->akap = 10.0; 00902 00903 // nuclear.viscosity.(beta) 00904 fiss->bet = 1.5; 00905 00906 // potential-curvature 00907 fiss->homega = 1.0; 00908 00909 // fission-barrier-coefficient 00910 fiss->koeff = 1.; 00911 00912 //collective enhancement switched on 1 or off 0 in densn (qr=val or =1.) 00913 fiss->optcol = 0; 00914 00915 // switch-for-low-energy-sys 00916 fiss->optles = 0; 00917 00918 opt->eefac = 2.; 00919 00920 ald->optafan = 0; 00921 00922 ald->av = 0.073e0; 00923 ald->as = 0.095e0; 00924 ald->ak = 0.0e0; 00925 00926 if(verboseLevel > 3) { 00927 G4cout <<"ifis " << fiss->ifis << G4endl; 00928 G4cout <<"optshp " << fiss->optshp << G4endl; 00929 G4cout <<"optemd " << opt->optemd << G4endl; 00930 G4cout <<"optcha " << opt->optcha << G4endl; 00931 G4cout <<"akap " << fiss->akap << G4endl; 00932 G4cout <<"bet " << fiss->bet << G4endl; 00933 G4cout <<"homega " << fiss->homega << G4endl; 00934 G4cout <<"koeff " << fiss->koeff << G4endl; 00935 G4cout <<"optcol " << fiss->optcol << G4endl; 00936 G4cout <<"optles " << fiss->optles << G4endl; 00937 G4cout <<"eefac " << opt->eefac << G4endl; 00938 G4cout <<"optafan " << ald->optafan << G4endl; 00939 G4cout <<"av " << ald->av << G4endl; 00940 G4cout <<"as " << ald->as << G4endl; 00941 G4cout <<"ak " << ald->ak << G4endl; 00942 } 00943 fiss->optxfis = 1; 00944 00945 G4InclAblaDataFile *dataInterface = new G4InclAblaDataFile(); 00946 if(dataInterface->readData() == true) { 00947 if(verboseLevel > 0) { 00948 G4cout <<"G4Abla: Datafiles read successfully." << G4endl; 00949 } 00950 } 00951 else { 00952 G4Exception("ERROR: Failed to read datafiles."); 00953 } 00954 00955 for(int z = 0; z < 98; z++) { //do 30 z = 0,98,1 00956 for(int n = 0; n < 154; n++) { //do 31 n = 0,153,1 00957 ecld->ecfnz[n][z] = 0.e0; 00958 ec2sub->ecnz[n][z] = dataInterface->getEcnz(n,z); 00959 ecld->ecgnz[n][z] = dataInterface->getEcnz(n,z); 00960 ecld->alpha[n][z] = dataInterface->getAlpha(n,z); 00961 ecld->vgsld[n][z] = dataInterface->getVgsld(n,z); 00962 } 00963 } 00964 00965 for(int z = 0; z < 500; z++) { 00966 for(int a = 0; a < 500; a++) { 00967 pace->dm[z][a] = dataInterface->getPace2(z,a); 00968 } 00969 } 00970 00971 delete dataInterface; 00972 }
| void G4Abla::lorab | ( | G4double | gam, | |
| G4double | eta, | |||
| G4double | ein, | |||
| G4double | pin[], | |||
| G4double * | eout, | |||
| G4double | pout[] | |||
| ) |
Definition at line 4977 of file G4Abla.cc.
04979 { 04980 // C Transformation de lorentz brute pour vérifs. 04981 // C P(3) = P_longitudinal (transformé) 04982 // C P(1) et P(2) = P_transvers (non transformés) 04983 // DIMENSION Pin(3),Pout(3) 04984 // REAL*8 GAM,ETA,Ein 04985 04986 pout[1] = pin[1]; 04987 pout[2] = pin[2]; 04988 (*eout) = gam*ein + eta*pin[3]; 04989 pout[3] = eta*ein + gam*pin[3]; 04990 }
This subroutine calculates the ordinary legendre polynomials of order 0 to n-1 of argument x and stores them in the vector pl. They are calculated by recursion relation from the first two polynomials. Written by A.J.Sierk LANL t-9 February, 1984
Definition at line 2524 of file G4Abla.cc.
Referenced by barfit().
02525 { 02526 // THIS SUBROUTINE CALCULATES THE ORDINARY LEGENDRE POLYNOMIALS OF 02527 // ORDER 0 TO N-1 OF ARGUMENT X AND STORES THEM IN THE VECTOR PL. 02528 // THEY ARE CALCULATED BY RECURSION RELATION FROM THE FIRST TWO 02529 // POLYNOMIALS. 02530 // WRITTEN BY A.J.SIERK LANL T-9 FEBRUARY, 1984 02531 02532 // NOTE: PL AND X MUST BE DOUBLE PRECISION ON 32-BIT COMPUTERS! 02533 02534 pl[0] = 1.0; 02535 pl[1] = x; 02536 02537 for(int i = 2; i < n; i++) { 02538 pl[i] = ((2*double(i+1) - 3.0)*x*pl[i-1] - (double(i+1) - 2.0)*pl[i-2])/(double(i+1)-1.0); 02539 } 02540 }
Mglms
Definition at line 1051 of file G4Abla.cc.
References G4Ec2sub::ecnz, eflmac(), and idnint().
Referenced by breakItUp(), and direct().
01052 { 01053 // USING FUNCTION EFLMAC(IA,IZ,0) 01054 // 01055 // REFOPT4 = 0 : WITHOUT MICROSCOPIC CORRECTIONS 01056 // REFOPT4 = 1 : WITH SHELL CORRECTION 01057 // REFOPT4 = 2 : WITH PAIRING CORRECTION 01058 // REFOPT4 = 3 : WITH SHELL- AND PAIRING CORRECTION 01059 01060 // 1839 C----------------------------------------------------------------------- 01061 // 1840 C A1 LOCAL MASS NUMBER (INTEGER VARIABLE OF A) 01062 // 1841 C Z1 LOCAL NUCLEAR CHARGE (INTEGER VARIABLE OF Z) 01063 // 1842 C REFOPT4 OPTION, SPECIFYING THE MASS FORMULA (SEE ABOVE) 01064 // 1843 C A MASS NUMBER 01065 // 1844 C Z NUCLEAR CHARGE 01066 // 1845 C DEL PAIRING CORRECTION 01067 // 1846 C EL BINDING ENERGY 01068 // 1847 C ECNZ( , ) TABLE OF SHELL CORRECTIONS 01069 // 1848 C----------------------------------------------------------------------- 01070 // 1849 C 01071 G4int a1 = idnint(a); 01072 G4int z1 = idnint(z); 01073 01074 if ( (a1 <= 0) || (z1 <= 0) || ((a1-z1) <= 0) ) { //then 01075 // modif pour récupérer une masse p et n correcte: 01076 (*el) = 0.0; 01077 return; 01078 // goto mglms50; 01079 } 01080 else { 01081 // binding energy incl. pairing contr. is calculated from 01082 // function eflmac 01083 assert(a1 != 0); 01084 (*el) = eflmac(a1,z1,0,refopt4); 01085 // assert(isnan((*el)) == false); 01086 if (refopt4 > 0) { 01087 if (refopt4 != 2) { 01088 (*el) = (*el) + ec2sub->ecnz[a1-z1][z1]; 01089 //(*el) = (*el) + ec2sub->ecnz[z1][a1-z1]; 01090 //assert(isnan((*el)) == false); 01091 } 01092 } 01093 } 01094 return; 01095 }
Model de la goutte liquide de c. f. weizsacker. usually an obsolete option
Definition at line 1020 of file G4Abla.cc.
References idnint().
Referenced by direct().
01021 { 01022 // MODEL DE LA GOUTTE LIQUIDE DE C. F. WEIZSACKER. 01023 // USUALLY AN OBSOLETE OPTION 01024 01025 G4int a1 = 0, z1 = 0; 01026 G4double xv = 0.0, xs = 0.0, xc = 0.0, xa = 0.0; 01027 01028 a1 = idnint(a); 01029 z1 = idnint(z); 01030 01031 if ((a <= 0.01) || (z < 0.01)) { 01032 (*el) = 1.0e38; 01033 } 01034 else { 01035 xv = -15.56*a; 01036 xs = 17.23*std::pow(a,(2.0/3.0)); 01037 01038 if (a > 1.0) { 01039 xc = 0.7*z*(z-1.0)*std::pow((a-1.0),(-1.e0/3.e0)); 01040 } 01041 else { 01042 xc = 0.0; 01043 } 01044 } 01045 01046 xa = 23.6*(std::pow((a-2.0*z),2)/a); 01047 (*el) = xv+xs+xc+xa; 01048 return; 01049 }
Definition at line 5143 of file G4Abla.cc.
Referenced by fmaxhaz(), haz(), and translab().
05144 { 05145 G4double intpart = 0.0; 05146 G4double fractpart = 0.0; 05147 fractpart = std::modf(number, &intpart); 05148 if(number == 0) { 05149 return 0; 05150 } 05151 if(number > 0) { 05152 if(fractpart < 0.5) { 05153 return int(std::floor(number)); 05154 } 05155 else { 05156 return int(std::ceil(number)); 05157 } 05158 } 05159 if(number < 0) { 05160 if(fractpart < -0.5) { 05161 return int(std::floor(number)); 05162 } 05163 else { 05164 return int(std::ceil(number)); 05165 } 05166 } 05167 05168 return int(std::floor(number)); 05169 }
Definition at line 3281 of file G4Abla.cc.
References G4Pace::dm, guet(), and idint().
Referenced by breakItUp().
03282 { 03283 // PACE2 03284 // Cette fonction retourne le defaut de masse du noyau A,Z en MeV 03285 // Révisée pour a, z flottants 25/4/2002 = 03286 03287 G4double pace2 = 0.0; 03288 03289 G4int ii = idint(a+0.5); 03290 G4int jj = idint(z+0.5); 03291 03292 if(ii <= 0 || jj < 0) { 03293 pace2=0.; 03294 return pace2; 03295 } 03296 03297 if(jj > 300) { 03298 pace2=0.0; 03299 } 03300 else { 03301 pace2=pace->dm[ii][jj]; 03302 } 03303 pace2=pace2/1000.; 03304 03305 if(pace->dm[ii][jj] == 0.) { 03306 if(ii < 12) { 03307 pace2=-500.; 03308 } 03309 else { 03310 guet(&a, &z, &pace2); 03311 pace2=pace2-ii*931.5; 03312 pace2=pace2/1000.; 03313 } 03314 } 03315 03316 return pace2; 03317 }
PROCEDURE FOR CALCULATING THE PARITY OF THE NUMBER N. RETURNS -1 IF N IS ODD AND +1 IF N IS EVEN
Definition at line 2739 of file G4Abla.cc.
References dint(), and idnint().
Referenced by appariem(), densniv(), and direct().
02740 { 02741 // CALCUL DE LA PARITE DU NOMBRE N 02742 // 02743 // PROCEDURE FOR CALCULATING THE PARITY OF THE NUMBER N. 02744 // RETURNS -1 IF N IS ODD AND +1 IF N IS EVEN 02745 02746 G4double n1 = 0.0, n2 = 0.0, n3 = 0.0; 02747 02748 // N NUMBER TO BE TESTED 02749 // N1,N2 HELP VARIABLES 02750 // PAR HELP VARIABLE FOR PARITY OF N 02751 02752 n3 = double(idnint(n)); 02753 n1 = n3/2.0; 02754 n2 = n1 - dint(n1); 02755 02756 if (n2 > 0.0) { 02757 (*par) = -1.0; 02758 } 02759 else { 02760 (*par) = 1.0; 02761 } 02762 }
Coefficient of collective enhancement including damping Input: z,a,bet,sig,u Output: qr - collective enhancement factor See junghans et al., nucl. phys. a 629 (1998) 635
| z | charge number | |
| a | mass number | |
| bet | beta deformation | |
| sig | perpendicular spin cut-off factor | |
| u | Energy |
Definition at line 974 of file G4Abla.cc.
References CLHEP::detail::n.
Referenced by densniv().
00975 { 00976 G4double ucr = 10.0; // Critical energy for damping. 00977 G4double dcr = 40.0; // Width of damping. 00978 G4double ponq = 0.0, dn = 0.0, n = 0.0, dz = 0.0; 00979 00980 if(((std::fabs(bet)-1.15) < 0) || ((std::fabs(bet)-1.15) == 0)) { 00981 goto qrot10; 00982 } 00983 00984 if((std::fabs(bet)-1.15) > 0) { 00985 goto qrot11; 00986 } 00987 00988 qrot10: 00989 n = a - z; 00990 dz = std::fabs(z - 82.0); 00991 if (n > 104) { 00992 dn = std::fabs(n-126.e0); 00993 } 00994 else { 00995 dn = std::fabs(n - 82.0); 00996 } 00997 00998 bet = 0.022 + 0.003*dn + 0.005*dz; 00999 01000 sig = 25.0*std::pow(bet,2) * sig; 01001 01002 qrot11: 01003 ponq = (u - ucr)/dcr; 01004 01005 if (ponq > 700.0) { 01006 ponq = 700.0; 01007 } 01008 if (sig < 1.0) { 01009 sig = 1.0; 01010 } 01011 (*qr) = 1.0/(1.0 + std::exp(ponq)) * (sig - 1.0) + 1.0; 01012 01013 if ((*qr) < 1.0) { 01014 (*qr) = 1.0; 01015 } 01016 01017 return; 01018 }
Definition at line 4993 of file G4Abla.cc.
04994 { 04995 // C Rotation d'un vecteur 04996 // DIMENSION Pin(3),Pout(3) 04997 // REAL*8 R(3,3) 04998 04999 for(G4int i = 1; i <= 3; i++) { // do i=1,3 05000 pout[i] = 0.0; 05001 for(G4int j = 1; j <= 3; j++) { //do j=1,3 05002 // pout[i] = pout[i] + R[i][j]*pin[j]; 05003 pout[i] = pout[i] + R[j][i]*pin[j]; 05004 } // enddo 05005 } //enddo 05006 }
Definition at line 5171 of file G4Abla.cc.
Referenced by haz().
05172 { 05173 time_t mytime; 05174 tm *mylocaltime; 05175 05176 time(&mytime); 05177 mylocaltime = localtime(&mytime); 05178 05179 if(x == 0) { 05180 return(mylocaltime->tm_hour*60*60 + mylocaltime->tm_min*60 + mylocaltime->tm_sec); 05181 } 05182 else { 05183 return(mytime - x); 05184 } 05185 }
| void G4Abla::setVerboseLevel | ( | G4int | level | ) |
Set verbosity level.
Definition at line 1097 of file G4Abla.cc.
References fissility(), and idnint().
Referenced by direct().
01098 { 01099 01100 // INPUT: A,Z,OPTXFIS MASS AND CHARGE OF A NUCLEUS, 01101 // OPTION FOR FISSILITY 01102 // OUTPUT: SPDEF 01103 01104 // ALPHA2 SADDLE POINT DEF. COHEN&SWIATECKI ANN.PHYS. 22 (1963) 406 01105 // RANGING FROM FISSILITY X=0.30 TO X=1.00 IN STEPS OF 0.02 01106 01107 G4int index = 0; 01108 G4double x = 0.0, v = 0.0, dx = 0.0; 01109 01110 const G4int alpha2Size = 37; 01111 // The value 0.0 at alpha2[0] added by PK. 01112 G4double alpha2[alpha2Size] = {0.0, 2.5464e0, 2.4944e0, 2.4410e0, 2.3915e0, 2.3482e0, 01113 2.3014e0, 2.2479e0, 2.1982e0, 2.1432e0, 2.0807e0, 2.0142e0, 1.9419e0, 01114 1.8714e0, 1.8010e0, 1.7272e0, 1.6473e0, 1.5601e0, 1.4526e0, 1.3164e0, 01115 1.1391e0, 0.9662e0, 0.8295e0, 0.7231e0, 0.6360e0, 0.5615e0, 0.4953e0, 01116 0.4354e0, 0.3799e0, 0.3274e0, 0.2779e0, 0.2298e0, 0.1827e0, 0.1373e0, 01117 0.0901e0, 0.0430e0, 0.0000e0}; 01118 01119 dx = 0.02; 01120 x = fissility(a,z,optxfis); 01121 01122 if (x > 1.0) { 01123 x = 1.0; 01124 } 01125 01126 if (x < 0.0) { 01127 x = 0.0; 01128 } 01129 01130 v = (x - 0.3)/dx + 1.0; 01131 index = idnint(v); 01132 01133 if (index < 1) { 01134 return(alpha2[1]); // alpha2[0] -> alpha2[1] 01135 } 01136 01137 if (index == 36) { //then // :::FIXME:: Possible off-by-one bug... 01138 return(alpha2[36]); 01139 } 01140 else { 01141 return(alpha2[index] + (alpha2[index+1] - alpha2[index]) / dx * ( x - (0.3e0 + dx*(index-1)))); //:::FIXME::: Possible off-by-one 01142 } 01143 01144 return alpha2[0]; // The algorithm is not supposed to reach this point. 01145 }
Definition at line 5012 of file G4Abla.cc.
References G4UniformRand.
Referenced by direct(), evapora(), fmaxhaz(), and haz().
05013 { 05014 (*seed) = (*seed); // Avoid warning during compilation. 05015 // Use Geant4 G4UniformRand 05016 (*rndm) = G4UniformRand(); 05017 }
RISE TIME IN WHICH THE FISSION WIDTH HAS REACHED 90 PERCENT OF ITS FINAL VALUE
Definition at line 2764 of file G4Abla.cc.
Referenced by direct().
02765 { 02766 // INPUT : BET, HOMEGA, EF, T 02767 // OUTPUT: TAU - RISE TIME IN WHICH THE FISSION WIDTH HAS REACHED 02768 // 90 PERCENT OF ITS FINAL VALUE 02769 // 02770 // BETA - NUCLEAR VISCOSITY 02771 // HOMEGA - CURVATURE OF POTENTIAL 02772 // EF - FISSION BARRIER 02773 // T - NUCLEAR TEMPERATURE 02774 02775 G4double tauResult = 0.0; 02776 02777 G4double tlim = 8.e0 * ef; 02778 if (t > tlim) { 02779 t = tlim; 02780 } 02781 02782 // modified bj and khs 6.1.2000: 02783 if (bet/(std::sqrt(2.0)*10.0*(homega/6.582122)) <= 1.0) { 02784 tauResult = std::log(10.0*ef/t)/(bet*1.0e21); 02785 // assert(isnan(tauResult) == false); 02786 } 02787 else { 02788 tauResult = std::log(10.0*ef/t)/ (2.0*std::pow((10.0*homega/6.582122),2))*(bet*1.0e-21); 02789 // assert(isnan(tauResult) == false); 02790 } //end if 02791 02792 return tauResult; 02793 }
| void G4Abla::translab | ( | G4double | gamrem, | |
| G4double | etrem, | |||
| G4double | csrem[4], | |||
| G4int | nopart, | |||
| G4int | ndec | |||
| ) |
Definition at line 4752 of file G4Abla.cc.
References G4Volant::acv, G4cout, G4endl, G4Volant::iv, nint(), and G4Volant::zpcv.
Referenced by breakItUp().
04753 { 04754 // c Ce subroutine transforme dans un repere 1 les impulsions pcv des 04755 // c particules acv, zcv et de cosinus directeurs xcv, ycv, zcv calculees 04756 // c dans un repere 2. 04757 // c La transformation de lorentz est definie par GAMREM (gamma) et 04758 // c ETREM (eta). La direction du repere 2 dans 1 est donnees par les 04759 // c cosinus directeurs ALREM,BEREM,GAREM (axe oz du repere 2). 04760 // c L'axe oy(2) est fixe par le produit vectoriel oz(1)*oz(2). 04761 // c Le calcul est fait pour les particules de NDEC a iv du common volant. 04762 // C Resultats dans le NTUPLE (common VAR_NTP) decale de NOPART (cascade). 04763 04764 // REAL*8 GAMREM,ETREM,ER,PLABI(3),PLABF(3),R(3,3) 04765 // real*8 MASSE,PTRAV2,CSREM(3),UMA,MELEC,EL 04766 // real*4 acv,zpcv,pcv,xcv,ycv,zcv 04767 // common/volant/acv(200),zpcv(200),pcv(200),xcv(200), 04768 // s ycv(200),zcv(200),iv 04769 04770 // parameter (max=250) 04771 // real*4 EXINI,ENERJ,BIMPACT,PLAB,TETLAB,PHILAB,ESTFIS 04772 // integer AVV,ZVV,JREMN,KFIS,IZFIS,IAFIS 04773 // common/VAR_NTP/MASSINI,MZINI,EXINI,MULNCASC,MULNEVAP, 04774 // +MULNTOT,BIMPACT,JREMN,KFIS,ESTFIS,IZFIS,IAFIS,NTRACK, 04775 // +ITYPCASC(max),AVV(max),ZVV(max),ENERJ(max),PLAB(max), 04776 // +TETLAB(max),PHILAB(max) 04777 04778 // DATA UMA,MELEC/931.4942,0.511/ 04779 04780 // C Matrice de rotation dans le labo: 04781 G4double sitet = std::sqrt(std::pow(csrem[1],2)+std::pow(csrem[2],2)); 04782 G4double cstet = 0.0, siphi = 0.0, csphi = 0.0; 04783 G4double R[4][4]; 04784 for(G4int init_i = 0; init_i < 4; init_i++) { 04785 for(G4int init_j = 0; init_j < 4; init_j++) { 04786 R[init_i][init_j] = 0.0; 04787 } 04788 } 04789 04790 if(sitet > 1.0e-6) { //then 04791 cstet = csrem[3]; 04792 siphi = csrem[2]/sitet; 04793 csphi = csrem[1]/sitet; 04794 04795 R[1][1] = cstet*csphi; 04796 R[1][2] = -siphi; 04797 R[1][3] = sitet*csphi; 04798 R[2][1] = cstet*siphi; 04799 R[2][2] = csphi; 04800 R[2][3] = sitet*siphi; 04801 R[3][1] = -sitet; 04802 R[3][2] = 0.0; 04803 R[3][3] = cstet; 04804 } 04805 else { 04806 R[1][1] = 1.0; 04807 R[1][2] = 0.0; 04808 R[1][3] = 0.0; 04809 R[2][1] = 0.0; 04810 R[2][2] = 1.0; 04811 R[2][3] = 0.0; 04812 R[3][1] = 0.0; 04813 R[3][2] = 0.0; 04814 R[3][3] = 1.0; 04815 } //endif 04816 04817 G4int intp = 0; 04818 G4double el = 0.0; 04819 G4double masse = 0.0; 04820 G4double er = 0.0; 04821 G4double plabi[4]; 04822 G4double ptrav2 = 0.0; 04823 G4double plabf[4]; 04824 G4double bidon = 0.0; 04825 for(G4int init_i = 0; init_i < 4; init_i++) { 04826 plabi[init_i] = 0.0; 04827 plabf[init_i] = 0.0; 04828 } 04829 04830 for(G4int i = ndec; i <= volant->iv; i++) { //do i=ndec,iv 04831 intp = i + nopart; 04832 varntp->ntrack = varntp->ntrack + 1; 04833 if(nint(volant->acv[i]) == 0 && nint(volant->zpcv[i]) == 0) { 04834 if(verboseLevel > 2) { 04835 G4cout <<"Error: Particles with A = 0 Z = 0 detected! " << G4endl; 04836 } 04837 continue; 04838 } 04839 if(varntp->ntrack >= VARNTPSIZE) { 04840 if(verboseLevel > 2) { 04841 G4cout <<"Error! Output data structure not big enough!" << G4endl; 04842 } 04843 } 04844 varntp->avv[intp] = nint(volant->acv[i]); 04845 varntp->zvv[intp] = nint(volant->zpcv[i]); 04846 varntp->itypcasc[intp] = 0; 04847 // transformation de lorentz remnan --> labo: 04848 if (varntp->avv[intp] == -1) { //then 04849 masse = 138.00; // cugnon 04850 // c if (avv(intp).eq.1) masse=938.2796 !cugnon 04851 // c if (avv(intp).eq.4) masse=3727.42 !ok 04852 } 04853 else { 04854 mglms(double(volant->acv[i]),double(volant->zpcv[i]),0, &el); 04855 // assert(isnan(el) == false); 04856 masse = volant->zpcv[i]*938.27 + (volant->acv[i] - volant->zpcv[i])*939.56 + el; 04857 } //end if 04858 04859 er = std::sqrt(std::pow(volant->pcv[i],2) + std::pow(masse,2)); 04860 // assert(isnan(er) == false); 04861 plabi[1] = volant->pcv[i]*(volant->xcv[i]); 04862 plabi[2] = volant->pcv[i]*(volant->ycv[i]); 04863 plabi[3] = er*etrem + gamrem*(volant->pcv[i])*(volant->zcv[i]); 04864 04865 ptrav2 = std::pow(plabi[1],2) + std::pow(plabi[2],2) + std::pow(plabi[3],2); 04866 // assert(isnan(ptrav2) == false); 04867 varntp->plab[intp] = std::sqrt(ptrav2); 04868 varntp->enerj[intp] = std::sqrt(ptrav2 + std::pow(masse,2)) - masse; 04869 04870 // Rotation dans le labo: 04871 for(G4int j = 1; j <= 3; j++) { //do j=1,3 04872 plabf[j] = 0.0; 04873 for(G4int k = 1; k <= 3; k++) { //do k=1,3 04874 plabf[j] = plabf[j] + R[k][j]*plabi[k]; // :::Fixme::: (indices?) 04875 } // end do 04876 } // end do 04877 // C impulsions dans le nouveau systeme copiees dans /volant/ 04878 volant->pcv[i] = varntp->plab[intp]; 04879 ptrav2 = std::sqrt(std::pow(plabf[1],2) + std::pow(plabf[2],2) + std::pow(plabf[3],2)); 04880 if(ptrav2 >= 1.0e-6) { //then 04881 volant->xcv[i] = plabf[1]/ptrav2; 04882 volant->ycv[i] = plabf[2]/ptrav2; 04883 volant->zcv[i] = plabf[3]/ptrav2; 04884 } 04885 else { 04886 volant->xcv[i] = 1.0; 04887 volant->ycv[i] = 0.0; 04888 volant->zcv[i] = 0.0; 04889 } //endif 04890 // impulsions dans le nouveau systeme copiees dans /VAR_NTP/ 04891 if(varntp->plab[intp] >= 1.0e-6) { //then 04892 bidon = plabf[3]/(varntp->plab[intp]); 04893 // assert(isnan(bidon) == false); 04894 if(bidon > 1.0) { 04895 bidon = 1.0; 04896 } 04897 if(bidon < -1.0) { 04898 bidon = -1.0; 04899 } 04900 varntp->tetlab[intp] = std::acos(bidon); 04901 sitet = std::sin(varntp->tetlab[intp]); 04902 varntp->philab[intp] = std::atan2(plabf[2],plabf[1]); 04903 varntp->tetlab[intp] = varntp->tetlab[intp]*57.2957795; 04904 varntp->philab[intp] = varntp->philab[intp]*57.2957795; 04905 } 04906 else { 04907 varntp->tetlab[intp] = 90.0; 04908 varntp->philab[intp] = 0.0; 04909 } // endif 04910 } // end do 04911 }
| void G4Abla::translabpf | ( | G4double | masse1, | |
| G4double | t1, | |||
| G4double | p1, | |||
| G4double | ctet1, | |||
| G4double | phi1, | |||
| G4double | gamrem, | |||
| G4double | etrem, | |||
| G4double | R[][4], | |||
| G4double * | plab1, | |||
| G4double * | gam1, | |||
| G4double * | eta1, | |||
| G4double | csdir[] | |||
| ) |
Definition at line 4916 of file G4Abla.cc.
04919 { 04920 // C Calcul de l'impulsion du PF (PLAB1, cos directeurs CSDIR(3)) dans le 04921 // C systeme remnant et des coefs de Lorentz GAM1,ETA1 de passage 04922 // c du systeme PF --> systeme remnant. 04923 // c 04924 // C Input: MASSE1, T1 (energie cinetique), CTET1,PHI1 (cosTHETA et PHI) 04925 // C (le PF dans le systeme du Noyau de Fission (NF)). 04926 // C GAMREM,ETREM les coefs de Lorentz systeme NF --> syst remnant, 04927 // C R(3,3) la matrice de rotation systeme NF--> systeme remnant. 04928 // C 04929 // C 04930 // REAL*8 MASSE1,T1,P1,CTET1,PHI1,GAMREM,ETREM,R(3,3), 04931 // s PLAB1,GAM1,ETA1,CSDIR(3),ER,SITET,PLABI(3),PLABF(3) 04932 04933 G4double er = t1 + masse1; 04934 04935 G4double sitet = std::sqrt(1.0 - std::pow(ctet1,2)); 04936 04937 G4double plabi[4]; 04938 G4double plabf[4]; 04939 for(G4int init_i = 0; init_i < 4; init_i++) { 04940 plabi[init_i] = 0.0; 04941 plabf[init_i] = 0.0; 04942 } 04943 04944 // C ----Transformation de Lorentz Noyau fissionnant --> Remnant: 04945 plabi[1] = p1*sitet*std::cos(phi1); 04946 plabi[2] = p1*sitet*std::sin(phi1); 04947 plabi[3] = er*etrem + gamrem*p1*ctet1; 04948 04949 // C ----Rotation du syst Noyaut Fissionant vers syst remnant: 04950 for(G4int j = 1; j <= 3; j++) { // do j=1,3 04951 plabf[j] = 0.0; 04952 for(G4int k = 1; k <= 3; k++) { //do k=1,3 04953 // plabf[j] = plabf[j] + R[j][k]*plabi[k]; 04954 plabf[j] = plabf[j] + R[k][j]*plabi[k]; 04955 } //end do 04956 } //end do 04957 // C ----Cosinus directeurs et coefs de la transf de Lorentz dans le 04958 // c nouveau systeme: 04959 (*plab1) = std::pow(plabf[1],2) + std::pow(plabf[2],2) + std::pow(plabf[3],2); 04960 (*gam1) = std::sqrt(std::pow(masse1,2) + (*plab1))/masse1; 04961 (*plab1) = std::sqrt((*plab1)); 04962 (*eta1) = (*plab1)/masse1; 04963 04964 if((*plab1) <= 1.0e-6) { //then 04965 csdir[1] = 0.0; 04966 csdir[2] = 0.0; 04967 csdir[3] = 1.0; 04968 } 04969 else { 04970 for(G4int i = 1; i <= 3; i++) { //do i=1,3 04971 csdir[i] = plabf[i]/(*plab1); 04972 } // end do 04973 } //endif 04974 }
Definition at line 3452 of file G4Abla.cc.
References G4InuclParticleNames::alpha, and G4INCL::Math::pi.
Referenced by fissionDistri().
03453 { 03454 // liquid-drop mass, Myers & Swiatecki, Lysekil, 1967 03455 // pure liquid drop, without pairing and shell effects 03456 03457 // On input: Z nuclear charge of nucleus 03458 // N number of neutrons in nucleus 03459 // beta deformation of nucleus 03460 // On output: binding energy of nucleus 03461 03462 G4double a = 0.0, umass = 0.0; 03463 G4double alpha = 0.0; 03464 G4double xcom = 0.0, xvs = 0.0, xe = 0.0; 03465 const G4double pi = 3.1416; 03466 03467 a = n + z; 03468 alpha = ( std::sqrt(5.0/(4.0*pi)) ) * beta; 03469 // assert(isnan(alpha) == false); 03470 03471 xcom = 1.0 - 1.7826 * ((a - 2.0*z)/a)*((a - 2.0*z)/a); 03472 // assert(isnan(xcom) == false); 03473 // factor for asymmetry dependence of surface and volume term 03474 xvs = - xcom * ( 15.4941 * a - 03475 17.9439 * std::pow(a,0.66667) * (1.0+0.4*alpha*alpha) ); 03476 // sum of volume and surface energy 03477 xe = z*z * (0.7053/(std::pow(a,0.33333)) * (1.0-0.2*alpha*alpha) - 1.1529/a); 03478 // assert(isnan(xe) == false); 03479 umass = xvs + xe; 03480 03481 return umass; 03482 }
1.4.7