Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ExcitationHandler.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4ExcitationHandler.cc 102726 2017-02-20 13:15:29Z gcosmo $
27 //
28 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara (May 1998)
30 //
31 //
32 // Modified:
33 // 30 June 1998 by V. Lara:
34 // -Modified the Transform method for use G4ParticleTable and
35 // therefore G4IonTable. It makes possible to convert all kind
36 // of fragments (G4Fragment) produced in deexcitation to
37 // G4DynamicParticle
38 // -It uses default algorithms for:
39 // Evaporation: G4Evaporation
40 // MultiFragmentation: G4StatMF
41 // Fermi Breakup model: G4FermiBreakUp
42 // 24 Jul 2008 by M. A. Cortes Giraldo:
43 // -Max Z,A for Fermi Break-Up turns to 9,17 by default
44 // -BreakItUp() reorganised and bug in Evaporation loop fixed
45 // -Transform() optimised
46 // (September 2008) by J. M. Quesada. External choices have been added for :
47 // -inverse cross section option (default OPTxs=3)
48 // -superimposed Coulomb barrier (if useSICB is set true, by default it is false)
49 // September 2009 by J. M. Quesada:
50 // -according to Igor Pshenichnov, SMM will be applied (just in case) only once.
51 // 27 Nov 2009 by V.Ivanchenko:
52 // -cleanup the logic, reduce number internal vectors, fixed memory leak.
53 // 11 May 2010 by V.Ivanchenko:
54 // -FermiBreakUp activated, used integer Z and A, used BreakUpFragment method for
55 // final photon deexcitation; used check on adundance of a fragment, decay
56 // unstable fragments with A <5
57 // 22 March 2011 by V.Ivanchenko: general cleanup and addition of a condition:
58 // products of Fermi Break Up cannot be further deexcited by this model
59 // 30 March 2011 by V.Ivanchenko removed private inline methods, moved Set methods
60 // to the source
61 // 23 January 2012 by V.Ivanchenko general cleanup including destruction of
62 // objects, propagate G4PhotonEvaporation pointer to G4Evaporation class and
63 // not delete it here
64 
65 #include "G4ExcitationHandler.hh"
66 #include "G4SystemOfUnits.hh"
67 #include "G4LorentzVector.hh"
68 #include "G4NistManager.hh"
69 #include "G4ParticleTable.hh"
70 #include "G4ParticleTypes.hh"
71 #include "G4Ions.hh"
72 
73 #include "G4VMultiFragmentation.hh"
74 #include "G4VFermiBreakUp.hh"
75 
76 #include "G4VEvaporation.hh"
77 #include "G4VEvaporationChannel.hh"
78 #include "G4Evaporation.hh"
79 #include "G4StatMF.hh"
80 #include "G4FermiBreakUp.hh"
81 #include "G4FermiBreakUpVI.hh"
82 #include "G4NuclearLevelData.hh"
83 #include "G4Pow.hh"
84 
86  : maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),
87  fVerbose(0),isInitialised(false),isEvapLocal(true)
88 {
89  theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
90  nist = G4NistManager::Instance();
91 
92  theMultiFragmentation = nullptr;
93  theFermiModel = nullptr;
95  theEvaporation = new G4Evaporation();
96  thePhotonEvaporation = theEvaporation->GetPhotonEvaporation();
97  theResults.reserve(60);
98  results.reserve(30);
99  theEvapList.reserve(30);
100  thePhotoEvapList.reserve(10);
101  SetParameters();
102  if(fVerbose > 0) { G4cout << "### New handler " << this << G4endl; }
103 }
104 
106 {
107  //G4cout << "### Delete handler " << this << G4endl;
108  delete theMultiFragmentation;
109  delete theFermiModel;
110  if(isEvapLocal) { delete theEvaporation; }
111 }
112 
113 void G4ExcitationHandler::SetParameters()
114 {
115  G4DeexPrecoParameters* param =
117  minEForMultiFrag = param->GetMinExPerNucleounForMF();
118  minExcitation = param->GetMinExcitation();
119  if(!theFermiModel) { theFermiModel = new G4FermiBreakUpVI(); }
120  theEvaporation->SetFermiBreakUp(theFermiModel);
121 }
122 
124 {
125  if(isInitialised) { return; }
126  if(fVerbose > 0) {
127  G4cout << "G4ExcitationHandler::Initialise() started " << this << G4endl;
128  }
130  isInitialised = true;
131  SetParameters();
132  theMultiFragmentation = new G4StatMF();
133  theFermiModel->Initialise();
134  theEvaporation->InitialiseChannels();
136  G4cout << "Number of de-excitation channels "
137  << theEvaporation->GetNumberOfChannels();
138  if(fVerbose > 0) { G4cout << " " << this; }
139  G4cout << G4endl;
140  }
141 }
142 
144 {
145  if(ptr && ptr != theEvaporation) {
146  delete theEvaporation;
147  theEvaporation = ptr;
148  thePhotonEvaporation = ptr->GetPhotonEvaporation();
149  theEvaporation->SetFermiBreakUp(theFermiModel);
150  isEvapLocal = flag;
151  }
152 }
153 
154 void
156 {
157  if(ptr && ptr != theMultiFragmentation) {
158  delete theMultiFragmentation;
159  theMultiFragmentation = ptr;
160  }
161 }
162 
164 {
165  if(ptr && ptr != theFermiModel) {
166  delete theFermiModel;
167  theFermiModel = ptr;
168  theEvaporation->SetFermiBreakUp(theFermiModel);
169  }
170 }
171 
172 void
174 {
175  if(ptr && ptr != thePhotonEvaporation) {
176  thePhotonEvaporation = ptr;
177  theEvaporation->SetPhotonEvaporation(ptr);
178  }
179 }
180 
182 {
183  G4Evaporation* evap = static_cast<G4Evaporation*>(theEvaporation);
184  if(!evap) { return; }
185  if(val == fEvaporation) {
186  evap->SetDefaultChannel();
187  } else if(val == fCombined) {
188  evap->SetCombinedChannel();
189  } else if(val == fGEM) {
190  evap->SetGEMChannel();
191  }
192  evap->InitialiseChannels();
194  G4cout << "Number of de-excitation channels is changed to "
195  << theEvaporation->GetNumberOfChannels();
196  if(fVerbose > 0) { G4cout << " " << this; }
197  G4cout << G4endl;
198  }
199 }
200 
203 {
204  // Variables existing until end of method
205  G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
206  if(fVerbose > 1) {
207  G4cout << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@ " << G4endl;
208  G4cout << theInitialState << G4endl;
209  }
210  if(!isInitialised) { Initialise(); }
211 
212  // pointer to fragment vector which receives temporal results
213  G4FragmentVector * theTempResult = nullptr;
214 
215  theResults.clear();
216  thePhotoEvapList.clear();
217  theEvapList.clear();
218 
219  // Variables to describe the excited configuration
220  G4double exEnergy = theInitialState.GetExcitationEnergy();
221  G4int A = theInitialState.GetA_asInt();
222  G4int Z = theInitialState.GetZ_asInt();
223 
224  // In case A <= 1 the fragment will not perform any nucleon emission
225  if (A <= 1) {
226  theResults.push_back( theInitialStatePtr );
227 
228  // check if a fragment is stable
229  } else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0) {
230  theResults.push_back( theInitialStatePtr );
231 
232  // JMQ 150909: first step in de-excitation is treated separately
233  // Fragments after the first step are stored in theEvapList
234  } else {
235  if((A<maxAForFermiBreakUp && Z<maxZForFermiBreakUp)
236  || exEnergy <= minEForMultiFrag*A) {
237  theEvapList.push_back(theInitialStatePtr);
238 
239  // Statistical Multifragmentation will take place only once
240  } else {
241  theTempResult = theMultiFragmentation->BreakItUp(theInitialState);
242  if(!theTempResult) {
243  theEvapList.push_back(theInitialStatePtr);
244  } else {
245  size_t nsec = theTempResult->size();
246 
247  // no fragmentation
248  if(0 == nsec) {
249  theEvapList.push_back(theInitialStatePtr);
250 
251  // secondary are produced - sort out secondary fragments
252  } else {
253  G4bool deletePrimary = true;
254  G4FragmentVector::iterator j;
255  for (j = theTempResult->begin(); j != theTempResult->end(); ++j) {
256 
257  if((*j) == theInitialStatePtr) { deletePrimary = false; }
258  A = (*j)->GetA_asInt();
259 
260  // gamma, p, n
261  if(A <= 1) {
262  theResults.push_back(*j);
263 
264  // Analyse fragment A > 1
265  } else {
266  G4double exEnergy1 = (*j)->GetExcitationEnergy();
267 
268  // cold fragments
269  if(exEnergy1 < minExcitation) {
270  Z = (*j)->GetZ_asInt();
271  if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
272  theResults.push_back(*j); // stable fragment
273  } else {
274  theEvapList.push_back(*j);
275  }
276  // hot fragments are unstable
277  } else {
278  theEvapList.push_back(*j);
279  }
280  }
281  }
282  if( deletePrimary ) { delete theInitialStatePtr; }
283  }
284  delete theTempResult; // end multifragmentation
285  }
286  }
287  }
288  if(fVerbose > 2) {
289  G4cout << "## After first step " << theEvapList.size() << " for evap; "
290  << thePhotoEvapList.size() << " for photo-evap; "
291  << theResults.size() << " results. " << G4endl;
292  }
293  // -----------------------------------
294  // FermiBreakUp and De-excitation loop
295  // -----------------------------------
296 
297  static const G4int countmax = 1000;
298  G4Fragment* frag;
299  size_t kk;
300  for (kk=0; kk<theEvapList.size(); ++kk) {
301  frag = theEvapList[kk];
302  if(fVerbose > 2) {
303  G4cout << "Next evaporate: " << G4endl;
304  G4cout << *frag << G4endl;
305  }
306  if(kk >= countmax) {
308  ed << "Infinite loop in the de-excitation module: " << kk
309  << " iterations \n"
310  << " Initial fragment: \n" << theInitialState
311  << "\n Current fragment: \n" << *frag;
312  G4Exception("G4ExcitationHandler::BreakItUp","had0333",FatalException,
313  ed,"Stop execution");
314 
315  }
316  A = frag->GetA_asInt();
317  Z = frag->GetZ_asInt();
318  results.clear();
319 
320  // Fermi Break-Up
321  if(theFermiModel->IsApplicable(Z, A, frag->GetExcitationEnergy())) {
322  theFermiModel->BreakFragment(&results, frag);
323  size_t nsec = results.size();
324  if(fVerbose > 2) { G4cout << "FermiBreakUp Nsec= " << nsec << G4endl; }
325 
326  // FBU takes care to delete input fragment or add it to the results
327  // The secondary may be excited - photo-evaporation should be applied
328  for(size_t j=0; j<nsec; ++j) {
329  exEnergy = results[j]->GetExcitationEnergy();
330  if(exEnergy < minExcitation) { theResults.push_back(results[j]); }
331  else { thePhotoEvapList.push_back(results[j]); }
332  }
333  continue;
334  }
335  // apply Evaporation, residual nucleus is always added to the results
336  theEvaporation->BreakFragment(&results, frag);
337  size_t nsec = results.size();
338  if(fVerbose > 2) { G4cout << "Evaporation Nsec= " << nsec << G4endl; }
339 
340  // no evaporation
341  if(1 >= nsec) {
342  theResults.push_back(frag);
343  continue;
344  }
345 
346  // Sort out secondary fragments
347  for (size_t j = 0; j<nsec; ++j) {
348  if(fVerbose > 3) {
349  G4cout << "Evaporated product #" << j << G4endl;
350  G4cout << results[j] << G4endl;
351  }
352  A = results[j]->GetA_asInt();
353  if(A <= 1) {
354  theResults.push_back(results[j]); // gamma, p, n
355  continue;
356  }
357  exEnergy = results[j]->GetExcitationEnergy();
358 
359  // hot fragment
360  if(exEnergy >= minExcitation) {
361  theEvapList.push_back(results[j]);
362 
363  // cold fragment
364  } else {
365  Z = results[j]->GetZ_asInt();
366 
367  // natural isotope
368  if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
369  theResults.push_back(results[j]); // stable fragment
370 
371  } else {
372  theEvapList.push_back(results[j]);
373  }
374  }
375  } // end of loop on secondary
376  } // end of the loop over theEvapList
377  if(fVerbose > 2) {
378  G4cout << "## After 2nd step " << theEvapList.size() << " was evap; "
379  << thePhotoEvapList.size() << " for photo-evap; "
380  << theResults.size() << " results. " << G4endl;
381  }
382  // -----------------------
383  // Photon-Evaporation loop
384  // -----------------------
385 
386  // at this point only photon evaporation is possible
387  size_t kkmax = thePhotoEvapList.size();
388  for (kk=0; kk<kkmax; ++kk) {
389  frag = thePhotoEvapList[kk];
390  if(fVerbose > 2) {
391  G4cout << "Next photon evaporate: " << thePhotonEvaporation << G4endl;
392  G4cout << *frag << G4endl;
393  }
394  exEnergy = frag->GetExcitationEnergy();
395 
396  // photon de-excitation only for hot fragments
397  if(exEnergy > minExcitation) {
398  thePhotonEvaporation->BreakUpChain(&theResults, frag);
399  }
400 
401  // primary fragment is kept
402  theResults.push_back(frag);
403 
404  } // end of photon-evaporation loop
405  if(fVerbose > 2) {
406  G4cout << "## After 3d step " << theEvapList.size() << " was evap; "
407  << thePhotoEvapList.size() << " was photo-evap; "
408  << theResults.size() << " results. " << G4endl;
409  }
410  G4ReactionProductVector * theReactionProductVector =
412 
413  // MAC (24/07/08)
414  // To optimise the storing speed, we reserve space in memory for the vector
415  theReactionProductVector->reserve( theResults.size() );
416 
417  G4int theFragmentA, theFragmentZ;
418 
419  if(fVerbose > 1) {
420  G4cout << "### ExcitationHandler provides " << theResults.size()
421  << " evaporated products:" << G4endl;
422  }
423  kkmax = theResults.size();
424  for (kk=0; kk<kkmax; ++kk) {
425  frag = theResults[kk];
426  if(fVerbose > 1) { G4cout << *frag << G4endl; }
427 
428  theFragmentA = frag->GetA_asInt();
429  theFragmentZ = frag->GetZ_asInt();
430  G4double etot= frag->GetMomentum().e();
431  G4double eexc = 0.0;
432  const G4ParticleDefinition* theKindOfFragment = nullptr;
433  if (theFragmentA == 0) { // photon or e-
434  theKindOfFragment = frag->GetParticleDefinition();
435  } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
436  theKindOfFragment = G4Neutron::NeutronDefinition();
437  } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
438  theKindOfFragment = G4Proton::ProtonDefinition();
439  } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
440  theKindOfFragment = G4Deuteron::DeuteronDefinition();
441  } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
442  theKindOfFragment = G4Triton::TritonDefinition();
443  } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
444  theKindOfFragment = G4He3::He3Definition();
445  } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
446  theKindOfFragment = G4Alpha::AlphaDefinition();;
447  } else {
448 
449  // fragment
450  eexc = frag->GetExcitationEnergy();
451  G4int idxf = frag->GetFloatingLevelNumber();
452  if(eexc < minExcitation) {
453  eexc = 0.0;
454  idxf = 0;
455  }
456 
457  theKindOfFragment = theTableOfIons->GetIon(theFragmentZ,theFragmentA,eexc,
458  G4Ions::FloatLevelBase(idxf));
459  if(fVerbose > 2) {
460  G4cout << "### EXCH: Find ion Z= " << theFragmentZ << " A= " << theFragmentA
461  << " Eexc(MeV)= " << eexc/MeV << " " << theKindOfFragment
462  << G4endl;
463  }
464  }
465  // fragment identified
466  if(theKindOfFragment) {
467  G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
468  theNew->SetMomentum(frag->GetMomentum().vect());
469  theNew->SetTotalEnergy(etot);
470  theNew->SetFormationTime(frag->GetCreationTime());
471  theReactionProductVector->push_back(theNew);
472 
473  // fragment not found out ground state is created
474  } else {
475  theKindOfFragment =
476  theTableOfIons->GetIon(theFragmentZ,theFragmentA,0.0,noFloat,0);
477  if(theKindOfFragment) {
478  G4ThreeVector mom(0.0,0.0,0.0);
479  G4double ionmass = theKindOfFragment->GetPDGMass();
480  if(etot <= ionmass) {
481  etot = ionmass;
482  } else {
483  G4double ptot = std::sqrt((etot - ionmass)*(etot + ionmass));
484  mom = (frag->GetMomentum().vect().unit())*ptot;
485  }
486  G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
487  theNew->SetMomentum(mom);
488  theNew->SetTotalEnergy(etot);
489  theNew->SetFormationTime(frag->GetCreationTime());
490  theReactionProductVector->push_back(theNew);
491  if(fVerbose > 2) {
492  G4cout << "### Find ion Z= " << theFragmentZ << " A= " << theFragmentA
493  << " ground state, energy corrected E(MeV)= " << etot << G4endl;
494  }
495  }
496  }
497  delete frag;
498  }
499  if(fVerbose > 2) {
500  G4cout << "@@@@@@@@@@ End G4Excitation Handler "<< G4endl;
501  }
502  return theReactionProductVector;
503 }
504 
505 void G4ExcitationHandler::ModelDescription(std::ostream& outFile) const
506 {
507  outFile << "G4ExcitationHandler description\n"
508  << "This class samples de-excitation of excited nucleus using\n"
509  << "Fermi Break-up model for light fragments (Z < 9, A < 17), "
510  << "evaporation, fission, and photo-evaporation models. Evaporated\n"
511  << "particle may be proton, neutron, and other light fragment \n"
512  << "(Z < 13, A < 29). During photon evaporation produced gamma \n"
513  << "or electrons due to internal conversion \n";
514 }
515 
516 
517 
518 
519 
520 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static G4Triton * TritonDefinition()
Definition: G4Triton.cc:90
void SetDefaultChannel()
G4int GetFloatingLevelNumber() const
Definition: G4Fragment.hh:427
static G4He3 * He3Definition()
Definition: G4He3.cc:89
void ModelDescription(std::ostream &outFile) const
virtual void BreakFragment(G4FragmentVector *results, G4Fragment *theNucleus)=0
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
virtual G4bool BreakUpChain(G4FragmentVector *theResult, G4Fragment *theNucleus)
G4VEvaporationChannel * GetPhotonEvaporation()
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
virtual void BreakFragment(G4FragmentVector *, G4Fragment *theNucleus)
const G4ParticleDefinition * GetParticleDefinition() const
Definition: G4Fragment.hh:438
virtual G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)=0
virtual G4bool IsApplicable(G4int Z, G4int A, G4double mass) const =0
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void SetGEMChannel()
G4IonTable * GetIonTable() const
static G4Ions::G4FloatLevelBase FloatLevelBase(char flbChar)
Definition: G4Ions.hh:189
Hep3Vector vect() const
virtual void InitialiseChannels() final
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:266
G4double GetCreationTime() const
Definition: G4Fragment.hh:448
bool G4bool
Definition: G4Types.hh:79
virtual void InitialiseChannels()
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:307
void SetFermiModel(G4VFermiBreakUp *ptr)
void SetTotalEnergy(const G4double en)
size_t GetNumberOfChannels() const
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
G4DeexPrecoParameters * GetParameters()
void SetMultiFragmentation(G4VMultiFragmentation *ptr)
G4double GetIsotopeAbundance(G4int Z, G4int N) const
void SetEvaporation(G4VEvaporation *ptr, G4bool isLocal=false)
void SetFermiBreakUp(G4VFermiBreakUp *ptr)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void SetCombinedChannel()
void SetDeexChannelsType(G4DeexChannelType val)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
virtual void SetPhotonEvaporation(G4VEvaporationChannel *ptr)
Hep3Vector unit() const
G4int GetZ_asInt() const
Definition: G4Fragment.hh:271
G4double GetMinExPerNucleounForMF() const
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
G4bool IsMasterThread()
Definition: G4Threading.cc:146
void SetFormationTime(G4double aTime)
virtual void Initialise()=0
double G4double
Definition: G4Types.hh:76
void SetPhotonEvaporation(G4VEvaporationChannel *ptr)
static G4Deuteron * DeuteronDefinition()
Definition: G4Deuteron.cc:89
static G4Alpha * AlphaDefinition()
Definition: G4Alpha.cc:84
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:99
static G4NuclearLevelData * GetInstance()
G4double GetMinExcitation() const
#define noFloat
Definition: G4Ions.hh:118
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:283