Geant4  10.02.p02
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 97619 2016-06-06 12:57:38Z 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 "G4FermiFragmentsPool.hh"
82 #include "G4Pow.hh"
83 
85  maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(400*GeV),
86  minExcitation(0.1*keV),OPTxs(3),useSICB(false),isEvapLocal(true)
87 {
90 
96  SetParameters();
98  theResults.reserve(60);
99  results.reserve(30);
100  theEvapList.reserve(30);
101  thePhotoEvapList.reserve(10);
102  //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 
115 {
116  //G4cout << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@" << G4endl;
117  // Variables existing until end of method
118  G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
119  //G4cout << theInitialState << G4endl;
120 
121  // pointer to fragment vector which receives temporal results
122  G4FragmentVector * theTempResult = 0;
123 
124  theResults.clear();
125  thePhotoEvapList.clear();
126  theEvapList.clear();
127 
128  // Variables to describe the excited configuration
129  G4double exEnergy = theInitialState.GetExcitationEnergy();
130  G4int A = theInitialState.GetA_asInt();
131  G4int Z = theInitialState.GetZ_asInt();
132 
133  // In case A <= 1 the fragment will not perform any nucleon emission
134  if (A <= 1) {
135  theResults.push_back( theInitialStatePtr );
136 
137  // check if a fragment is stable
138  } else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0) {
139  theResults.push_back( theInitialStatePtr );
140 
141  // JMQ 150909: first step in de-excitation is treated separately
142  // Fragments after the first step are stored in theEvapList
143  } else {
145  || exEnergy <= minEForMultiFrag*A) {
146  theEvapList.push_back(theInitialStatePtr);
147 
148  // Statistical Multifragmentation will take place only once
149  } else {
150  theTempResult = theMultiFragmentation->BreakItUp(theInitialState);
151  if(!theTempResult) {
152  theEvapList.push_back(theInitialStatePtr);
153  } else {
154  size_t nsec = theTempResult->size();
155 
156  // no fragmentation
157  if(0 == nsec) {
158  theEvapList.push_back(theInitialStatePtr);
159 
160  // secondary are produced - sort out secondary fragments
161  } else {
162  G4bool deletePrimary = true;
163  G4FragmentVector::iterator j;
164  for (j = theTempResult->begin(); j != theTempResult->end(); ++j) {
165  if((*j) == theInitialStatePtr) { deletePrimary = false; }
166  A = (*j)->GetA_asInt();
167 
168  // gamma, p, n
169  if(A <= 1) {
170  theResults.push_back(*j);
171 
172  // Analyse fragment A > 1
173  } else {
174  G4double exEnergy1 = (*j)->GetExcitationEnergy();
175 
176  // cold fragments
177  if(exEnergy1 < minExcitation) {
178  Z = (*j)->GetZ_asInt();
179  if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
180  theResults.push_back(*j); // stable fragment
181  } else {
182  theEvapList.push_back(*j);
183  }
184  // hot fragments are unstable
185  } else {
186  theEvapList.push_back(*j);
187  }
188  }
189  }
190  if( deletePrimary ) { delete theInitialStatePtr; }
191  }
192  delete theTempResult; // end multifragmentation
193  }
194  }
195  }
196  /*
197  G4cout << "## After first step " << theEvapList.size() << " for evap; "
198  << thePhotoEvapList.size() << " for photo-evap; "
199  << theResults.size() << " results. " << G4endl;
200  */
201  // -----------------------------------
202  // FermiBreakUp and De-excitation loop
203  // -----------------------------------
204 
205  static const G4int countmax = 1000;
206  G4Fragment* frag;
207  size_t kk;
208  for (kk=0; kk<theEvapList.size(); ++kk) {
209  frag = theEvapList[kk];
210  //G4cout << "Next evaporate: " << G4endl;
211  //G4cout << *frag << G4endl;
212  if(kk >= countmax) {
214  ed << "Infinite loop in the de-excitation module: " << kk
215  << " iterations \n"
216  << " Initial fragment: \n" << theInitialState
217  << "\n Current fragment: \n" << *frag;
218  G4Exception("G4ExcitationHandler::BreakItUp","had0333",FatalException,
219  ed,"Stop execution");
220 
221  }
222  A = frag->GetA_asInt();
223  Z = frag->GetZ_asInt();
224  results.clear();
225 
226  // Fermi Break-Up
227  if(A < maxAForFermiBreakUp && Z < maxZForFermiBreakUp && Z > 0 && A > Z) {
228  G4double etot = frag->GetExcitationEnergy() + frag->GetGroundStateMass();
229  if(thePool->IsApplicable(Z, A, etot)) {
231  size_t nsec = results.size();
232  //G4cout << "FermiBreakUp Nsec= " << nsec << G4endl;
233 
234  // FBU takes care to delete input fragment or add it to the results
235  // results may be excited - photo-evaporation should be applied
236  // If no final products then the fragment should be de-excited
237  // by evaporation
238  if(0 < nsec) {
239  for(size_t j=0; j<nsec; ++j) {
240  exEnergy = results[j]->GetExcitationEnergy();
241  if(exEnergy < minExcitation) { theResults.push_back(results[j]); }
242  else { thePhotoEvapList.push_back(results[j]); }
243  }
244  continue;
245  }
246  }
247  }
248  // apply Evaporation, residual nucleus is always added to the results
250  size_t nsec = results.size();
251  //G4cout << "Evaporation Nsec= " << nsec << G4endl;
252 
253  // no evaporation
254  if(1 >= nsec) {
255  theResults.push_back(frag);
256  continue;
257  }
258 
259  // Sort out secondary fragments
260  for (size_t j = 0; j<nsec; ++j) {
261  //G4cout << "Evaporated product #" << j << G4endl;
262  //G4cout << results[j] << G4endl;
263  A = results[j]->GetA_asInt();
264  //G4cout << "A= " << A << G4endl;
265  if(A <= 1) {
266  theResults.push_back(results[j]); // gamma, p, n
267  continue;
268  }
269  exEnergy = results[j]->GetExcitationEnergy();
270 
271  // hot fragment
272  if(exEnergy >= minExcitation) {
273  theEvapList.push_back(results[j]);
274 
275  // cold fragment
276  } else {
277  Z = results[j]->GetZ_asInt();
278 
279  // natural isotope
280  if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
281  theResults.push_back(results[j]); // stable fragment
282 
283  } else {
284  theEvapList.push_back(results[j]);
285  }
286  }
287  } // end of loop on secondary
288  } // end of the loop over theEvapList
289  /*
290  G4cout << "## After 2nd step " << theEvapList.size() << " was evap; "
291  << thePhotoEvapList.size() << " for photo-evap; "
292  << theResults.size() << " results. " << G4endl;
293  */
294  // -----------------------
295  // Photon-Evaporation loop
296  // -----------------------
297 
298  // at this point only photon evaporation is possible
299  size_t kkmax = thePhotoEvapList.size();
300  for (kk=0; kk<kkmax; ++kk) {
301  frag = thePhotoEvapList[kk];
302  //G4cout << "Next photon evaporate: " << thePhotonEvaporation << G4endl;
303  //G4cout << *frag << G4endl;
304  exEnergy = frag->GetExcitationEnergy();
305 
306  // photon de-excitation only for hot fragments
307  if(exEnergy > minExcitation) {
309  }
310 
311  // priamry fragment is kept
312  theResults.push_back(frag);
313 
314  } // end of photon-evaporation loop
315  /*
316  G4cout << "## After 3d step " << theEvapList.size() << " was evap; "
317  << thePhotoEvapList.size() << " was photo-evap; "
318  << theResults.size() << " results. " << G4endl;
319  */
320  G4ReactionProductVector * theReactionProductVector =
322 
323  // MAC (24/07/08)
324  // To optimise the storing speed, we reserve space in memory for the vector
325  theReactionProductVector->reserve( theResults.size() );
326 
327  G4int theFragmentA, theFragmentZ;
328 
329  //G4cout << "### ExcitationHandler provides " << theResults.size()
330  // << " evaporated products:" << G4endl;
331  kkmax = theResults.size();
332  for (kk=0; kk<kkmax; ++kk) {
333  frag = theResults[kk];
334  //G4cout << *frag << G4endl;
335 
336  theFragmentA = frag->GetA_asInt();
337  theFragmentZ = frag->GetZ_asInt();
338  G4double etot= frag->GetMomentum().e();
339  G4double eexc = 0.0;
340  const G4ParticleDefinition* theKindOfFragment = 0;
341  if (theFragmentA == 0) { // photon or e-
342  theKindOfFragment = frag->GetParticleDefinition();
343  } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
344  theKindOfFragment = G4Neutron::NeutronDefinition();
345  } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
346  theKindOfFragment = G4Proton::ProtonDefinition();
347  } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
348  theKindOfFragment = G4Deuteron::DeuteronDefinition();
349  } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
350  theKindOfFragment = G4Triton::TritonDefinition();
351  } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
352  theKindOfFragment = G4He3::He3Definition();
353  } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
354  theKindOfFragment = G4Alpha::AlphaDefinition();;
355  } else {
356 
357  // fragment
358  eexc = frag->GetExcitationEnergy();
359  if(eexc < minExcitation) { eexc = 0.0; }
360  theKindOfFragment = theTableOfIons->GetIon(theFragmentZ,theFragmentA,eexc);
361  /*
362  G4cout << "### Find ion Z= " << theFragmentZ << " A= " << theFragmentA
363  << " Eexc(MeV)= " << eexc/MeV << " " << theKindOfFragment
364  << G4endl;
365  */
366  }
367  // fragment identified
368  if(theKindOfFragment) {
369  G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
370  theNew->SetMomentum(frag->GetMomentum().vect());
371  theNew->SetTotalEnergy(etot);
372  theNew->SetFormationTime(frag->GetCreationTime());
373  theReactionProductVector->push_back(theNew);
374 
375  // fragment not found out ground state is created
376  } else {
377  theKindOfFragment = theTableOfIons->GetIon(theFragmentZ,theFragmentA,0.0);
378  if(theKindOfFragment) {
379  G4ThreeVector mom(0.0,0.0,0.0);
380  G4double ionmass = theKindOfFragment->GetPDGMass();
381  if(etot <= ionmass) {
382  etot = ionmass;
383  } else {
384  G4double ptot = std::sqrt((etot - ionmass)*(etot + ionmass));
385  mom = (frag->GetMomentum().vect().unit())*ptot;
386  }
387  G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
388  theNew->SetMomentum(mom);
389  theNew->SetTotalEnergy(etot);
390  theNew->SetFormationTime(frag->GetCreationTime());
391  theReactionProductVector->push_back(theNew);
392  /*
393  G4cout << "### Find ion Z= " << theFragmentZ << " A= " << theFragmentA
394  << " ground state, energy corrected " << theKindOfFragment << G4endl;
395  } else {
396  G4cout << "### Find ion Z= " << theFragmentZ
397  << " A= " << theFragmentA << " failed " << G4endl;
398  */
399  }
400  }
401  delete frag;
402  }
403  return theReactionProductVector;
404 }
405 
407 {
408  //for inverse cross section choice
410  //for the choice of superimposed Coulomb Barrier for inverse cross sections
413 }
414 
416 {
417  if(ptr && ptr != theEvaporation) {
418  delete theEvaporation;
419  theEvaporation = ptr;
421  SetParameters();
422  isEvapLocal = false;
423  }
424 }
425 
426 void
428 {
429  if(ptr && ptr != theMultiFragmentation) {
430  delete theMultiFragmentation;
431  theMultiFragmentation = ptr;
432  }
433 }
434 
436 {
437  if(ptr && ptr != theFermiModel) {
438  delete theFermiModel;
439  theFermiModel = ptr;
440  }
441 }
442 
443 void
445 {
446  if(ptr && ptr != thePhotonEvaporation) {
447  thePhotonEvaporation = ptr;
449  ptr->Initialise();
450  }
451 }
452 
454 {
455  maxZForFermiBreakUp = aZ;
456 }
457 
459 {
460  maxAForFermiBreakUp = anA;
461 }
462 
464 {
467 }
468 
470 {
471  minEForMultiFrag = anE;
472 }
473 void G4ExcitationHandler::ModelDescription(std::ostream& outFile) const
474 {
475  outFile << "G4ExcitationHandler description\n"
476  << "This class samples de-excitation of excited nucleus using\n"
477  << "Fermi Break-up model for light fragments (Z < 9, A < 17), "
478  << "evaporation, fission, and photo-evaporation models. Evaporated\n"
479  << "particle may be proton, neutron, and other light fragment \n"
480  << "(Z < 13, A < 29). During photon evaporation produced gamma \n"
481  << "or electrons due to internal conversion \n";
482 }
483 
484 
485 
486 
487 
488 
void SetOPTxs(G4int opt)
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static G4Triton * TritonDefinition()
Definition: G4Triton.cc:90
static G4He3 * He3Definition()
Definition: G4He3.cc:89
void ModelDescription(std::ostream &outFile) const
std::vector< G4Fragment * > theResults
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
CLHEP::Hep3Vector G4ThreeVector
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:491
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
virtual void BreakFragment(G4FragmentVector *, G4Fragment *theNucleus)
const G4ParticleDefinition * GetParticleDefinition() const
Definition: G4Fragment.hh:410
void SetMinEForMultiFrag(G4double anE)
virtual G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)=0
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState)
std::vector< G4Fragment * > results
std::vector< G4ReactionProduct * > G4ReactionProductVector
std::vector< G4Fragment * > thePhotoEvapList
G4bool IsApplicable(G4int Z, G4int A, G4double mass) const
G4IonTable * GetIonTable() const
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:256
G4double GetCreationTime() const
Definition: G4Fragment.hh:420
void UseSICB(G4bool use)
std::vector< G4Fragment * > theEvapList
bool G4bool
Definition: G4Types.hh:79
virtual void InitialiseChannels()
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:289
void SetFermiModel(G4VFermiBreakUp *ptr)
void SetTotalEnergy(const G4double en)
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
void SetMultiFragmentation(G4VMultiFragmentation *ptr)
G4VEvaporationChannel * thePhotonEvaporation
G4FermiFragmentsPool * thePool
static const double GeV
Definition: G4SIunits.hh:214
G4double GetIsotopeAbundance(G4int Z, G4int N) const
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:278
G4VMultiFragmentation * theMultiFragmentation
G4VFermiBreakUp * theFermiModel
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void BreakFragment(G4FragmentVector *, G4Fragment *theNucleus)=0
void SetMaxZForFermiBreakUp(G4int aZ)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
virtual void SetPhotonEvaporation(G4VEvaporationChannel *ptr)
void SetMaxAandZForFermiBreakUp(G4int anA, G4int aZ)
G4int GetZ_asInt() const
Definition: G4Fragment.hh:261
void SetEvaporation(G4VEvaporation *ptr)
G4VEvaporation * theEvaporation
void SetMaxAForFermiBreakUp(G4int anA)
static const double keV
Definition: G4SIunits.hh:213
void SetFormationTime(G4double aTime)
double G4double
Definition: G4Types.hh:76
static G4FermiFragmentsPool * Instance()
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
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:273