#include <G4WilsonAblationModel.hh>
Inheritance diagram for G4WilsonAblationModel:
Public Types | |
typedef std::vector< G4ParticleDefinition * > | VectorOfFragmentTypes |
Public Member Functions | |
G4WilsonAblationModel () | |
~G4WilsonAblationModel () | |
G4FragmentVector * | BreakItUp (const G4Fragment &theNucleus) |
void | SetProduceSecondaries (G4bool) |
G4bool | GetProduceSecondaries () |
void | SetVerboseLevel (G4int) |
G4int | GetVerboseLevel () |
Definition at line 83 of file G4WilsonAblationModel.hh.
typedef std::vector<G4ParticleDefinition*> G4WilsonAblationModel::VectorOfFragmentTypes |
Definition at line 89 of file G4WilsonAblationModel.hh.
G4WilsonAblationModel::G4WilsonAblationModel | ( | ) |
Definition at line 105 of file G4WilsonAblationModel.cc.
References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4VEvaporationFactory::GetChannel(), G4He3::He3(), G4Neutron::Neutron(), G4VEvaporation::OPTxs, G4Proton::Proton(), G4Triton::Triton(), and G4VEvaporation::useSICB.
00106 { 00107 // 00108 // 00109 // Send message to stdout to advise that the G4Abrasion model is being used. 00110 // 00111 PrintWelcomeMessage(); 00112 // 00113 // 00114 // Set the default verbose level to 0 - no output. 00115 // 00116 verboseLevel = 0; 00117 // 00118 // 00119 // Set the binding energy per nucleon .... did I mention that this is a crude 00120 // model for nuclear de-excitation? 00121 // 00122 B = 10.0 * MeV; 00123 // 00124 // 00125 // It is possuble to switch off secondary particle production (other than the 00126 // final nuclear fragment). The default is on. 00127 // 00128 produceSecondaries = true; 00129 // 00130 // 00131 // Now we need to define the decay modes. We're using the G4Evaporation model 00132 // to help determine the kinematics of the decay. 00133 // 00134 nFragTypes = 6; 00135 fragType[0] = G4Alpha::Alpha(); 00136 fragType[1] = G4He3::He3(); 00137 fragType[2] = G4Triton::Triton(); 00138 fragType[3] = G4Deuteron::Deuteron(); 00139 fragType[4] = G4Proton::Proton(); 00140 fragType[5] = G4Neutron::Neutron(); 00141 // 00142 // 00143 // Set verboseLevel default to no output. 00144 // 00145 verboseLevel = 0; 00146 theChannelFactory = new G4EvaporationFactory(); 00147 theChannels = theChannelFactory->GetChannel(); 00148 // 00149 // 00150 // Set defaults for evaporation classes. These can be overridden by user 00151 // "set" methods. 00152 // 00153 OPTxs = 3; 00154 useSICB = false; 00155 fragmentVector = 0; 00156 }
G4WilsonAblationModel::~G4WilsonAblationModel | ( | ) |
Definition at line 159 of file G4WilsonAblationModel.cc.
00160 { 00161 if (theChannels != 0) theChannels = 0; 00162 if (theChannelFactory != 0) delete theChannelFactory; 00163 }
G4FragmentVector * G4WilsonAblationModel::BreakItUp | ( | const G4Fragment & | theNucleus | ) | [virtual] |
Implements G4VEvaporation.
Definition at line 167 of file G4WilsonAblationModel.cc.
References G4cout, G4endl, G4UniformRand, G4Fragment::GetA_asInt(), G4ParticleDefinition::GetBaryonNumber(), G4Fragment::GetExcitationEnergy(), G4Fragment::GetGroundStateMass(), G4ParticleTable::GetIonTable(), G4Fragment::GetMomentum(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGCharge(), G4Fragment::GetZ_asInt(), CLHEP::detail::n, G4Neutron::Neutron(), and G4Proton::Proton().
00168 { 00169 // 00170 // 00171 // Initilise the pointer to the G4FragmentVector used to return the information 00172 // about the breakup. 00173 // 00174 fragmentVector = new G4FragmentVector; 00175 fragmentVector->clear(); 00176 // 00177 // 00178 // Get the A, Z and excitation of the nucleus. 00179 // 00180 G4int A = theNucleus.GetA_asInt(); 00181 G4int Z = theNucleus.GetZ_asInt(); 00182 G4double ex = theNucleus.GetExcitationEnergy(); 00183 if (verboseLevel >= 2) 00184 { 00185 G4cout <<"oooooooooooooooooooooooooooooooooooooooo" 00186 <<"oooooooooooooooooooooooooooooooooooooooo" 00187 <<G4endl; 00188 G4cout.precision(6); 00189 G4cout <<"IN G4WilsonAblationModel" <<G4endl; 00190 G4cout <<"Initial prefragment A=" <<A 00191 <<", Z=" <<Z 00192 <<", excitation energy = " <<ex/MeV <<" MeV" 00193 <<G4endl; 00194 } 00195 // 00196 // 00197 // Check that there is a nucleus to speak of. It's possible there isn't one 00198 // or its just a proton or neutron. In either case, the excitation energy 00199 // (from the Lorentz vector) is not used. 00200 // 00201 if (A == 0) 00202 { 00203 if (verboseLevel >= 2) 00204 { 00205 G4cout <<"No nucleus to decay" <<G4endl; 00206 G4cout <<"oooooooooooooooooooooooooooooooooooooooo" 00207 <<"oooooooooooooooooooooooooooooooooooooooo" 00208 <<G4endl; 00209 } 00210 return fragmentVector; 00211 } 00212 else if (A == 1) 00213 { 00214 G4LorentzVector lorentzVector = theNucleus.GetMomentum(); 00215 lorentzVector.setE(lorentzVector.e()-ex+10.0*eV); 00216 if (Z == 0) 00217 { 00218 G4Fragment *fragment = new G4Fragment(lorentzVector,G4Neutron::Neutron()); 00219 fragmentVector->push_back(fragment); 00220 } 00221 else 00222 { 00223 G4Fragment *fragment = new G4Fragment(lorentzVector,G4Proton::Proton()); 00224 fragmentVector->push_back(fragment); 00225 } 00226 if (verboseLevel >= 2) 00227 { 00228 G4cout <<"Final fragment is in fact only a nucleon) :" <<G4endl; 00229 G4cout <<(*fragmentVector)[0] <<G4endl; 00230 G4cout <<"oooooooooooooooooooooooooooooooooooooooo" 00231 <<"oooooooooooooooooooooooooooooooooooooooo" 00232 <<G4endl; 00233 } 00234 return fragmentVector; 00235 } 00236 // 00237 // 00238 // Then the number of nucleons ablated (either as nucleons or light nuclear 00239 // fragments) is based on a simple argument for the binding energy per nucleon. 00240 // 00241 G4int DAabl = (G4int) (ex / B); 00242 if (DAabl > A) DAabl = A; 00243 // The following lines are no longer accurate given we now treat the final fragment 00244 // if (verboseLevel >= 2) 00245 // G4cout <<"Number of nucleons ejected = " <<DAabl <<G4endl; 00246 00247 // 00248 // 00249 // Determine the nuclear fragment from the ablation process by sampling the 00250 // Rudstam equation. 00251 // 00252 G4int AF = A - DAabl; 00253 G4int ZF = 0; 00254 if (AF > 0) 00255 { 00256 G4double AFd = (G4double) AF; 00257 G4double R = 11.8 / std::pow(AFd, 0.45); 00258 G4int minZ = Z - DAabl; 00259 if (minZ <= 0) minZ = 1; 00260 // 00261 // 00262 // Here we define an integral probability distribution based on the Rudstam 00263 // equation assuming a constant AF. 00264 // 00265 G4double sig[100]; 00266 G4double sum = 0.0; 00267 for (G4int ii=minZ; ii<= Z; ii++) 00268 { 00269 sum += std::exp(-R*std::pow(std::abs(ii - 0.486*AFd + 3.8E-04*AFd*AFd),1.5)); 00270 sig[ii] = sum; 00271 } 00272 // 00273 // 00274 // Now sample that distribution to determine a value for ZF. 00275 // 00276 G4double xi = G4UniformRand(); 00277 G4int iz = minZ; 00278 G4bool found = false; 00279 while (iz <= Z && !found) 00280 { 00281 found = (xi <= sig[iz]/sum); 00282 if (!found) iz++; 00283 } 00284 if (iz > Z) 00285 ZF = Z; 00286 else 00287 ZF = iz; 00288 } 00289 G4int DZabl = Z - ZF; 00290 // 00291 // 00292 // Now determine the nucleons or nuclei which have bee ablated. The preference 00293 // is for the production of alphas, then other nuclei in order of decreasing 00294 // binding energy. The energies assigned to the products of the decay are 00295 // provisional for the moment (the 10eV is just to avoid errors with negative 00296 // excitation energies due to rounding). 00297 // 00298 G4double totalEpost = 0.0; 00299 evapType.clear(); 00300 for (G4int ift=0; ift<nFragTypes; ift++) 00301 { 00302 G4ParticleDefinition *type = fragType[ift]; 00303 G4double n = std::floor((G4double) DAabl / type->GetBaryonNumber() + 1.0E-10); 00304 G4double n1 = 1.0E+10; 00305 if (fragType[ift]->GetPDGCharge() > 0.0) 00306 n1 = std::floor((G4double) DZabl / type->GetPDGCharge() + 1.0E-10); 00307 if (n > n1) n = n1; 00308 if (n > 0.0) 00309 { 00310 G4double mass = type->GetPDGMass(); 00311 for (G4int j=0; j<(G4int) n; j++) 00312 { 00313 totalEpost += mass; 00314 evapType.push_back(type); 00315 } 00316 DAabl -= (G4int) (n * type->GetBaryonNumber() + 1.0E-10); 00317 DZabl -= (G4int) (n * type->GetPDGCharge() + 1.0E-10); 00318 } 00319 } 00320 // 00321 // 00322 // Determine the properties of the final nuclear fragment. Note that if 00323 // the final fragment is predicted to have a nucleon number of zero, then 00324 // really it's the particle last in the vector evapType which becomes the 00325 // final fragment. Therefore delete this from the vector if this is the 00326 // case. 00327 // 00328 G4double massFinalFrag = 0.0; 00329 if (AF > 0) 00330 massFinalFrag = G4ParticleTable::GetParticleTable()->GetIonTable()-> 00331 GetIonMass(ZF,AF); 00332 else 00333 { 00334 G4ParticleDefinition *type = evapType[evapType.size()-1]; 00335 AF = type->GetBaryonNumber(); 00336 ZF = (G4int) (type->GetPDGCharge() + 1.0E-10); 00337 evapType.erase(evapType.end()-1); 00338 } 00339 totalEpost += massFinalFrag; 00340 // 00341 // 00342 // Provide verbose output on the nuclear fragment if requested. 00343 // 00344 if (verboseLevel >= 2) 00345 { 00346 G4cout <<"Final fragment A=" <<AF 00347 <<", Z=" <<ZF 00348 <<G4endl; 00349 for (G4int ift=0; ift<nFragTypes; ift++) 00350 { 00351 G4ParticleDefinition *type = fragType[ift]; 00352 G4int n = std::count(evapType.begin(),evapType.end(),type); 00353 if (n > 0) 00354 G4cout <<"Particle type: " <<std::setw(10) <<type->GetParticleName() 00355 <<", number of particles emitted = " <<n <<G4endl; 00356 } 00357 } 00358 // 00359 // Add the total energy from the fragment. Note that the fragment is assumed 00360 // to be de-excited and does not undergo photo-evaporation .... I did mention 00361 // this is a bit of a crude model? 00362 // 00363 G4double massPreFrag = theNucleus.GetGroundStateMass(); 00364 G4double totalEpre = massPreFrag + ex; 00365 G4double excess = totalEpre - totalEpost; 00366 // G4Fragment *resultNucleus(theNucleus); 00367 G4Fragment *resultNucleus = new G4Fragment(A, Z, theNucleus.GetMomentum()); 00368 G4ThreeVector boost(0.0,0.0,0.0); 00369 G4int nEvap = 0; 00370 if (produceSecondaries && evapType.size()>0) 00371 { 00372 if (excess > 0.0) 00373 { 00374 SelectSecondariesByEvaporation (resultNucleus); 00375 nEvap = fragmentVector->size(); 00376 boost = resultNucleus->GetMomentum().findBoostToCM(); 00377 if (evapType.size() > 0) 00378 SelectSecondariesByDefault (boost); 00379 } 00380 else 00381 SelectSecondariesByDefault(G4ThreeVector(0.0,0.0,0.0)); 00382 } 00383 00384 if (AF > 0) 00385 { 00386 G4double mass = G4ParticleTable::GetParticleTable()->GetIonTable()-> 00387 GetIonMass(ZF,AF); 00388 G4double e = mass + 10.0*eV; 00389 G4double p = std::sqrt(e*e-mass*mass); 00390 G4ThreeVector direction(0.0,0.0,1.0); 00391 G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e); 00392 lorentzVector.boost(-boost); 00393 G4Fragment* frag = new G4Fragment(AF, ZF, lorentzVector); 00394 fragmentVector->push_back(frag); 00395 } 00396 delete resultNucleus; 00397 // 00398 // 00399 // Provide verbose output on the ablation products if requested. 00400 // 00401 if (verboseLevel >= 2) 00402 { 00403 if (nEvap > 0) 00404 { 00405 G4cout <<"----------------------" <<G4endl; 00406 G4cout <<"Evaporated particles :" <<G4endl; 00407 G4cout <<"----------------------" <<G4endl; 00408 } 00409 G4int ie = 0; 00410 G4FragmentVector::iterator iter; 00411 for (iter = fragmentVector->begin(); iter != fragmentVector->end(); iter++) 00412 { 00413 if (ie == nEvap) 00414 { 00415 // G4cout <<*iter <<G4endl; 00416 G4cout <<"---------------------------------" <<G4endl; 00417 G4cout <<"Particles from default emission :" <<G4endl; 00418 G4cout <<"---------------------------------" <<G4endl; 00419 } 00420 G4cout <<*iter <<G4endl; 00421 } 00422 G4cout <<"oooooooooooooooooooooooooooooooooooooooo" 00423 <<"oooooooooooooooooooooooooooooooooooooooo" 00424 <<G4endl; 00425 } 00426 00427 return fragmentVector; 00428 }
G4bool G4WilsonAblationModel::GetProduceSecondaries | ( | ) | [inline] |
G4int G4WilsonAblationModel::GetVerboseLevel | ( | ) | [inline] |
void G4WilsonAblationModel::SetProduceSecondaries | ( | G4bool | ) | [inline] |
void G4WilsonAblationModel::SetVerboseLevel | ( | G4int | ) | [inline] |
Definition at line 126 of file G4WilsonAblationModel.hh.
Referenced by G4WilsonAbrasionModel::G4WilsonAbrasionModel(), G4WilsonAbrasionModel::SetUseAblation(), and G4WilsonAbrasionModel::SetVerboseLevel().