Geant4 10.7.0
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 ()
 
void SetOPTxs (G4int opt)
 
void UseSICB (G4bool use)
 
size_t GetNumberOfChannels () const
 

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 116 of file G4WilsonAblationModel.cc.

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

◆ ~G4WilsonAblationModel()

G4WilsonAblationModel::~G4WilsonAblationModel ( )
virtual

Definition at line 171 of file G4WilsonAblationModel.cc.

172{}

Member Function Documentation

◆ BreakItUp()

G4FragmentVector * G4WilsonAblationModel::BreakItUp ( const G4Fragment theNucleus)

Definition at line 176 of file G4WilsonAblationModel.cc.

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

121 {return produceSecondaries;}

◆ GetVerboseLevel()

G4int G4WilsonAblationModel::GetVerboseLevel ( )
inline

Definition at line 128 of file G4WilsonAblationModel.hh.

129 {return verboseLevel;}

◆ SetProduceSecondaries()

void G4WilsonAblationModel::SetProduceSecondaries ( G4bool  produceSecondaries1)
inline

Definition at line 115 of file G4WilsonAblationModel.hh.

117 {produceSecondaries = produceSecondaries1;}

◆ SetVerboseLevel()

void G4WilsonAblationModel::SetVerboseLevel ( G4int  verboseLevel1)
inline

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