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