Geant4  10.00.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 79200 2014-02-20 11:22:39Z 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 {
94 
100  SetParameters();
102 }
103 
105 {
106  if(isEvapLocal) { delete theEvaporation; }
107  delete theMultiFragmentation;
108  delete theFermiModel;
109 }
110 
112 {
113  //for inverse cross section choice
115  //for the choice of superimposed Coulomb Barrier for inverse cross sections
118 }
119 
121 G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const
122 {
123  //G4cout << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@" << G4endl;
124 
125  // Variables existing until end of method
126  G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
127 
128  // pointer to fragment vector which receives temporal results
129  G4FragmentVector * theTempResult = 0;
130 
131  // list of fragments to apply Evaporation or Fermi Break-Up
132  std::list<G4Fragment*> theEvapList;
133 
134  // list of fragments to apply PhotonEvaporation
135  std::list<G4Fragment*> thePhotoEvapList;
136 
137  // list of fragments to store final result
138  std::list<G4Fragment*> theResults;
139  //
140  //G4cout << theInitialState << G4endl;
141 
142  // Variables to describe the excited configuration
143  G4double exEnergy = theInitialState.GetExcitationEnergy();
144  G4int A = theInitialState.GetA_asInt();
145  G4int Z = theInitialState.GetZ_asInt();
146 
148 
149  // In case A <= 1 the fragment will not perform any nucleon emission
150  if (A <= 1)
151  {
152  theResults.push_back( theInitialStatePtr );
153  }
154  // check if a fragment is stable
155  else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0)
156  {
157  theResults.push_back( theInitialStatePtr );
158  }
159  else
160  {
161  // JMQ 150909: first step in de-excitation is treated separately
162  // Fragments after the first step are stored in theEvapList
163  // Statistical Multifragmentation will take place only once
164  //
165  // move to evaporation loop
167  || exEnergy <= minEForMultiFrag*A)
168  {
169  theEvapList.push_back(theInitialStatePtr);
170  }
171  else
172  {
173  theTempResult = theMultiFragmentation->BreakItUp(theInitialState);
174  if(!theTempResult) { theEvapList.push_back(theInitialStatePtr); }
175  else {
176  size_t nsec = theTempResult->size();
177  if(0 == nsec) { theEvapList.push_back(theInitialStatePtr); }
178  else {
179  // secondary are produced
180  // Sort out secondary fragments
181  G4bool deletePrimary = true;
182  G4FragmentVector::iterator j;
183  for (j = theTempResult->begin(); j != theTempResult->end(); ++j) {
184  if((*j) == theInitialStatePtr) { deletePrimary = false; }
185  A = (*j)->GetA_asInt();
186 
187  // gamma, p, n
188  if(A <= 1) { theResults.push_back(*j); }
189 
190  // Analyse fragment A > 1
191  else {
192  G4double exEnergy1 = (*j)->GetExcitationEnergy();
193 
194  // cold fragments
195  if(exEnergy1 < minExcitation) {
196  Z = (*j)->GetZ_asInt();
197  if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
198  theResults.push_back(*j); // stable fragment
199 
200  } else {
201 
202  // check if the cold fragment is from FBU pool
203  const G4VFermiFragment* ffrag = thePool->GetFragment(Z, A);
204  if(ffrag) {
205  if(ffrag->IsStable()) { theResults.push_back(*j); }
206  else { theEvapList.push_back(*j); }
207 
208  // cold fragment may be unstable
209  } else {
210  theEvapList.push_back(*j);
211  }
212  }
213 
214  // hot fragments are unstable
215  } else { theEvapList.push_back(*j); }
216  }
217  }
218  if( deletePrimary ) { delete theInitialStatePtr; }
219  }
220  delete theTempResult;
221  }
222  }
223  }
224  /*
225  G4cout << "## After first step " << theEvapList.size() << " for evap; "
226  << thePhotoEvapList.size() << " for photo-evap; "
227  << theResults.size() << " results. " << G4endl;
228  */
229  // -----------------------------------
230  // FermiBreakUp and De-excitation loop
231  // -----------------------------------
232 
233  std::list<G4Fragment*>::iterator iList;
234  for (iList = theEvapList.begin(); iList != theEvapList.end(); ++iList)
235  {
236  //G4cout << "Next evaporate: " << G4endl;
237  //G4cout << *iList << G4endl;
238  A = (*iList)->GetA_asInt();
239  Z = (*iList)->GetZ_asInt();
240 
241  // Fermi Break-Up
242  G4bool wasFBU = false;
244  {
245  theTempResult = theFermiModel->BreakItUp(*(*iList));
246  wasFBU = true;
247  // if initial fragment returned unchanged try to evaporate it
248  if(1 == theTempResult->size()) {
249  delete *(theTempResult->begin());
250  delete theTempResult;
251  theTempResult = theEvaporation->BreakItUp(*(*iList));
252  }
253  }
254  else // apply Evaporation in another case
255  {
256  theTempResult = theEvaporation->BreakItUp(*(*iList));
257  }
258 
259  G4bool deletePrimary = true;
260  size_t nsec = theTempResult->size();
261  //G4cout << "Nproducts= " << nsec << G4endl;
262 
263  // Sort out secondary fragments
264  if ( nsec > 0 ) {
265  G4FragmentVector::iterator j;
266  for (j = theTempResult->begin(); j != theTempResult->end(); ++j) {
267  if((*j) == (*iList)) { deletePrimary = false; }
268 
269  //G4cout << *j << G4endl;
270  A = (*j)->GetA_asInt();
271  exEnergy = (*j)->GetExcitationEnergy();
272 
273  if(A <= 1) { theResults.push_back(*j); } // gamma, p, n
274 
275  // evaporation is not possible
276  else if(1 == nsec) {
277  if(exEnergy < minExcitation) { theResults.push_back(*j); }
278  else { thePhotoEvapList.push_back(*j); }
279 
280  } else { // Analyse fragment
281 
282  // cold fragment
283  if(exEnergy < minExcitation) {
284  Z = (*j)->GetZ_asInt();
285 
286  // natural isotope
287  if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
288  theResults.push_back(*j); // stable fragment
289 
290  } else {
291  const G4VFermiFragment* ffrag = thePool->GetFragment(Z, A);
292 
293  // isotope from FBU pool
294  if(ffrag) {
295  if(ffrag->IsStable()) { theResults.push_back(*j); }
296  else { theEvapList.push_back(*j); }
297 
298  // isotope may be unstable
299  } else {
300  theEvapList.push_back(*j);
301  }
302  }
303 
304  // hot fragment
305  } else if (wasFBU) {
306  thePhotoEvapList.push_back(*j); // FBU applied only once
307  } else {
308  theEvapList.push_back(*j);
309  }
310  }
311  }
312  }
313  if( deletePrimary ) { delete (*iList); }
314  delete theTempResult;
315  } // end of the loop over theEvapList
316 
317  //G4cout << "## After 2nd step " << theEvapList.size() << " was evap; "
318  // << thePhotoEvapList.size() << " for photo-evap; "
319  // << theResults.size() << " results. " << G4endl;
320 
321  // -----------------------
322  // Photon-Evaporation loop
323  // -----------------------
324 
325  // at this point only photon evaporation is possible
326  for(iList = thePhotoEvapList.begin(); iList != thePhotoEvapList.end(); ++iList)
327  {
328  //G4cout << "Next photon evaporate: " << thePhotonEvaporation << G4endl;
329  //G4cout << *iList << G4endl;
330  exEnergy = (*iList)->GetExcitationEnergy();
331 
332  // only hot fragments
333  if(exEnergy > minExcitation) {
334  theTempResult = thePhotonEvaporation->BreakUpFragment(*iList);
335  size_t nsec = theTempResult->size();
336  //G4cout << "Nproducts= " << nsec << G4endl;
337 
338  // if there is a gamma emission then
339  if (nsec > 0)
340  {
341  G4FragmentVector::iterator j;
342  for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
343  {
344  theResults.push_back(*j);
345  }
346  }
347  delete theTempResult;
348  }
349 
350  // priamry fragment is kept
351  theResults.push_back(*iList);
352 
353  } // end of photon-evaporation loop
354 
355  //G4cout << "## After 3d step " << theEvapList.size() << " was evap; "
356  // << thePhotoEvapList.size() << " was photo-evap; "
357  // << theResults.size() << " results. " << G4endl;
358 
359  G4ReactionProductVector * theReactionProductVector =
361 
362  // MAC (24/07/08)
363  // To optimise the storing speed, we reserve space in memory for the vector
364  theReactionProductVector->reserve( theResults.size() );
365 
366  G4int theFragmentA, theFragmentZ;
367 
368  std::list<G4Fragment*>::iterator i;
369  for (i = theResults.begin(); i != theResults.end(); ++i)
370  {
371  theFragmentA = (*i)->GetA_asInt();
372  theFragmentZ = (*i)->GetZ_asInt();
373  G4double etot= (*i)->GetMomentum().e();
374  G4ParticleDefinition* theKindOfFragment = 0;
375  if (theFragmentA == 0) { // photon or e-
376  theKindOfFragment = (*i)->GetParticleDefinition();
377  } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
378  theKindOfFragment = G4Neutron::NeutronDefinition();
379  } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
380  theKindOfFragment = G4Proton::ProtonDefinition();
381  } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
382  theKindOfFragment = G4Deuteron::DeuteronDefinition();
383  } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
384  theKindOfFragment = G4Triton::TritonDefinition();
385  } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
386  theKindOfFragment = G4He3::He3Definition();
387  } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
388  theKindOfFragment = G4Alpha::AlphaDefinition();;
389  } else {
390 
391  // ground state by default
392  G4double eexc = (*i)->GetExcitationEnergy();
393  G4double excitation = eexc;
394 
395  G4int level = 0;
396  theKindOfFragment =
397  theTableOfIons->GetIon(theFragmentZ,theFragmentA,level);
398  /*
399  G4cout << "### Find ion Z= " << theFragmentZ << " A= " << theFragmentA
400  << " Eexc(MeV)= " << excitation/MeV << " "
401  << theKindOfFragment << G4endl;
402  */
403  // production of an isomer
404  if(eexc > minExcitation) {
405  G4double elevel1 = 0.0;
406  G4double elevel2 = 0.0;
408  for(level=1; level<9; ++level) {
409  ion = theTableOfIons->GetIon(theFragmentZ,theFragmentA,level);
410  //G4cout << level << " " << ion << G4endl;
411  if(ion) {
412  G4Ions* ip = dynamic_cast<G4Ions*>(ion);
413  if(ip) {
414  elevel2 = ip->GetExcitationEnergy();
415  //G4cout<<" Level "<<level<<" E(MeV)= "<<elevel2/MeV<<G4endl;
416  // close level
417  if(std::fabs(eexc - elevel2) < minExcitation) {
418  excitation = eexc - elevel2;
419  theKindOfFragment = ion;
420  break;
421  // previous level was closer
422  } else if(elevel2 - eexc >= eexc - elevel1) {
423  excitation = eexc - elevel1;
424  break;
425  // will check next level and save current
426  } else {
427  theKindOfFragment = ion;
428  excitation = eexc - elevel2;
429  elevel1 = elevel2;
430  }
431  }
432  } else {
433  break;
434  }
435  }
436  }
437  // correction of total energy for ground state isotopes
438  etot += excitation;
439  G4double ionmass = theKindOfFragment->GetPDGMass();
440  if(etot < ionmass) { etot = ionmass; }
441  }
442  if (theKindOfFragment != 0)
443  {
444  G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
445  theNew->SetMomentum((*i)->GetMomentum().vect());
446  theNew->SetTotalEnergy(etot);
447  theNew->SetFormationTime((*i)->GetCreationTime());
448  theReactionProductVector->push_back(theNew);
449  }
450  delete (*i);
451  }
452 
453  return theReactionProductVector;
454 }
455 
457 {
458  if(ptr && ptr != theEvaporation) {
459  delete theEvaporation;
460  theEvaporation = ptr;
462  SetParameters();
463  isEvapLocal = false;
464  }
465 }
466 
467 void
469 {
470  if(ptr && ptr != theMultiFragmentation) {
471  delete theMultiFragmentation;
472  theMultiFragmentation = ptr;
473  }
474 }
475 
477 {
478  if(ptr && ptr != theFermiModel) {
479  delete theFermiModel;
480  theFermiModel = ptr;
481  }
482 }
483 
484 void
486 {
487  if(ptr && ptr != thePhotonEvaporation) {
488  thePhotonEvaporation = ptr;
490  }
491 }
492 
494 {
495  maxZForFermiBreakUp = aZ;
496 }
497 
499 {
500  maxAForFermiBreakUp = anA;
501 }
502 
504 {
507 }
508 
510 {
511  minEForMultiFrag = anE;
512 }
513 
514 
515 
516 
517 
518 
void SetOPTxs(G4int opt)
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
static G4Triton * TritonDefinition()
Definition: G4Triton.cc:90
static G4He3 * He3Definition()
Definition: G4He3.cc:89
virtual G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)=0
G4double GetExcitationEnergy() const
Definition: G4Ions.hh:113
virtual G4FragmentVector * BreakUpFragment(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:449
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
void SetMinEForMultiFrag(G4double anE)
virtual G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)=0
virtual G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)=0
virtual void Initialise()
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState) const
G4bool IsStable() const
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4IonTable * GetIonTable() const
G4int GetA_asInt() const
Definition: G4Fragment.hh:238
Definition: G4Ions.hh:51
void UseSICB(G4bool use)
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
G4VMultiFragmentation * theMultiFragmentation
static const G4double A[nN]
G4VFermiBreakUp * theFermiModel
void SetMaxZForFermiBreakUp(G4int aZ)
const G4VFermiFragment * GetFragment(G4int Z, G4int A)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
virtual void SetPhotonEvaporation(G4VEvaporationChannel *ptr)
void SetMaxAandZForFermiBreakUp(G4int anA, G4int aZ)
G4int GetZ_asInt() const
Definition: G4Fragment.hh:243
void SetEvaporation(G4VEvaporation *ptr)
G4double GetExcitationEnergy(void) const
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:255