#include <G4EMDissociation.hh>
Inheritance diagram for G4EMDissociation:
Public Member Functions | |
G4EMDissociation () | |
G4EMDissociation (const G4EMDissociation &emd) | |
G4EMDissociation (G4ExcitationHandler *) | |
~G4EMDissociation () | |
const G4EMDissociation & | operator= (G4EMDissociation &right) |
virtual G4HadFinalState * | ApplyYourself (const G4HadProjectile &, G4Nucleus &) |
Definition at line 79 of file G4EMDissociation.hh.
G4EMDissociation::G4EMDissociation | ( | ) |
Definition at line 83 of file G4EMDissociation.cc.
References G4ExcitationHandler::SetEvaporation(), G4ExcitationHandler::SetFermiModel(), G4ExcitationHandler::SetMaxAandZForFermiBreakUp(), G4HadronicInteraction::SetMaxEnergy(), G4ExcitationHandler::SetMinEForMultiFrag(), G4HadronicInteraction::SetMinEnergy(), G4ExcitationHandler::SetMultiFragmentation(), and G4HadronicInteraction::verboseLevel.
00083 :G4HadronicInteraction("EMDissociation") { 00084 00085 // Send message to stdout to advise that the G4EMDissociation model is being 00086 // used. 00087 PrintWelcomeMessage(); 00088 00089 // No de-excitation handler has been supplied - define the default handler. 00090 theExcitationHandler = new G4ExcitationHandler; 00091 G4Evaporation* theEvaporation = new G4Evaporation; 00092 G4FermiBreakUp* theFermiBreakUp = new G4FermiBreakUp; 00093 G4StatMF* theMF = new G4StatMF; 00094 theExcitationHandler->SetEvaporation(theEvaporation); 00095 theExcitationHandler->SetFermiModel(theFermiBreakUp); 00096 theExcitationHandler->SetMultiFragmentation(theMF); 00097 theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6); 00098 theExcitationHandler->SetMinEForMultiFrag(5.0*MeV); 00099 handlerDefinedInternally = true; 00100 00101 // This EM dissociation model needs access to the cross-sections held in 00102 // G4EMDissociationCrossSection. 00103 dissociationCrossSection = new G4EMDissociationCrossSection; 00104 thePhotonSpectrum = new G4EMDissociationSpectrum; 00105 00106 // Set the minimum and maximum range for the model (despite nomanclature, this 00107 // is in energy per nucleon number). 00108 SetMinEnergy(100.0*MeV); 00109 SetMaxEnergy(500.0*GeV); 00110 00111 // Set the default verbose level to 0 - no output. 00112 verboseLevel = 0; 00113 }
G4EMDissociation::G4EMDissociation | ( | const G4EMDissociation & | emd | ) |
G4EMDissociation::G4EMDissociation | ( | G4ExcitationHandler * | ) |
Definition at line 137 of file G4EMDissociation.cc.
References G4HadronicInteraction::SetMaxEnergy(), G4HadronicInteraction::SetMinEnergy(), and G4HadronicInteraction::verboseLevel.
00138 { 00139 // 00140 // 00141 // Send message to stdout to advise that the G4EMDissociation model is being 00142 // used. 00143 // 00144 PrintWelcomeMessage(); 00145 00146 theExcitationHandler = aExcitationHandler; 00147 handlerDefinedInternally = false; 00148 // 00149 // 00150 // This EM dissociation model needs access to the cross-sections held in 00151 // G4EMDissociationCrossSection. 00152 // 00153 dissociationCrossSection = new G4EMDissociationCrossSection; 00154 thePhotonSpectrum = new G4EMDissociationSpectrum; 00155 // 00156 // 00157 // Set the minimum and maximum range for the model (despite nomanclature, this 00158 // is in energy per nucleon number). 00159 // 00160 SetMinEnergy(100.0*MeV); 00161 SetMaxEnergy(500.0*GeV); 00162 // 00163 // 00164 // Set the default verbose level to 0 - no output. 00165 // 00166 verboseLevel = 0; 00167 }
G4EMDissociation::~G4EMDissociation | ( | ) |
Definition at line 170 of file G4EMDissociation.cc.
00170 { 00171 if (handlerDefinedInternally) delete theExcitationHandler; 00172 // delete dissociationCrossSection; 00173 // Cross section deleted by G4CrossSectionRegistry; don't do it here 00174 // Bug reported by Gong Ding in Bug Report #1339 00175 delete thePhotonSpectrum; 00176 }
G4HadFinalState * G4EMDissociation::ApplyYourself | ( | const G4HadProjectile & | , | |
G4Nucleus & | ||||
) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 180 of file G4EMDissociation.cc.
References G4HadFinalState::Clear(), G4DynamicParticle::DumpInfo(), G4cout, G4endl, G4UniformRand, G4DynamicParticle::Get4Momentum(), G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4ParticleDefinition::GetBaryonNumber(), G4HadProjectile::GetDefinition(), G4HadProjectile::GetKineticEnergy(), G4NucleiProperties::GetNuclearMass(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4HadProjectile::GetTotalEnergy(), G4Nucleus::GetZ_asInt(), G4Neutron::NeutronDefinition(), G4INCL::Math::pi, G4InuclParticleNames::pp, G4Proton::ProtonDefinition(), G4DynamicParticle::Set4Momentum(), and stopAndKill.
00181 { 00182 // 00183 // 00184 // The secondaries will be returned in G4HadFinalState &theParticleChange - 00185 // initialise this. 00186 // 00187 theParticleChange.Clear(); 00188 theParticleChange.SetStatusChange(stopAndKill); 00189 // 00190 // 00191 // Get relevant information about the projectile and target (A, Z) and 00192 // energy/nuc, momentum, velocity, Lorentz factor and rest-mass of the 00193 // projectile. 00194 // 00195 const G4ParticleDefinition *definitionP = theTrack.GetDefinition(); 00196 const G4double AP = definitionP->GetBaryonNumber(); 00197 const G4double ZP = definitionP->GetPDGCharge(); 00198 G4LorentzVector pP = theTrack.Get4Momentum(); 00199 G4double E = theTrack.GetKineticEnergy()/AP; 00200 G4double MP = theTrack.GetTotalEnergy() - E*AP; 00201 G4double b = pP.beta(); 00202 G4double AT = theTarget.GetA_asInt(); 00203 G4double ZT = theTarget.GetZ_asInt(); 00204 G4double MT = G4NucleiProperties::GetNuclearMass(AT,ZT); 00205 // 00206 // 00207 // Depending upon the verbosity level, output the initial information on the 00208 // projectile and target. 00209 // 00210 if (verboseLevel >= 2) 00211 { 00212 G4cout.precision(6); 00213 G4cout <<"########################################" 00214 <<"########################################" 00215 <<G4endl; 00216 G4cout <<"IN G4EMDissociation" <<G4endl; 00217 G4cout <<"Initial projectile A=" <<AP 00218 <<", Z=" <<ZP 00219 <<G4endl; 00220 G4cout <<"Initial target A=" <<AT 00221 <<", Z=" <<ZT 00222 <<G4endl; 00223 G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl; 00224 } 00225 // 00226 // 00227 // Initialise the variables which will be used with the phase-space decay and 00228 // to boost the secondaries from the interaction. 00229 // 00230 G4ParticleDefinition *typeNucleon = NULL; 00231 G4ParticleDefinition *typeDaughter = NULL; 00232 G4double Eg = 0.0; 00233 G4double mass = 0.0; 00234 G4ThreeVector boost = G4ThreeVector(0.0, 0.0, 0.0); 00235 // 00236 // 00237 // Determine the cross-sections at the giant dipole and giant quadrupole 00238 // resonance energies for the projectile and then target. The information is 00239 // initially provided in the G4PhysicsFreeVector individually for the E1 00240 // and E2 fields. These are then summed. 00241 // 00242 G4double bmin = thePhotonSpectrum->GetClosestApproach(AP, ZP, AT, ZT, b); 00243 G4PhysicsFreeVector *crossSectionP = dissociationCrossSection-> 00244 GetCrossSectionForProjectile(AP, ZP, AT, ZT, b, bmin); 00245 G4PhysicsFreeVector *crossSectionT = dissociationCrossSection-> 00246 GetCrossSectionForTarget(AP, ZP, AT, ZT, b, bmin); 00247 00248 G4double totCrossSectionP = (*crossSectionP)[0]+(*crossSectionP)[1]; 00249 G4double totCrossSectionT = (*crossSectionT)[0]+(*crossSectionT)[1]; 00250 // 00251 // 00252 // Now sample whether the interaction involved EM dissociation of the projectile 00253 // or the target. 00254 // 00255 if (G4UniformRand() < 00256 totCrossSectionP / (totCrossSectionP + totCrossSectionT)) 00257 { 00258 // 00259 // 00260 // It was the projectile which underwent EM dissociation. Define the Lorentz 00261 // boost to be applied to the secondaries, and sample whether a proton or a 00262 // neutron was ejected. Then determine the energy of the virtual gamma ray 00263 // which passed from the target nucleus ... this will be used to define the 00264 // excitation of the projectile. 00265 // 00266 mass = MP; 00267 if (G4UniformRand() < dissociationCrossSection-> 00268 GetWilsonProbabilityForProtonDissociation (AP, ZP)) 00269 { 00270 if (verboseLevel >= 2) 00271 G4cout <<"Projectile underwent EM dissociation producing a proton" 00272 <<G4endl; 00273 typeNucleon = G4Proton::ProtonDefinition(); 00274 typeDaughter = G4ParticleTable::GetParticleTable()-> 00275 GetIon((G4int) ZP-1, (G4int) AP-1, 0.0); 00276 } 00277 else 00278 { 00279 if (verboseLevel >= 2) 00280 G4cout <<"Projectile underwent EM dissociation producing a neutron" 00281 <<G4endl; 00282 typeNucleon = G4Neutron::NeutronDefinition(); 00283 typeDaughter = G4ParticleTable::GetParticleTable()-> 00284 GetIon((G4int) ZP, (G4int) AP-1, 0.0); 00285 } 00286 if (G4UniformRand() < (*crossSectionP)[0]/totCrossSectionP) 00287 { 00288 Eg = crossSectionP->GetLowEdgeEnergy(0); 00289 if (verboseLevel >= 2) 00290 G4cout <<"Transition type was E1" <<G4endl; 00291 } 00292 else 00293 { 00294 Eg = crossSectionP->GetLowEdgeEnergy(1); 00295 if (verboseLevel >= 2) 00296 G4cout <<"Transition type was E2" <<G4endl; 00297 } 00298 // 00299 // 00300 // We need to define a Lorentz vector with the original momentum, but total 00301 // energy includes the projectile and virtual gamma. This is then used 00302 // to calculate the boost required for the secondaries. 00303 // 00304 pP.setE(pP.e()+Eg); 00305 boost = pP.findBoostToCM(); 00306 } 00307 else 00308 { 00309 // 00310 // 00311 // It was the target which underwent EM dissociation. Sample whether a 00312 // proton or a neutron was ejected. Then determine the energy of the virtual 00313 // gamma ray which passed from the projectile nucleus ... this will be used to 00314 // define the excitation of the target. 00315 // 00316 mass = MT; 00317 if (G4UniformRand() < dissociationCrossSection-> 00318 GetWilsonProbabilityForProtonDissociation (AT, ZT)) 00319 { 00320 if (verboseLevel >= 2) 00321 G4cout <<"Target underwent EM dissociation producing a proton" 00322 <<G4endl; 00323 typeNucleon = G4Proton::ProtonDefinition(); 00324 typeDaughter = G4ParticleTable::GetParticleTable()-> 00325 GetIon((G4int) ZT-1, (G4int) AT-1, 0.0); 00326 } 00327 else 00328 { 00329 if (verboseLevel >= 2) 00330 G4cout <<"Target underwent EM dissociation producing a neutron" 00331 <<G4endl; 00332 typeNucleon = G4Neutron::NeutronDefinition(); 00333 typeDaughter = G4ParticleTable::GetParticleTable()-> 00334 GetIon((G4int) ZT, (G4int) AT-1, 0.0); 00335 } 00336 if (G4UniformRand() < (*crossSectionT)[0]/totCrossSectionT) 00337 { 00338 Eg = crossSectionT->GetLowEdgeEnergy(0); 00339 if (verboseLevel >= 2) 00340 G4cout <<"Transition type was E1" <<G4endl; 00341 } 00342 else 00343 { 00344 Eg = crossSectionT->GetLowEdgeEnergy(1); 00345 if (verboseLevel >= 2) 00346 G4cout <<"Transition type was E2" <<G4endl; 00347 } 00348 // 00349 // 00350 // Add the projectile to theParticleChange, less the energy of the 00351 // not-so-virtual gamma-ray. Not that at the moment, no lateral momentum 00352 // is transferred between the projectile and target nuclei. 00353 // 00354 G4ThreeVector v = pP.vect(); 00355 v.setMag(1.0); 00356 G4DynamicParticle *changedP = new G4DynamicParticle 00357 (const_cast<G4ParticleDefinition*>(definitionP), v, E*AP-Eg); 00358 theParticleChange.AddSecondary (changedP); 00359 if (verboseLevel >= 2) 00360 { 00361 G4cout <<"Projectile change:" <<G4endl; 00362 changedP->DumpInfo(); 00363 } 00364 } 00365 // 00366 // 00367 // Perform a two-body decay based on the restmass energy of the parent and 00368 // gamma-ray, and the masses of the daughters. In the frame of reference of 00369 // the nucles, the angular distribution is sampled isotropically, but the 00370 // the nucleon and secondary nucleus are boosted if they've come from the 00371 // projectile. 00372 // 00373 G4double e = mass + Eg; 00374 G4double mass1 = typeNucleon->GetPDGMass(); 00375 G4double mass2 = typeDaughter->GetPDGMass(); 00376 G4double pp = (e+mass1+mass2)*(e+mass1-mass2)* 00377 (e-mass1+mass2)*(e-mass1-mass2)/(4.0*e*e); 00378 if (pp < 0.0) 00379 { 00380 pp = 1.0*eV; 00381 // if (verboseLevel >`= 1) 00382 // { 00383 // G4cout <<"IN G4EMDissociation::ApplyYoursef" <<G4endl; 00384 // G4cout <<"Error in mass of secondaries compared with primary:" <<G4endl; 00385 // G4cout <<"Rest mass of primary = " <<mass <<" MeV" <<G4endl; 00386 // G4cout <<"Virtual gamma energy = " <<Eg <<" MeV" <<G4endl; 00387 // G4cout <<"Rest mass of secondary #1 = " <<mass1 <<" MeV" <<G4endl; 00388 // G4cout <<"Rest mass of secondary #2 = " <<mass2 <<" MeV" <<G4endl; 00389 // } 00390 } 00391 else 00392 pp = std::sqrt(pp); 00393 G4double costheta = 2.*G4UniformRand()-1.0; 00394 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta)); 00395 G4double phi = 2.0*pi*G4UniformRand()*rad; 00396 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta); 00397 G4DynamicParticle *dynamicNucleon = 00398 new G4DynamicParticle(typeNucleon, direction*pp); 00399 dynamicNucleon->Set4Momentum(dynamicNucleon->Get4Momentum().boost(-boost)); 00400 G4DynamicParticle *dynamicDaughter = 00401 new G4DynamicParticle(typeDaughter, -direction*pp); 00402 dynamicDaughter->Set4Momentum(dynamicDaughter->Get4Momentum().boost(-boost)); 00403 // 00404 // 00405 // The "decay" products have to be transferred to the G4HadFinalState object. 00406 // Furthermore, the residual nucleus should be de-excited. 00407 // 00408 theParticleChange.AddSecondary (dynamicNucleon); 00409 if (verboseLevel >= 2) 00410 { 00411 G4cout <<"Nucleon from the EMD process:" <<G4endl; 00412 dynamicNucleon->DumpInfo(); 00413 } 00414 00415 G4Fragment *theFragment = new 00416 G4Fragment((G4int) typeDaughter->GetBaryonNumber(), 00417 (G4int) typeDaughter->GetPDGCharge(), dynamicDaughter->Get4Momentum()); 00418 if (verboseLevel >= 2) 00419 { 00420 G4cout <<"Dynamic properties of the prefragment:" <<G4endl; 00421 G4cout.precision(6); 00422 dynamicDaughter->DumpInfo(); 00423 G4cout <<"Nuclear properties of the prefragment:" <<G4endl; 00424 G4cout <<theFragment <<G4endl; 00425 } 00426 G4ReactionProductVector *products = 00427 theExcitationHandler->BreakItUp(*theFragment); 00428 delete theFragment; 00429 theFragment = NULL; 00430 00431 G4ReactionProductVector::iterator iter; 00432 for (iter = products->begin(); iter != products->end(); ++iter) 00433 { 00434 G4DynamicParticle *secondary = 00435 new G4DynamicParticle((*iter)->GetDefinition(), 00436 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum()); 00437 theParticleChange.AddSecondary (secondary); 00438 } 00439 00440 if (verboseLevel >= 2) 00441 G4cout <<"########################################" 00442 <<"########################################" 00443 <<G4endl; 00444 00445 return &theParticleChange; 00446 }
const G4EMDissociation& G4EMDissociation::operator= | ( | G4EMDissociation & | right | ) |