#include <G4PreCompoundTransitions.hh>
Inheritance diagram for G4PreCompoundTransitions:
Public Member Functions | |
G4PreCompoundTransitions () | |
virtual | ~G4PreCompoundTransitions () |
virtual G4double | CalculateProbability (const G4Fragment &aFragment) |
virtual void | PerformTransition (G4Fragment &aFragment) |
Definition at line 52 of file G4PreCompoundTransitions.hh.
G4PreCompoundTransitions::G4PreCompoundTransitions | ( | ) |
Definition at line 57 of file G4PreCompoundTransitions.cc.
References G4PreCompoundParameters::GetAddress(), G4PreCompoundParameters::GetFermiEnergy(), G4Pow::GetInstance(), G4PreCompoundParameters::GetLevelDensity(), G4PreCompoundParameters::GetTransitionsr0(), and G4Proton::Proton().
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 }
G4PreCompoundTransitions::~G4PreCompoundTransitions | ( | ) | [virtual] |
G4double G4PreCompoundTransitions::CalculateProbability | ( | const G4Fragment & | aFragment | ) | [virtual] |
Implements G4VPreCompoundTransitions.
Definition at line 72 of file G4PreCompoundTransitions.cc.
References G4UniformRand, GE, G4Fragment::GetA_asInt(), G4Fragment::GetExcitationEnergy(), G4Fragment::GetNumberOfCharged(), G4Fragment::GetNumberOfHoles(), G4Fragment::GetNumberOfParticles(), G4Fragment::GetZ_asInt(), G4INCL::Math::pi, G4Pow::powN(), G4VPreCompoundTransitions::TransitionProb1, G4VPreCompoundTransitions::TransitionProb2, G4VPreCompoundTransitions::TransitionProb3, G4VPreCompoundTransitions::useCEMtr, and G4VPreCompoundTransitions::useNGB.
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 }
void G4PreCompoundTransitions::PerformTransition | ( | G4Fragment & | aFragment | ) | [virtual] |
Implements G4VPreCompoundTransitions.
Definition at line 236 of file G4PreCompoundTransitions.cc.
References G4UniformRand, G4Fragment::GetA_asInt(), G4Fragment::GetNumberOfCharged(), G4Fragment::GetNumberOfHoles(), G4Fragment::GetNumberOfParticles(), G4Fragment::GetZ_asInt(), G4Fragment::SetNumberOfCharged(), G4Fragment::SetNumberOfHoles(), G4Fragment::SetNumberOfParticles(), G4VPreCompoundTransitions::TransitionProb1, G4VPreCompoundTransitions::TransitionProb2, and G4VPreCompoundTransitions::TransitionProb3.
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 }