Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4WilsonAblationModel Class Reference

#include <G4WilsonAblationModel.hh>

+ Inheritance diagram for G4WilsonAblationModel:

Public Types

typedef std::vector< G4ParticleDefinition * > VectorOfFragmentTypes
 

Public Member Functions

 G4WilsonAblationModel ()
 
virtual ~G4WilsonAblationModel ()
 
G4FragmentVectorBreakItUp (const G4Fragment &theNucleus)
 
void SetProduceSecondaries (G4bool)
 
G4bool GetProduceSecondaries ()
 
void SetVerboseLevel (G4int)
 
G4int GetVerboseLevel ()
 
- Public Member Functions inherited from G4VEvaporation
 G4VEvaporation ()
 
virtual ~G4VEvaporation ()
 
virtual void BreakFragment (G4FragmentVector *, G4Fragment *theNucleus)
 
virtual void InitialiseChannels ()
 
virtual void SetPhotonEvaporation (G4VEvaporationChannel *ptr)
 
void SetFermiBreakUp (G4VFermiBreakUp *ptr)
 
G4VFermiBreakUpGetFermiBreakUp () const
 
G4VEvaporationChannelGetPhotonEvaporation ()
 
G4VEvaporationChannelGetFissionChannel ()
 
G4VEvaporationChannelGetChannel (size_t idx)
 
void SetOPTxs (G4int opt)
 
void UseSICB (G4bool use)
 
size_t GetNumberOfChannels () const
 
 G4VEvaporation (const G4VEvaporation &right)=delete
 
const G4VEvaporationoperator= (const G4VEvaporation &right)=delete
 
G4bool operator== (const G4VEvaporation &right) const =delete
 
G4bool operator!= (const G4VEvaporation &right) const =delete
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEvaporation
void CleanChannels ()
 
- Protected Attributes inherited from G4VEvaporation
G4VEvaporationChannelthePhotonEvaporation
 
G4VFermiBreakUptheFBU
 
G4int OPTxs
 
G4bool useSICB
 
std::vector< G4VEvaporationChannel * > * theChannels
 
G4VEvaporationFactorytheChannelFactory
 

Detailed Description

Definition at line 84 of file G4WilsonAblationModel.hh.

Member Typedef Documentation

◆ VectorOfFragmentTypes

Constructor & Destructor Documentation

◆ G4WilsonAblationModel()

G4WilsonAblationModel::G4WilsonAblationModel ( )

Definition at line 117 of file G4WilsonAblationModel.cc.

118{
119//
120//
121// Send message to stdout to advise that the G4Abrasion model is being used.
122//
123 PrintWelcomeMessage();
124//
125//
126// Set the default verbose level to 0 - no output.
127//
128 verboseLevel = 0;
129//
130//
131// Set the binding energy per nucleon .... did I mention that this is a crude
132// model for nuclear de-excitation?
133//
134 B = 10.0 * MeV;
135//
136//
137// It is possuble to switch off secondary particle production (other than the
138// final nuclear fragment). The default is on.
139//
140 produceSecondaries = true;
141//
142//
143// Now we need to define the decay modes. We're using the G4Evaporation model
144// to help determine the kinematics of the decay.
145//
146 nFragTypes = 6;
147 fragType[0] = G4Alpha::Alpha();
148 fragType[1] = G4He3::He3();
149 fragType[2] = G4Triton::Triton();
150 fragType[3] = G4Deuteron::Deuteron();
151 fragType[4] = G4Proton::Proton();
152 fragType[5] = G4Neutron::Neutron();
153 for(G4int i=0; i<200; ++i) { fSig[i] = 0.0; }
154//
155//
156// Set verboseLevel default to no output.
157//
158 verboseLevel = 0;
161//
162//
163// Set defaults for evaporation classes. These can be overridden by user
164// "set" methods.
165//
166 OPTxs = 3;
167 useSICB = false;
168 fragmentVector = 0;
169
170 secID = G4PhysicsModelCatalog::GetModelID("model_G4WilsonAblationModel");
171}
int G4int
Definition: G4Types.hh:85
static G4Alpha * Alpha()
Definition: G4Alpha.cc:88
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:93
static G4He3 * He3()
Definition: G4He3.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:103
static G4int GetModelID(const G4int modelIndex)
static G4Proton * Proton()
Definition: G4Proton.cc:92
static G4Triton * Triton()
Definition: G4Triton.cc:93
virtual std::vector< G4VEvaporationChannel * > * GetChannel()=0
G4VEvaporationFactory * theChannelFactory
std::vector< G4VEvaporationChannel * > * theChannels

