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 // ------------------------------------------------------------------- 00029 // 00030 // GEANT4 Class file 00031 // 00032 // 00033 // File name: G4PreCompoundTransitions 00034 // 00035 // Author: V.Lara 00036 // 00037 // Modified: 00038 // 16.02.2008 J.M.Quesada fixed bugs 00039 // 06.09.2008 J.M.Quesada added external choices for: 00040 // - "never go back" hipothesis (useNGB=true) 00041 // - CEM transition probabilities (useCEMtr=true) 00042 // 30.10.2009 J.M.Quesada: CEM transition probabilities have been renormalized 00043 // (IAEA benchmark) 00044 // 20.08.2010 V.Ivanchenko move constructor and destructor to the source and 00045 // optimise the code 00046 // 30.08.2011 M.Kelsey - Skip CalculateProbability if no excitons 00047 00048 #include "G4PreCompoundTransitions.hh" 00049 #include "G4PhysicalConstants.hh" 00050 #include "G4SystemOfUnits.hh" 00051 #include "Randomize.hh" 00052 #include "G4Pow.hh" 00053 #include "G4HadronicException.hh" 00054 #include "G4PreCompoundParameters.hh" 00055 #include "G4Proton.hh" 00056 00057 G4PreCompoundTransitions::G4PreCompoundTransitions() 00058 { 00059 proton = G4Proton::Proton(); 00060 FermiEnergy = G4PreCompoundParameters::GetAddress()->GetFermiEnergy(); 00061 r0 = G4PreCompoundParameters::GetAddress()->GetTransitionsr0(); 00062 aLDP = G4PreCompoundParameters::GetAddress()->GetLevelDensity(); 00063 g4pow = G4Pow::GetInstance(); 00064 } 00065 00066 G4PreCompoundTransitions::~G4PreCompoundTransitions() 00067 {} 00068 00069 // Calculates transition probabilities with 00070 // DeltaN = +2 (Trans1) -2 (Trans2) and 0 (Trans3) 00071 G4double G4PreCompoundTransitions:: 00072 CalculateProbability(const G4Fragment & aFragment) 00073 { 00074 // Number of holes 00075 G4int H = aFragment.GetNumberOfHoles(); 00076 // Number of Particles 00077 G4int P = aFragment.GetNumberOfParticles(); 00078 // Number of Excitons 00079 G4int N = P+H; 00080 // Nucleus 00081 G4int A = aFragment.GetA_asInt(); 00082 G4int Z = aFragment.GetZ_asInt(); 00083 G4double U = aFragment.GetExcitationEnergy(); 00084 00085 //G4cout << aFragment << G4endl; 00086 00087 if(U < 10*eV || 0==N) { return 0.0; } 00088 00089 //J. M. Quesada (Feb. 08) new physics 00090 // OPT=1 Transitions are calculated according to Gudima's paper 00091 // (original in G4PreCompound from VL) 00092 // OPT=2 Transitions are calculated according to Gupta's formulae 00093 // 00094 if (useCEMtr){ 00095 00096 // Relative Energy (T_{rel}) 00097 G4double RelativeEnergy = 1.6*FermiEnergy + U/G4double(N); 00098 00099 // Sample kind of nucleon-projectile 00100 G4bool ChargedNucleon(false); 00101 G4double chtest = 0.5; 00102 if (P > 0) { 00103 chtest = G4double(aFragment.GetNumberOfCharged())/G4double(P); 00104 } 00105 if (G4UniformRand() < chtest) { ChargedNucleon = true; } 00106 00107 // Relative Velocity: 00108 // <V_{rel}>^2 00109 G4double RelativeVelocitySqr(0.0); 00110 if (ChargedNucleon) { 00111 RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::proton_mass_c2; 00112 } else { 00113 RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::neutron_mass_c2; 00114 } 00115 00116 // <V_{rel}> 00117 G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr); 00118 00119 // Proton-Proton Cross Section 00120 G4double ppXSection = 00121 (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9) 00122 * CLHEP::millibarn; 00123 // Proton-Neutron Cross Section 00124 G4double npXSection = 00125 (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2) 00126 * CLHEP::millibarn; 00127 00128 // Averaged Cross Section: \sigma(V_{rel}) 00129 // G4double AveragedXSection = (ppXSection+npXSection)/2.0; 00130 G4double AveragedXSection(0.0); 00131 if (ChargedNucleon) 00132 { 00133 //JMQ: small bug fixed 00134 //AveragedXSection=((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection)/(A-1.0); 00135 AveragedXSection = ((Z-1)*ppXSection + (A-Z)*npXSection)/G4double(A-1); 00136 } 00137 else 00138 { 00139 AveragedXSection = ((A-Z-1)*ppXSection + Z*npXSection)/G4double(A-1); 00140 //AveragedXSection = ((A-Z-1)*npXSection + Z*ppXSection)/G4double(A-1); 00141 } 00142 00143 // Fermi relative energy ratio 00144 G4double FermiRelRatio = FermiEnergy/RelativeEnergy; 00145 00146 // This factor is introduced to take into account the Pauli principle 00147 G4double PauliFactor = 1.0 - 1.4*FermiRelRatio; 00148 if (FermiRelRatio > 0.5) { 00149 G4double x = 2.0 - 1.0/FermiRelRatio; 00150 PauliFactor += 0.4*FermiRelRatio*x*x*std::sqrt(x); 00151 //PauliFactor += 00152 //(2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0); 00153 } 00154 // Interaction volume 00155 // G4double Vint = (4.0/3.0) 00156 //*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0); 00157 G4double xx = 2.0*r0 + hbarc/(CLHEP::proton_mass_c2*RelativeVelocity); 00158 // G4double Vint = (4.0/3.0)*CLHEP::pi*xx*xx*xx; 00159 G4double Vint = CLHEP::pi*xx*xx*xx/0.75; 00160 00161 // Transition probability for \Delta n = +2 00162 00163 TransitionProb1 = AveragedXSection*PauliFactor 00164 *std::sqrt(2.0*RelativeEnergy/CLHEP::proton_mass_c2)/Vint; 00165 00166 //JMQ 281009 phenomenological factor in order to increase 00167 // equilibrium contribution 00168 // G4double factor=5.0; 00169 // TransitionProb1 *= factor; 00170 // 00171 if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; } 00172 00173 // GE = g*E where E is Excitation Energy 00174 G4double GE = (6.0/pi2)*aLDP*A*U; 00175 00176 //G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0); 00177 G4double Fph = G4double(P*P+H*H+P-3*H)/4.0; 00178 00179 G4bool NeverGoBack(false); 00180 if(useNGB) { NeverGoBack=true; } 00181 00182 //JMQ/AH bug fixed: if (U-Fph < 0.0) NeverGoBack = true; 00183 if (GE-Fph < 0.0) { NeverGoBack = true; } 00184 00185 // F(p+1,h+1) 00186 G4double Fph1 = Fph + N/2.0; 00187 00188 G4double ProbFactor = g4pow->powN((GE-Fph)/(GE-Fph1),N+1); 00189 00190 if (NeverGoBack) 00191 { 00192 TransitionProb2 = 0.0; 00193 TransitionProb3 = 0.0; 00194 } 00195 else 00196 { 00197 // Transition probability for \Delta n = -2 (at F(p,h) = 0) 00198 TransitionProb2 = 00199 TransitionProb1 * ProbFactor * (P*H*(N+1)*(N-2))/((GE-Fph)*(GE-Fph)); 00200 if (TransitionProb2 < 0.0) { TransitionProb2 = 0.0; } 00201 00202 // Transition probability for \Delta n = 0 (at F(p,h) = 0) 00203 TransitionProb3 = TransitionProb1*(N+1)* ProbFactor 00204 * (P*(P-1) + 4.0*P*H + H*(H-1))/(N*(GE-Fph)); 00205 if (TransitionProb3 < 0.0) { TransitionProb3 = 0.0; } 00206 } 00207 } else { 00208 //JMQ: Transition probabilities from Gupta's work 00209 // GE = g*E where E is Excitation Energy 00210 G4double GE = (6.0/pi2)*aLDP*A*U; 00211 00212 G4double Kmfp=2.; 00213 00214 //TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U); 00215 TransitionProb1 = 3.0e-9*(1.4e+21*U - 1.2e+19*U*U/G4double(N+1)) 00216 /(8*Kmfp*CLHEP::c_light); 00217 if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; } 00218 00219 TransitionProb2=0.; 00220 TransitionProb3=0.; 00221 00222 if (!useNGB && N > 1) { 00223 // TransitionProb2=1./Kmfp*3./8.*1./c_light*1.0e-9*(N-1.)*(N-2.)*P*H/(GE*GE)*(1.4e+21*U - 2./(N-1)*6.0e+18*U*U); 00224 TransitionProb2 = 00225 3.0e-9*(N-2)*P*H*(1.4e+21*U*(N-1) - 1.2e+19*U*U)/(8*Kmfp*c_light*GE*GE); 00226 if (TransitionProb2 < 0.0) TransitionProb2 = 0.0; 00227 } 00228 } 00229 // G4cout<<"U = "<<U<<G4endl; 00230 // G4cout<<"N="<<N<<" P="<<P<<" H="<<H<<G4endl; 00231 // G4cout<<"l+ ="<<TransitionProb1<<" l- ="<< TransitionProb2 00232 // <<" l0 ="<< TransitionProb3<<G4endl; 00233 return TransitionProb1 + TransitionProb2 + TransitionProb3; 00234 } 00235 00236 void G4PreCompoundTransitions::PerformTransition(G4Fragment & result) 00237 { 00238 G4double ChosenTransition = 00239 G4UniformRand()*(TransitionProb1 + TransitionProb2 + TransitionProb3); 00240 G4int deltaN = 0; 00241 // G4int Nexcitons = result.GetNumberOfExcitons(); 00242 G4int Npart = result.GetNumberOfParticles(); 00243 G4int Ncharged = result.GetNumberOfCharged(); 00244 G4int Nholes = result.GetNumberOfHoles(); 00245 if (ChosenTransition <= TransitionProb1) 00246 { 00247 // Number of excitons is increased on \Delta n = +2 00248 deltaN = 2; 00249 } 00250 else if (ChosenTransition <= TransitionProb1+TransitionProb2) 00251 { 00252 // Number of excitons is increased on \Delta n = -2 00253 deltaN = -2; 00254 } 00255 00256 // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and 00257 // in proportion to the number charges w.r.t. number of particles, 00258 // PROVIDED that there are charged particles 00259 deltaN /= 2; 00260 00261 //G4cout << "deltaN= " << deltaN << G4endl; 00262 00263 // JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on 00264 // number of charged vs. number of particles fails 00265 result.SetNumberOfParticles(Npart+deltaN); 00266 result.SetNumberOfHoles(Nholes+deltaN); 00267 00268 if(deltaN < 0) { 00269 if( Ncharged >= 1 && G4int(Npart*G4UniformRand()) <= Ncharged) 00270 { 00271 result.SetNumberOfCharged(Ncharged+deltaN); // deltaN is negative! 00272 } 00273 00274 } else if ( deltaN > 0 ) { 00275 // With weight Z/A, number of charged particles is increased with +1 00276 G4int A = result.GetA_asInt(); 00277 G4int Z = result.GetZ_asInt(); 00278 if( G4int(std::max(1, A - Npart)*G4UniformRand()) <= Z) 00279 { 00280 result.SetNumberOfCharged(Ncharged+deltaN); 00281 } 00282 } 00283 00284 // Number of charged can not be greater that number of particles 00285 if ( Npart < Ncharged ) 00286 { 00287 result.SetNumberOfCharged(Npart); 00288 } 00289 //G4cout << "### After transition" << G4endl; 00290 //G4cout << result << G4endl; 00291 } 00292