00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 #include "G4Fissioner.hh"
00047 #include "G4CollisionOutput.hh"
00048 #include "G4HadTmpUtil.hh"
00049 #include "G4InuclNuclei.hh"
00050 #include "G4InuclParticle.hh"
00051 #include "G4FissionStore.hh"
00052 #include "G4FissionConfiguration.hh"
00053 #include "G4InuclSpecialFunctions.hh"
00054
00055 using namespace G4InuclSpecialFunctions;
00056
00057
00058 G4Fissioner::G4Fissioner() : G4CascadeColliderBase("G4Fissioner") {}
00059
00060 void G4Fissioner::collide(G4InuclParticle* ,
00061 G4InuclParticle* target,
00062 G4CollisionOutput& output) {
00063 if (verboseLevel) {
00064 G4cout << " >>> G4Fissioner::collide" << G4endl;
00065 }
00066
00067
00068
00069 G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target);
00070 if (!nuclei_target) {
00071 G4cerr << " >>> G4Fissioner -> target is not nuclei " << G4endl;
00072 return;
00073 }
00074
00075 if (verboseLevel > 1)
00076 G4cout << " Fissioner input\n" << *nuclei_target << G4endl;
00077
00078
00079 fissionStore.setVerboseLevel(verboseLevel);
00080 fissionStore.clear();
00081
00082 G4int A = nuclei_target->getA();
00083 G4int Z = nuclei_target->getZ();
00084
00085 G4double EEXS = nuclei_target->getExitationEnergy();
00086 G4double mass_in = nuclei_target->getMass();
00087 G4double e_in = mass_in;
00088 G4double PARA = 0.055 * G4cbrt(A*A) * (G4cbrt(A-Z) + G4cbrt(Z));
00089 G4double TEM = std::sqrt(EEXS / PARA);
00090 G4double TETA = 0.494 * G4cbrt(A) * TEM;
00091
00092 TETA = TETA / std::sinh(TETA);
00093
00094 if (A < 246) PARA += (nucleiLevelDensity(A) - PARA) * TETA;
00095
00096 G4int A1 = A/2 + 1;
00097 G4int Z1;
00098 G4int A2 = A - A1;
00099
00100 G4double ALMA = -1000.0;
00101 G4double DM1 = bindingEnergy(A,Z);
00102 G4double EVV = EEXS - DM1;
00103 G4double DM2 = bindingEnergyAsymptotic(A, Z);
00104 G4double DTEM = (A < 220 ? 0.5 : 1.15);
00105
00106 TEM += DTEM;
00107
00108 G4double AL1[2] = { -0.15, -0.15 };
00109 G4double BET1[2] = { 0.05, 0.05 };
00110
00111 G4double R12 = G4cbrt(A1) + G4cbrt(A2);
00112
00113 for (G4int i = 0; i < 50 && A1 > 30; i++) {
00114 A1--;
00115 A2 = A - A1;
00116 G4double X3 = 1.0 / G4cbrt(A1);
00117 G4double X4 = 1.0 / G4cbrt(A2);
00118 Z1 = G4lrint(getZopt(A1, A2, Z, X3, X4, R12) - 1.);
00119 G4double EDEF1[2];
00120 G4int Z2 = Z - Z1;
00121 G4double VPOT, VCOUL;
00122
00123 potentialMinimization(VPOT, EDEF1, VCOUL, A1, A2, Z1, Z2, AL1, BET1, R12);
00124
00125 G4double DM3 = bindingEnergy(A1,Z1);
00126 G4double DM4 = bindingEnergyAsymptotic(A1, Z1);
00127 G4double DM5 = bindingEnergy(A2,Z2);
00128 G4double DM6 = bindingEnergyAsymptotic(A2, Z2);
00129 G4double DMT1 = DM4 + DM6 - DM2;
00130 G4double DMT = DM3 + DM5 - DM1;
00131 G4double EZL = EEXS + DMT - VPOT;
00132
00133 if(EZL > 0.0) {
00134
00135 G4double C1 = std::sqrt(getC2(A1, A2, X3, X4, R12) / TEM);
00136 G4double DZ = randomGauss(C1);
00137
00138 DZ = DZ > 0.0 ? DZ + 0.5 : -std::fabs(DZ - 0.5);
00139 Z1 += G4int(DZ);
00140 Z2 -= G4int(DZ);
00141
00142 G4double DEfin = randomGauss(TEM);
00143 G4double EZ = (DMT1 + (DMT - DMT1) * TETA - VPOT + DEfin) / TEM;
00144
00145 if (EZ >= ALMA) ALMA = EZ;
00146 G4double EK = VCOUL + DEfin + 0.5 * TEM;
00147 G4double EV = EVV + bindingEnergy(A1,Z1) + bindingEnergy(A2,Z2) - EK;
00148
00149 if (EV > 0.0) fissionStore.addConfig(A1, Z1, EZ, EK, EV);
00150 };
00151 };
00152
00153 G4int store_size = fissionStore.size();
00154 if (store_size == 0) return;
00155
00156 G4FissionConfiguration config =
00157 fissionStore.generateConfiguration(ALMA, inuclRndm());
00158
00159 A1 = G4int(config.afirst);
00160 A2 = A - A1;
00161 Z1 = G4int(config.zfirst);
00162
00163 G4int Z2 = Z - Z1;
00164
00165 G4double mass1 = G4InuclNuclei::getNucleiMass(A1,Z1);
00166 G4double mass2 = G4InuclNuclei::getNucleiMass(A2,Z2);
00167 G4double EK = config.ekin;
00168 G4double pmod = std::sqrt(0.001 * EK * mass1 * mass2 / mass_in);
00169
00170 G4LorentzVector mom1 = generateWithRandomAngles(pmod, mass1);
00171 G4LorentzVector mom2; mom2.setVectM(-mom1.vect(), mass2);
00172
00173 G4double e_out = mom1.e() + mom2.e();
00174 G4double EV = 1000.0 * (e_in - e_out) / A;
00175 if (EV <= 0.0) return;
00176
00177 G4double EEXS1 = EV*A1;
00178 G4double EEXS2 = EV*A2;
00179
00180 G4InuclNuclei nuclei1(mom1, A1, Z1, EEXS1, G4InuclParticle::Fissioner);
00181 G4InuclNuclei nuclei2(mom2, A2, Z2, EEXS2, G4InuclParticle::Fissioner);
00182
00183
00184 static std::vector<G4InuclNuclei> frags(2);
00185 frags[0] = nuclei1;
00186 frags[1] = nuclei2;
00187 validateOutput(0, target, frags);
00188
00189 output.addOutgoingNuclei(frags);
00190 }
00191
00192 G4double G4Fissioner::getC2(G4int A1,
00193 G4int A2,
00194 G4double X3,
00195 G4double X4,
00196 G4double R12) const {
00197
00198 if (verboseLevel > 3) {
00199 G4cout << " >>> G4Fissioner::getC2" << G4endl;
00200 }
00201
00202 G4double C2 = 124.57 * (1.0 / A1 + 1.0 / A2) + 0.78 * (X3 + X4) - 176.9 *
00203 ((X3*X3*X3*X3) + (X4*X4*X4*X4)) + 219.36 * (1.0 / (A1 * A1) + 1.0 / (A2 * A2)) - 1.108 / R12;
00204
00205 return C2;
00206 }
00207
00208 G4double G4Fissioner::getZopt(G4int A1,
00209 G4int A2,
00210 G4int ZT,
00211 G4double X3,
00212 G4double X4,
00213 G4double R12) const {
00214
00215 if (verboseLevel > 3) {
00216 G4cout << " >>> G4Fissioner::getZopt" << G4endl;
00217 }
00218
00219 G4double Zopt = (87.7 * (X4 - X3) * (1.0 - 1.25 * (X4 + X3)) +
00220 ZT * ((124.57 / A2 + 0.78 * X4 - 176.9 * (X4*X4*X4*X4) + 219.36 / (A2 * A2)) - 0.554 / R12)) /
00221 getC2(A1, A2, X3, X4, R12);
00222
00223 return Zopt;
00224 }
00225
00226 void G4Fissioner::potentialMinimization(G4double& VP,
00227 G4double( &ED)[2],
00228 G4double& VC,
00229 G4int AF,
00230 G4int AS,
00231 G4int ZF,
00232 G4int ZS,
00233 G4double AL1[2],
00234 G4double BET1[2],
00235 G4double& R12) const {
00236
00237 if (verboseLevel > 3) {
00238 G4cout << " >>> G4Fissioner::potentialMinimization" << G4endl;
00239 }
00240
00241 const G4double huge_num = 2.0e35;
00242 const G4int itry_max = 2000;
00243 const G4double DSOL1 = 1.0e-6;
00244 const G4double DS1 = 0.3;
00245 const G4double DS2 = 1.0 / DS1 / DS1;
00246 G4int A1[2] = { AF, AS };
00247 G4int Z1[2] = { ZF, ZS };
00248 G4double D = 1.01844 * ZF * ZS;
00249 G4double D0 = 1.0e-3 * D;
00250 G4double R[2];
00251 R12 = 0.0;
00252 G4double C[2];
00253 G4double F[2];
00254 G4double Y1;
00255 G4double Y2;
00256 G4int i;
00257
00258 for (i = 0; i < 2; i++) {
00259 R[i] = G4cbrt(A1[i]);
00260 Y1 = R[i] * R[i];
00261 Y2 = Z1[i] * Z1[i] / R[i];
00262 C[i] = 6.8 * Y1 - 0.142 * Y2;
00263 F[i] = 12.138 * Y1 - 0.145 * Y2;
00264 };
00265
00266 G4double SAL[2];
00267 G4double SBE[2];
00268 G4double X[2];
00269 G4double X1[2];
00270 G4double X2[2];
00271 G4double RAL[2];
00272 G4double RBE[2];
00273 G4double A[4][4];
00274 G4double B[4];
00275 G4int itry = 0;
00276
00277 while (itry < itry_max) {
00278 itry++;
00279 G4double S = 0.0;
00280
00281 for (i = 0; i < 2; i++) {
00282 S += R[i] * (1.0 + AL1[i] + BET1[i] - 0.257 * AL1[i] * BET1[i]);
00283 };
00284 R12 = 0.0;
00285 Y1 = 0.0;
00286 Y2 = 0.0;
00287
00288 for (i = 0; i < 2; i++) {
00289 SAL[i] = R[i] * (1.0-0.257 * BET1[i]);
00290 SBE[i] = R[i] * (1.0-0.257 * AL1[i]);
00291 X[i] = R[i] / S;
00292 X1[i] = X[i] * X[i];
00293 X2[i] = X[i] * X1[i];
00294 Y1 += AL1[i] * X1[i];
00295 Y2 += BET1[i] * X2[i];
00296 R12 += R[i] * (1.0 - AL1[i] * (1.0 - 0.6 * X[i]) + BET1[i] * (1.0 - 0.429 * X1[i]));
00297 };
00298
00299 G4double Y3 = -0.6 * Y1 + 0.857 * Y2;
00300 G4double Y4 = (1.2 * Y1 - 2.571 * Y2) / S;
00301 G4double R2 = D0 / (R12 * R12);
00302 G4double R3 = 2.0 * R2 / R12;
00303
00304 for (i = 0; i < 2; i++) {
00305 RAL[i] = -R[i] * (1.0 - 0.6 * X[i]) + SAL[i] * Y3;
00306 RBE[i] = R[i] * (1.0 - 0.429 * X1[i]) + SBE[i] * Y3;
00307 };
00308
00309 G4double DX1;
00310 G4double DX2;
00311
00312 for (i = 0; i < 2; i++) {
00313
00314 for (G4int j = 0; j < 2; j++) {
00315 G4double DEL1 = i == j ? 1.0 : 0.0;
00316 DX1 = 0.0;
00317 DX2 = 0.0;
00318
00319 if (std::fabs(AL1[i]) >= DS1) {
00320 G4double XXX = AL1[i] * AL1[i] * DS2;
00321 G4double DEX = XXX > 100.0 ? huge_num : std::exp(XXX);
00322 DX1 = 2.0 * (1.0 + 2.0 * AL1[i] * AL1[i] * DS2) * DEX * DS2;
00323 };
00324
00325 if (std::fabs(BET1[i]) >= DS1) {
00326 G4double XXX = BET1[i] * BET1[i] * DS2;
00327 G4double DEX = XXX > 100.0 ? huge_num : std::exp(XXX);
00328 DX2 = 2.0 * (1.+2.0 * BET1[i] * BET1[i] * DS2) * DEX * DS2;
00329 };
00330
00331 G4double DEL = 2.0e-3 * DEL1;
00332 A[i][j] = R3 * RBE[i] * RBE[j] -
00333 R2 * (-0.6 *
00334 (X1[i] * SAL[j] +
00335 X1[j] * SAL[i]) + SAL[i] * SAL[j] * Y4) +
00336 DEL * C[i] + DEL1 * DX1;
00337 G4int i1 = i + 2;
00338 G4int j1 = j + 2;
00339 A[i1][j1] = R3 * RBE[i] * RBE[j]
00340 - R2 * (0.857 *
00341 (X2[i] * SBE[j] +
00342 X2[j] * SBE[i]) + SBE[i] * SBE[j] * Y4) +
00343 DEL * F[i] + DEL1 * DX2;
00344 A[i][j1] = R3 * RAL[i] * RBE[j] -
00345 R2 * (0.857 *
00346 (X2[j] * SAL[i] -
00347 0.6 * X1[i] * SBE[j]) + SBE[j] * SAL[i] * Y4 -
00348 0.257 * R[i] * Y3 * DEL1);
00349 A[j1][i] = A[i][j1];
00350 };
00351 };
00352
00353 for (i = 0; i < 2; i++) {
00354 DX1 = 0.0;
00355 DX2 = 0.0;
00356
00357 if (std::fabs(AL1[i]) >= DS1) DX1 = 2.0 * AL1[i] * DS2 * std::exp(AL1[i] * AL1[i] * DS2);
00358
00359 if (std::fabs(BET1[i]) >= DS1) DX2 = 2.0 * BET1[i] * DS2 * std::exp(BET1[i] * BET1[i] * DS2);
00360 B[i] = R2 * RAL[i] - 2.0e-3 * C[i] * AL1[i] + DX1;
00361 B[i + 2] = R2 * RBE[i] - 2.0e-3 * F[i] * BET1[i] + DX2;
00362 };
00363
00364 G4double ST = 0.0;
00365 G4double ST1 = 0.0;
00366
00367 for (i = 0; i < 4; i++) {
00368 ST += B[i] * B[i];
00369
00370 for (G4int j = 0; j < 4; j++) ST1 += A[i][j] * B[i] * B[j];
00371 };
00372
00373 G4double STEP = ST / ST1;
00374 G4double DSOL = 0.0;
00375
00376 for (i = 0; i < 2; i++) {
00377 AL1[i] += B[i] * STEP;
00378 BET1[i] += B[i + 2] * STEP;
00379 DSOL += B[i] * B[i] + B[i + 2] * B[i + 2];
00380 };
00381 DSOL = std::sqrt(DSOL);
00382
00383 if (DSOL < DSOL1) break;
00384 };
00385
00386 if (verboseLevel > 3) {
00387 if (itry == itry_max)
00388 G4cout << " maximal number of iterations in potentialMinimization " << G4endl
00389 << " A1 " << AF << " Z1 " << ZF << G4endl;
00390
00391 };
00392
00393 for (i = 0; i < 2; i++) ED[i] = F[i] * BET1[i] * BET1[i] + C[i] * AL1[i] * AL1[i];
00394
00395 VC = D / R12;
00396 VP = VC + ED[0] + ED[1];
00397 }