◆ ~G4WilsonAblationModel()

G4WilsonAblationModel::~G4WilsonAblationModel ( )
virtual

Definition at line 174 of file G4WilsonAblationModel.cc.

175{}

Member Function Documentation

◆ BreakItUp()

G4FragmentVector * G4WilsonAblationModel::BreakItUp ( const G4Fragment theNucleus)

Definition at line 179 of file G4WilsonAblationModel.cc.

181{
182//
183//
184// Initilise the pointer to the G4FragmentVector used to return the information
185// about the breakup.
186//
187 fragmentVector = new G4FragmentVector;
188 fragmentVector->clear();
189//
190//
191// Get the A, Z and excitation of the nucleus.
192//
193 G4int A = theNucleus.GetA_asInt();
194 G4int Z = theNucleus.GetZ_asInt();
195 G4double ex = theNucleus.GetExcitationEnergy();
196 if (verboseLevel >= 2)
197 {
198 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
199 <<"oooooooooooooooooooooooooooooooooooooooo"
200 <<G4endl;
201 G4cout.precision(6);
202 G4cout <<"IN G4WilsonAblationModel" <<G4endl;
203 G4cout <<"Initial prefragment A=" <<A
204 <<", Z=" <<Z
205 <<", excitation energy = " <<ex/MeV <<" MeV"
206 <<G4endl;
207 }
208//
209//
210// Check that there is a nucleus to speak of. It's possible there isn't one
211// or its just a proton or neutron. In either case, the excitation energy
212// (from the Lorentz vector) is not used.
213//
214 if (A == 0)
215 {
216 if (verboseLevel >= 2)
217 {
218 G4cout <<"No nucleus to decay" <<G4endl;
219 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
220 <<"oooooooooooooooooooooooooooooooooooooooo"
221 <<G4endl;
222 }
223 return fragmentVector;
224 }
225 else if (A == 1)
226 {
227 G4LorentzVector lorentzVector = theNucleus.GetMomentum();
228 lorentzVector.setE(lorentzVector.e()-ex+10.0*eV);
229 if (Z == 0)
230 {
231 G4Fragment *fragment = new G4Fragment(lorentzVector,G4Neutron::Neutron());
232 if (fragment != nullptr) { fragment->SetCreatorModelID(secID); }
233 fragmentVector->push_back(fragment);
234 }
235 else
236 {
237 G4Fragment *fragment = new G4Fragment(lorentzVector,G4Proton::Proton());
238 if (fragment != nullptr) { fragment->SetCreatorModelID(secID); }
239 fragmentVector->push_back(fragment);
240 }
241 if (verboseLevel >= 2)
242 {
243 G4cout <<"Final fragment is in fact only a nucleon) :" <<G4endl;
244 G4cout <<(*fragmentVector)[0] <<G4endl;
245 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
246 <<"oooooooooooooooooooooooooooooooooooooooo"
247 <<G4endl;
248 }
249 return fragmentVector;
250 }
251//
252//
253// Then the number of nucleons ablated (either as nucleons or light nuclear
254// fragments) is based on a simple argument for the binding energy per nucleon.
255//
256 G4int DAabl = (G4int) (ex / B);
257 if (DAabl > A) DAabl = A;
258// The following lines are no longer accurate given we now treat the final fragment
259// if (verboseLevel >= 2)
260// G4cout <<"Number of nucleons ejected = " <<DAabl <<G4endl;
261
262//
263//
264// Determine the nuclear fragment from the ablation process by sampling the
265// Rudstam equation.
266//
267 G4int AF = A - DAabl;
268 G4int ZF = 0;
269
270 if (AF > 0)
271 {
272 G4Pow* g4calc = G4Pow::GetInstance();
273 G4double AFd = (G4double) AF;
274 G4double R = 11.8 / g4calc->powZ(AF, 0.45);
275 G4int minZ = std::max(1, Z - DAabl);
276//
277//
278// Here we define an integral probability distribution based on the Rudstam
279// equation assuming a constant AF.
280//
281 G4int zmax = std::min(199, Z);
282 G4double sum = 0.0;
283 for (ZF=minZ; ZF<=zmax; ++ZF)
284 {
285 sum += G4Exp(-R*g4calc->powA(std::abs(ZF - 0.486*AFd + 3.8E-04*AFd*AFd),1.5));
286 fSig[ZF] = sum;
287 }
288//
289//
290// Now sample that distribution to determine a value for ZF.
291//
292 sum *= G4UniformRand();
293 for (ZF=minZ; ZF<=zmax; ++ZF) {
294 if(sum <= fSig[ZF]) { break; }
295 }
296 }
297 G4int DZabl = Z - ZF;
298//
299//
300// Now determine the nucleons or nuclei which have bee ablated. The preference
301// is for the production of alphas, then other nuclei in order of decreasing
302// binding energy. The energies assigned to the products of the decay are
303// provisional for the moment (the 10eV is just to avoid errors with negative
304// excitation energies due to rounding).
305//
306 G4double totalEpost = 0.0;
307 evapType.clear();
308 for (G4int ift=0; ift<nFragTypes; ift++)
309 {
310 G4ParticleDefinition *type = fragType[ift];
311 G4double n = std::floor((G4double) DAabl / type->GetBaryonNumber() + 1.0E-10);
312 G4double n1 = 1.0E+10;
313 if (fragType[ift]->GetPDGCharge() > 0.0)
314 n1 = std::floor((G4double) DZabl / type->GetPDGCharge() + 1.0E-10);
315 if (n > n1) n = n1;
316 if (n > 0.0)
317 {
318 G4double mass = type->GetPDGMass();
319 for (G4int j=0; j<(G4int) n; j++)
320 {
321 totalEpost += mass;
322 evapType.push_back(type);
323 }
324 DAabl -= (G4int) (n * type->GetBaryonNumber() + 1.0E-10);
325 DZabl -= (G4int) (n * type->GetPDGCharge() + 1.0E-10);
326 }
327 }
328//
329//
330// Determine the properties of the final nuclear fragment. Note that if
331// the final fragment is predicted to have a nucleon number of zero, then
332// really it's the particle last in the vector evapType which becomes the
333// final fragment. Therefore delete this from the vector if this is the
334// case.
335//
336 G4double massFinalFrag = 0.0;
337 if (AF > 0)
339 GetIonMass(ZF,AF);
340 else
341 {
342 G4ParticleDefinition *type = evapType[evapType.size()-1];
343 AF = type->GetBaryonNumber();
344 ZF = (G4int) (type->GetPDGCharge() + 1.0E-10);
345 evapType.erase(evapType.end()-1);
346 }
347 totalEpost += massFinalFrag;
348//
349//
350// Provide verbose output on the nuclear fragment if requested.
351//
352 if (verboseLevel >= 2)
353 {
354 G4cout <<"Final fragment A=" <<AF
355 <<", Z=" <<ZF
356 <<G4endl;
357 for (G4int ift=0; ift<nFragTypes; ift++)
358 {
359 G4ParticleDefinition *type = fragType[ift];
360 G4long n = std::count(evapType.cbegin(),evapType.cend(),type);
361 if (n > 0)
362 G4cout <<"Particle type: " <<std::setw(10) <<type->GetParticleName()
363 <<", number of particles emitted = " <<n <<G4endl;
364 }
365 }
366//
367// Add the total energy from the fragment. Note that the fragment is assumed
368// to be de-excited and does not undergo photo-evaporation .... I did mention
369// this is a bit of a crude model?
370//
371 G4double massPreFrag = theNucleus.GetGroundStateMass();
372 G4double totalEpre = massPreFrag + ex;
373 G4double excess = totalEpre - totalEpost;
374// G4Fragment *resultNucleus(theNucleus);
375 G4Fragment *resultNucleus = new G4Fragment(A, Z, theNucleus.GetMomentum());
376 G4ThreeVector boost(0.0,0.0,0.0);
377 std::size_t nEvap = 0;
378 if (produceSecondaries && evapType.size()>0)
379 {
380 if (excess > 0.0)
381 {
382 SelectSecondariesByEvaporation (resultNucleus);
383 nEvap = fragmentVector->size();
384 boost = resultNucleus->GetMomentum().findBoostToCM();
385 if (evapType.size() > 0)
386 SelectSecondariesByDefault (boost);
387 }
388 else
389 SelectSecondariesByDefault(G4ThreeVector(0.0,0.0,0.0));
390 }
391
392 if (AF > 0)
393 {
395 GetIonMass(ZF,AF);
396 G4double e = mass + 10.0*eV;
397 G4double p = std::sqrt(e*e-mass*mass);
398 G4ThreeVector direction(0.0,0.0,1.0);
399 G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
400 lorentzVector.boost(-boost);
401 G4Fragment* frag = new G4Fragment(AF, ZF, lorentzVector);
402 if (frag != nullptr) { frag->SetCreatorModelID(secID); }
403 fragmentVector->push_back(frag);
404 }
405 delete resultNucleus;
406//
407//
408// Provide verbose output on the ablation products if requested.
409//
410 if (verboseLevel >= 2)
411 {
412 if (nEvap > 0)
413 {
414 G4cout <<"----------------------" <<G4endl;
415 G4cout <<"Evaporated particles :" <<G4endl;
416 G4cout <<"----------------------" <<G4endl;
417 }
418 std::size_t ie = 0;
419 for (auto iter = fragmentVector->cbegin();
420 iter != fragmentVector->cend(); ++iter)
421 {
422 if (ie == nEvap)
423 {
424// G4cout <<*iter <<G4endl;
425 G4cout <<"---------------------------------" <<G4endl;
426 G4cout <<"Particles from default emission :" <<G4endl;
427 G4cout <<"---------------------------------" <<G4endl;
428 }
429 G4cout <<*iter <<G4endl;
430 }
431 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
432 <<"oooooooooooooooooooooooooooooooooooooooo"
433 <<G4endl;
434 }
435
436 return fragmentVector;
437}
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:65
CLHEP::HepLorentzVector G4LorentzVector
CLHEP::Hep3Vector G4ThreeVector
double G4double
Definition: G4Types.hh:83
long G4long
Definition: G4Types.hh:87
const G4int Z[17]
const G4double A[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
HepLorentzVector & boost(double, double, double)
Hep3Vector findBoostToCM() const
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:317
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:312
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:322
void SetCreatorModelID(G4int value)
Definition: G4Fragment.hh:433
G4int GetZ_asInt() const
Definition: G4Fragment.hh:289
G4int GetA_asInt() const
Definition: G4Fragment.hh:284
G4double GetPDGCharge() const
const G4String & GetParticleName() const
G4IonTable * GetIonTable() const
static G4ParticleTable * GetParticleTable()
Definition: G4Pow.hh:49
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:225
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230

◆ GetProduceSecondaries()

G4bool G4WilsonAblationModel::GetProduceSecondaries ( )
inline

Definition at line 121 of file G4WilsonAblationModel.hh.

122 {return produceSecondaries;}

◆ GetVerboseLevel()

G4int G4WilsonAblationModel::GetVerboseLevel ( )
inline

Definition at line 129 of file G4WilsonAblationModel.hh.

130 {return verboseLevel;}

◆ SetProduceSecondaries()

void G4WilsonAblationModel::SetProduceSecondaries ( G4bool  produceSecondaries1)
inline

Definition at line 116 of file G4WilsonAblationModel.hh.

118 {produceSecondaries = produceSecondaries1;}

◆ SetVerboseLevel()

void G4WilsonAblationModel::SetVerboseLevel ( G4int  verboseLevel1)
inline

The documentation for this class was generated from the following files: