Geant4  10.02.p01
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 94381 2015-11-13 10:17:06Z 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.resize(60,0);
99  results.resize(30,0);
100  theEvapList.resize(30,0);
101  thePhotoEvapList.resize(10,0);
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  std::vector<G4Fragment*>::iterator iList;
206  for (iList = theEvapList.begin(); iList != theEvapList.end(); ++iList) {
207  //G4cout << "Next evaporate: " << G4endl;
208  //G4cout << *iList << G4endl;
209  G4Fragment* frag = *iList;
210  A = frag->GetA_asInt();
211  Z = frag->GetZ_asInt();
212  results.clear();
213 
214  // Fermi Break-Up
215  if(A < maxAForFermiBreakUp && Z < maxZForFermiBreakUp && Z > 0 && A > Z) {
216  G4double etot = frag->GetExcitationEnergy() + frag->GetGroundStateMass();
217  if(thePool->IsApplicable(Z, A, etot)) {
219  size_t nsec = results.size();
220  //G4cout << "FermiBreakUp Nsec= " << nsec << G4endl;
221 
222  // FBU takes care to delete input fragment or add it to the results
223  // results may be excited - photo-evaporation should be applied
224  // If no final products then the fragment should be de-excited
225  // by evaporation
226  if(0 < nsec) {
227  for(size_t j=0; j<nsec; ++j) {
228  exEnergy = results[j]->GetExcitationEnergy();
229  if(exEnergy < minExcitation) { theResults.push_back(results[j]); }
230  else { thePhotoEvapList.push_back(results[j]); }
231  }
232  continue;
233  }
234  }
235  }
236  // apply Evaporation, residual nucleus is always added to the results
238  size_t nsec = results.size();
239  //G4cout << "Evaporation Nsec= " << nsec << G4endl;
240 
241  // no evaporation
242  if(1 >= nsec) {
243  theResults.push_back(frag);
244  continue;
245  }
246 
247  // Sort out secondary fragments
248  for (size_t j = 0; j<nsec; ++j) {
249  //G4cout << "Evaporated product #" << j << G4endl;
250  //G4cout << results[j] << G4endl;
251  A = results[j]->GetA_asInt();
252  //G4cout << "A= " << A << G4endl;
253  if(A <= 1) {
254  theResults.push_back(results[j]); // gamma, p, n
255  continue;
256  }
257  exEnergy = results[j]->GetExcitationEnergy();
258 
259  // hot fragment
260  if(exEnergy >= minExcitation) {
261  theEvapList.push_back(results[j]);
262 
263  // cold fragment
264  } else {
265  Z = results[j]->GetZ_asInt();
266 
267  // natural isotope
268  if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
269  theResults.push_back(results[j]); // stable fragment
270 
271  } else {
272  theEvapList.push_back(results[j]);
273  }
274  }
275  } // end of loop on secondary
276  } // end of the loop over theEvapList
277  /*
278  G4cout << "## After 2nd step " << theEvapList.size() << " was evap; "
279  << thePhotoEvapList.size() << " for photo-evap; "
280  << theResults.size() << " results. " << G4endl;
281  */
282  // -----------------------
283  // Photon-Evaporation loop
284  // -----------------------
285 
286  // at this point only photon evaporation is possible
287  for(iList = thePhotoEvapList.begin(); iList != thePhotoEvapList.end(); ++iList) {
288  //G4cout << "Next photon evaporate: " << thePhotonEvaporation << G4endl;
289  //G4cout << *iList << G4endl;
290  exEnergy = (*iList)->GetExcitationEnergy();
291 
292  // photon de-excitation only for hot fragments
293  if(exEnergy > minExcitation) {
295  }
296 
297  // priamry fragment is kept
298  theResults.push_back(*iList);
299 
300  } // end of photon-evaporation loop
301  /*
302  G4cout << "## After 3d step " << theEvapList.size() << " was evap; "
303  << thePhotoEvapList.size() << " was photo-evap; "
304  << theResults.size() << " results. " << G4endl;
305  */
306  G4ReactionProductVector * theReactionProductVector =
308 
309  // MAC (24/07/08)
310  // To optimise the storing speed, we reserve space in memory for the vector
311  theReactionProductVector->reserve( theResults.size() );
312 
313  G4int theFragmentA, theFragmentZ;
314 
315  //G4cout << "### ExcitationHandler provides " << theResults.size()
316  // << " evaporated products:" << G4endl;
317  for (iList = theResults.begin(); iList != theResults.end(); ++iList) {
318  //G4cout << (*iList) << G4endl;
319 
320  theFragmentA = (*iList)->GetA_asInt();
321  theFragmentZ = (*iList)->GetZ_asInt();
322  G4double etot= (*iList)->GetMomentum().e();
323  G4double eexc = 0.0;
324  const G4ParticleDefinition* theKindOfFragment = 0;
325  if (theFragmentA == 0) { // photon or e-
326  theKindOfFragment = (*iList)->GetParticleDefinition();
327  } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
328  theKindOfFragment = G4Neutron::NeutronDefinition();
329  } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
330  theKindOfFragment = G4Proton::ProtonDefinition();
331  } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
332  theKindOfFragment = G4Deuteron::DeuteronDefinition();
333  } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
334  theKindOfFragment = G4Triton::TritonDefinition();
335  } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
336  theKindOfFragment = G4He3::He3Definition();
337  } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
338  theKindOfFragment = G4Alpha::AlphaDefinition();;
339  } else {
340 
341  // fragment
342  eexc = (*iList)->GetExcitationEnergy();
343  if(eexc < minExcitation) { eexc = 0.0; }
344  theKindOfFragment = theTableOfIons->GetIon(theFragmentZ,theFragmentA,eexc);
345  /*
346  G4cout << "### Find ion Z= " << theFragmentZ << " A= " << theFragmentA
347  << " Eexc(MeV)= " << eexc/MeV << " " << theKindOfFragment << G4endl;
348  */
349  }
350  // fragment identified
351  if(theKindOfFragment) {
352  G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
353  theNew->SetMomentum((*iList)->GetMomentum().vect());
354  theNew->SetTotalEnergy(etot);
355  theNew->SetFormationTime((*iList)->GetCreationTime());
356  theReactionProductVector->push_back(theNew);
357 
358  // fragment not found out ground state is created
359  } else {
360  theKindOfFragment = theTableOfIons->GetIon(theFragmentZ,theFragmentA,0.0);
361  if(theKindOfFragment) {
362  G4ThreeVector mom(0.0,0.0,0.0);
363  G4double ionmass = theKindOfFragment->GetPDGMass();
364  if(etot <= ionmass) {
365  etot = ionmass;
366  } else {
367  G4double ptot = std::sqrt((etot - ionmass)*(etot + ionmass));
368  mom = ((*iList)->GetMomentum().vect().unit())*ptot;
369  }
370  G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
371  theNew->SetMomentum(mom);
372  theNew->SetTotalEnergy(etot);
373  theNew->SetFormationTime((*iList)->GetCreationTime());
374  theReactionProductVector->push_back(theNew);
375  /*
376  G4cout << "### Find ion Z= " << theFragmentZ << " A= " << theFragmentA
377  << " ground state, energy corrected " << theKindOfFragment << G4endl;
378  } else {
379  G4cout << "### Find ion Z= " << theFragmentZ
380  << " A= " << theFragmentA << " failed " << G4endl;
381  */
382  }
383  }
384  delete (*iList);
385  }
386  return theReactionProductVector;
387 }
388 
390 {
391  //for inverse cross section choice
393  //for the choice of superimposed Coulomb Barrier for inverse cross sections
396 }
397 
399 {
400  if(ptr && ptr != theEvaporation) {
401  delete theEvaporation;
402  theEvaporation = ptr;
404  SetParameters();
405  isEvapLocal = false;
406  }
407 }
408 
409 void
411 {
412  if(ptr && ptr != theMultiFragmentation) {
413  delete theMultiFragmentation;
414  theMultiFragmentation = ptr;
415  }
416 }
417 
419 {
420  if(ptr && ptr != theFermiModel) {
421  delete theFermiModel;
422  theFermiModel = ptr;
423  }
424 }
425 
426 void
428 {
429  if(ptr && ptr != thePhotonEvaporation) {
430  thePhotonEvaporation = ptr;
432  ptr->Initialise();
433  }
434 }
435 
437 {
438  maxZForFermiBreakUp = aZ;
439 }
440 
442 {
443  maxAForFermiBreakUp = anA;
444 }
445 
447 {
450 }
451 
453 {
454  minEForMultiFrag = anE;
455 }
456 void G4ExcitationHandler::ModelDescription(std::ostream& outFile) const
457 {
458  outFile << "G4ExcitationHandler description\n"
459  << "This class samples de-excitation of excited nucleus using\n"
460  << "Fermi Break-up model for light fragments (Z < 9, A < 17), "
461  << "evaporation, fission, and photo-evaporation models. Evaporated\n"
462  << "particle may be proton, neutron, and other light fragment \n"
463  << "(Z < 13, A < 29). During photon evaporation produced gamma \n"
464  << "or electrons due to internal conversion \n";
465 }
466 
467 
468 
469 
470 
471 
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:491
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
virtual void BreakFragment(G4FragmentVector *, G4Fragment *theNucleus)
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:251
void UseSICB(G4bool use)
std::vector< G4Fragment * > theEvapList
bool G4bool
Definition: G4Types.hh:79
virtual void InitialiseChannels()
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:273
G4VMultiFragmentation * theMultiFragmentation
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:256
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:268