Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4WilsonAblationModel Class Reference

#include <G4WilsonAblationModel.hh>

Inheritance diagram for G4WilsonAblationModel:
Collaboration 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

Constructor & Destructor Documentation

G4WilsonAblationModel::G4WilsonAblationModel ( )

Definition at line 113 of file G4WilsonAblationModel.cc.

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

Here is the call graph for this function:

G4WilsonAblationModel::~G4WilsonAblationModel ( )
virtual

Definition at line 168 of file G4WilsonAblationModel.cc.

169 {}

Member Function Documentation

G4FragmentVector * G4WilsonAblationModel::BreakItUp ( const G4Fragment theNucleus)

Definition at line 174 of file G4WilsonAblationModel.cc.

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

Here is the call graph for this function:

G4bool G4WilsonAblationModel::GetProduceSecondaries ( )
inline

Definition at line 120 of file G4WilsonAblationModel.hh.

121  {return produceSecondaries;}
G4int G4WilsonAblationModel::GetVerboseLevel ( )
inline

Definition at line 128 of file G4WilsonAblationModel.hh.

129  {return verboseLevel;}
void G4WilsonAblationModel::SetProduceSecondaries ( G4bool  produceSecondaries1)
inline

Definition at line 116 of file G4WilsonAblationModel.hh.

117  {produceSecondaries = produceSecondaries1;}
void G4WilsonAblationModel::SetVerboseLevel ( G4int  verboseLevel1)
inline

Definition at line 124 of file G4WilsonAblationModel.hh.

125  {verboseLevel = verboseLevel1;}

Here is the caller graph for this function:


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