Geant4  10.00.p01
G4Evaporation.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 //
27 // $Id: G4Evaporation.cc 71639 2013-06-19 15:07:54Z gcosmo $
28 //
29 // Hadronic Process: Nuclear De-excitations
30 // by V. Lara (Oct 1998)
31 //
32 // Alex Howard - added protection for negative probabilities in the sum, 14/2/07
33 //
34 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
35 // cross section option
36 // JMQ (06 September 2008) Also external choices have been added for
37 // superimposed Coulomb barrier (if useSICBis set true, by default is false)
38 //
39 // V.Ivanchenko (27 July 2009) added G4EvaporationDefaultGEMFactory option
40 // V.Ivanchenko (10 May 2010) rewrited BreakItUp method: do not make new/delete
41 // photon channel first, fission second,
42 // added G4UnstableFragmentBreakUp to decay
43 // unphysical fragments (like 2n or 2p),
44 // use Z and A integer
45 // V.Ivanchenko (22 April 2011) added check if a fragment can be deexcited
46 // by the FermiBreakUp model
47 // V.Ivanchenko (23 January 2012) added pointer of G4VPhotonEvaporation
48 // V.Ivanchenko (6 May 2013) added check of existence of residual ion
49 // in the ion table
50 
51 #include "G4Evaporation.hh"
52 #include "G4SystemOfUnits.hh"
53 #include "G4EvaporationFactory.hh"
56 #include "G4HadronicException.hh"
57 #include "G4NistManager.hh"
58 #include "G4FermiFragmentsPool.hh"
59 #include "G4PhotonEvaporation.hh"
60 
61 #include "G4ParticleTable.hh"
62 #include "G4IonTable.hh"
63 
65  : theChannels(0),nChannels(0)
66 {
69  SetParameters();
72 }
73 
75  : theChannels(0),nChannels(0)
76 {
77  if(photoEvaporation) { SetPhotonEvaporation(photoEvaporation); }
79 
81  SetParameters();
84 }
85 
87 {
88  CleanChannels();
89  delete thePhotonEvaporation;
90  delete theChannelFactory;
91 }
92 
94 {
95  for (size_t i=1; i<nChannels; ++i) { delete (*theChannels)[i]; }
96  delete theChannels;
97 }
98 
100 {
105  probabilities.reserve(68);
106 }
107 
109 {
110  CleanChannels();
112  nChannels = theChannels->size();
113  probabilities.resize(nChannels, 0.0);
114  Initialise();
115 }
116 
118 {
119  for(size_t i=0; i<nChannels; ++i) {
120  (*theChannels)[i]->SetOPTxs(OPTxs);
121  (*theChannels)[i]->UseSICB(useSICB);
122  }
123 }
124 
126 {
127  delete theChannelFactory;
130 }
131 
133 {
134  delete theChannelFactory;
137 }
138 
140 {
141  delete theChannelFactory;
144 }
145 
147 {
148  if(ptr) { G4VEvaporation::SetPhotonEvaporation(ptr); }
149  if(0 < nChannels) { (*theChannels)[0] = ptr; }
150 }
151 
153 {
154  G4FragmentVector * theResult = new G4FragmentVector;
155  G4FragmentVector * theTempResult;
156  static const G4double Elimit = 3*MeV;
157 
158  // The residual nucleus (after evaporation of each fragment)
159  G4Fragment* theResidualNucleus = new G4Fragment(theNucleus);
160 
161  G4double totprob, prob, oldprob = 0.0;
162  size_t maxchannel, i;
163 
164  G4int Amax = theResidualNucleus->GetA_asInt();
165 
166  // Starts loop over evaporated particles, loop is limited by number
167  // of nucleons
168  for(G4int ia=0; ia<Amax; ++ia) {
169 
170  // g,n,p and light fragments - evaporation is finished
171  G4int Z = theResidualNucleus->GetZ_asInt();
172  G4int A = theResidualNucleus->GetA_asInt();
173 
174  // stop deecitation loop if can be deexcited by FBU
175  if(maxZforFBU > Z && maxAforFBU >= A) {
176  theResult->push_back(theResidualNucleus);
177  return theResult;
178  }
179 
180  // check if it is stable, then finish evaporation
181  G4double abun = nist->GetIsotopeAbundance(Z, A);
182  /*
183  G4cout << "### G4Evaporation::BreakItUp step " << ia << " Z= " << Z
184  << " A= " << A << " Eex(MeV)= "
185  << theResidualNucleus->GetExcitationEnergy()
186  << " aban= " << abun << G4endl;
187  */
188  // stop deecitation loop in the case of the cold stable fragment
189  G4double Eex = theResidualNucleus->GetExcitationEnergy();
190  if(Eex <= minExcitation && abun > 0.0) {
191  theResult->push_back(theResidualNucleus);
192  return theResult;
193  }
194 
195  totprob = 0.0;
196  maxchannel = nChannels;
197  /*
198  G4cout << "### Evaporation loop #" << ia
199  << " Fragment: " << theResidualNucleus << G4endl;
200  */
201  // loop over evaporation channels
202  for(i=0; i<nChannels; ++i) {
203  prob = (*theChannels)[i]->GetEmissionProbability(theResidualNucleus);
204  //G4cout << " Channel# " << i << " prob= " << prob << G4endl;
205 
206  totprob += prob;
207  probabilities[i] = totprob;
208  // if two recent probabilities are near zero stop computations
209  if(i>=8) {
210  if(prob <= totprob*1.e-8 && oldprob <= totprob*1.e-8) {
211  maxchannel = i+1;
212  break;
213  }
214  }
215  oldprob = prob;
216  // protection for very excited fragment - avoid GEM
217  if(7 == i && Eex > Elimit*A) {
218  maxchannel = 8;
219  break;
220  }
221  }
222 
223  // photon evaporation in the case of no other channels available
224  // do evaporation chain and reset total probability
225  if(0.0 < totprob && probabilities[0] == totprob) {
226  //G4cout << "Start gamma evaporation" << G4endl;
227  theTempResult = (*theChannels)[0]->BreakUpFragment(theResidualNucleus);
228  if(theTempResult) {
229  size_t nsec = theTempResult->size();
230  for(size_t j=0; j<nsec; ++j) {
231  theResult->push_back((*theTempResult)[j]);
232  }
233  delete theTempResult;
234  }
235  totprob = 0.0;
236  }
237 
238  // stable fragnent - evaporation is finished
239  if(0.0 == totprob) {
240 
241  // if fragment is exotic, then force its decay
242  if(0.0 == abun && Z < 20) {
243  //G4cout << "$$$ Decay exotic fragment" << G4endl;
244  theTempResult = unstableBreakUp.BreakUpFragment(theResidualNucleus);
245  if(theTempResult) {
246  size_t nsec = theTempResult->size();
247  for(size_t j=0; j<nsec; ++j) {
248  theResult->push_back((*theTempResult)[j]);
249  }
250  delete theTempResult;
251  }
252  }
253 
254  // save residual fragment
255  theResult->push_back(theResidualNucleus);
256  return theResult;
257  }
258 
259  // select channel
260  totprob *= G4UniformRand();
261  // loop over evaporation channels
262  for(i=0; i<maxchannel; ++i) { if(probabilities[i] >= totprob) { break; } }
263 
264  // this should not happen
265  if(i >= nChannels) { i = nChannels - 1; }
266 
267 
268  // single photon evaporation, primary pointer is kept
269  if(0 == i) {
270  //G4cout << "Single gamma" << G4endl;
271  G4Fragment* gamma = (*theChannels)[0]->EmittedFragment(theResidualNucleus);
272  if(gamma) { theResult->push_back(gamma); }
273 
274  // fission, return results to the main loop if fission is succesful
275  } else if(1 == i) {
276  //G4cout << "Fission" << G4endl;
277  theTempResult = (*theChannels)[1]->BreakUp(*theResidualNucleus);
278  if(theTempResult) {
279  size_t nsec = theTempResult->size();
280  G4bool deletePrimary = true;
281  for(size_t j=0; j<nsec; ++j) {
282  if(theResidualNucleus == (*theTempResult)[j]) { deletePrimary = false; }
283  theResult->push_back((*theTempResult)[j]);
284  }
285  if(deletePrimary) { delete theResidualNucleus; }
286  delete theTempResult;
287  return theResult;
288  }
289 
290  // other channels
291  } else {
292  //G4cout << "Channel # " << i << G4endl;
293  theTempResult = (*theChannels)[i]->BreakUp(*theResidualNucleus);
294  if(theTempResult) {
295  size_t nsec = theTempResult->size();
296  if(nsec > 0) {
297  --nsec;
298  for(size_t j=0; j<nsec; ++j) {
299  theResult->push_back((*theTempResult)[j]);
300  }
301  // if the residual change its pointer
302  // then delete previous residual fragment and update to the new
303  if(theResidualNucleus != (*theTempResult)[nsec] ) {
304  delete theResidualNucleus;
305  theResidualNucleus = (*theTempResult)[nsec];
306  }
307  }
308  delete theTempResult;
309  }
310  }
311  }
312 
313  // loop is stopped, save residual
314  theResult->push_back(theResidualNucleus);
315 
316  return theResult;
317 }
318 
G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)
void SetDefaultChannel()
virtual std::vector< G4VEvaporationChannel * > * GetChannel()=0
static const double MeV
Definition: G4SIunits.hh:193
virtual void SetPhotonEvaporation(G4VEvaporationChannel *ptr)
virtual void Initialise()
std::vector< G4double > probabilities
void CleanChannels()
G4VEvaporationChannel * thePhotonEvaporation
std::vector< G4VEvaporationChannel * > * theChannels
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
void SetGEMChannel()
virtual G4FragmentVector * BreakUpFragment(G4Fragment *fragment)
G4IonTable * GetIonTable() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4int GetA_asInt() const
Definition: G4Fragment.hh:238
bool G4bool
Definition: G4Types.hh:79
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:65
G4IonTable * theTableOfIons
G4double GetIsotopeAbundance(G4int Z, G4int N) const
G4NistManager * nist
static const G4double A[nN]
void SetCombinedChannel()
G4VEvaporationFactory * theChannelFactory
static G4ParticleTable * GetParticleTable()
G4double minExcitation
void InitialiseEvaporation()
virtual void SetPhotonEvaporation(G4VEvaporationChannel *ptr)
G4UnstableFragmentBreakUp unstableBreakUp
virtual ~G4Evaporation()
G4int GetZ_asInt() const
Definition: G4Fragment.hh:243
void SetParameters()
static const double keV
Definition: G4SIunits.hh:195
double G4double
Definition: G4Types.hh:76
static G4FermiFragmentsPool * Instance()
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:255