G4Fissioner.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
00029 // 20100318  M. Kelsey -- Bug fix setting mass with G4LV
00030 // 20100319  M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
00031 //              eliminate some unnecessary std::pow()
00032 // 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
00033 // 20100517  M. Kelsey -- Inherit from common base class
00034 // 20100622  M. Kelsey -- Use local "bindingEnergy()" to call through
00035 // 20100711  M. Kelsey -- Add energy-conservation checking, reduce if-cascades
00036 // 20100713  M. Kelsey -- Don't add excitation energy to mass (already there)
00037 // 20100714  M. Kelsey -- Move conservation checking to base class
00038 // 20100728  M. Kelsey -- Make fixed arrays static, move G4FissionStore to data
00039 //              member and reuse
00040 // 20100914  M. Kelsey -- Migrate to integer A and Z
00041 // 20110214  M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
00042 // 20110801  M. Kelsey -- Replace constant-size std::vector's w/C arrays
00043 // 20110922  M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
00044 // 20120517  A. Ribon  -- Removed static in local vectors for reproducibility
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* /*bullet*/,
00061                           G4InuclParticle* target,
00062                           G4CollisionOutput& output) {
00063   if (verboseLevel) {
00064     G4cout << " >>> G4Fissioner::collide" << G4endl;
00065   }
00066 
00067   //  const G4int itry_max = 1000;
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   // Initialize buffer for fission possibilities
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;              // Mass includes excitation
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) { // generate fluctuations
00134       //  faster, using randomGauss
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;          // No fission products
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;                // No fission energy
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   // Pass only last two nuclear fragments
00184   static std::vector<G4InuclNuclei> frags(2);           // Always the same size!
00185   frags[0] = nuclei1;
00186   frags[1] = nuclei2;
00187   validateOutput(0, target, frags);             // Check energy conservation
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 }

Generated on Mon May 27 17:48:16 2013 for Geant4 by  doxygen 1.4.7