Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 102590 2017-02-09 10:13:30Z 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 "G4NistManager.hh"
57 #include "G4VFermiBreakUp.hh"
58 #include "G4PhotonEvaporation.hh"
59 #include "G4VEvaporationChannel.hh"
60 #include "G4ParticleTable.hh"
61 #include "G4IonTable.hh"
62 #include "G4NuclearLevelData.hh"
64 #include "Randomize.hh"
65 
67  : G4VEvaporation(),fVerbose(0),nChannels(0),minExcitation(0.1*keV),
68  isInitialised(false)
69 {
70  if(photoEvaporation) { SetPhotonEvaporation(photoEvaporation); }
72 
75  InitialiseChannelFactory();
76 
77  theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
78  nist = G4NistManager::Instance();
79  unstableBreakUp = new G4UnstableFragmentBreakUp();
80 }
81 
83 {
84  delete unstableBreakUp;
85 }
86 
88 {
89  if(isInitialised) { return; }
90 
91  G4DeexPrecoParameters* param =
93  minExcitation = param->GetMinExcitation();
94 
95  G4DeexChannelType type = param->GetDeexChannelsType();
96  if(type == fCombined) { SetCombinedChannel(); }
97  else if(type == fGEM) { SetGEMChannel(); }
98 
99  isInitialised = true;
100 }
101 
102 void G4Evaporation::InitialiseChannelFactory()
103 {
105  nChannels = theChannels->size();
106  probabilities.resize(nChannels, 0.0);
107 
108  if(fVerbose > 0) {
109  G4cout << "### G4Evaporation::InitialiseChannelFactory for "
110  << nChannels << " channels " << this << G4endl;
111  }
112  for(size_t i=0; i<nChannels; ++i) {
113  (*theChannels)[i]->SetOPTxs(OPTxs);
114  (*theChannels)[i]->Initialise();
115  }
116 }
117 
119 {
120  if(fEvaporation != channelType) {
122  CleanChannels();
123  delete theChannelFactory;
125  InitialiseChannelFactory();
126  }
127 }
128 
130 {
131  if(fGEM != channelType) {
133  CleanChannels();
134  delete theChannelFactory;
136  InitialiseChannelFactory();
137  }
138 }
139 
141 {
142  if(fCombined != channelType) {
144  CleanChannels();
145  delete theChannelFactory;
147  InitialiseChannelFactory();
148  }
149 }
150 
152  G4Fragment* theResidualNucleus)
153 {
154  if(!isInitialised) { InitialiseChannels(); }
155 
156  G4double totprob, prob, oldprob = 0.0;
157  size_t maxchannel, i;
158 
159  G4int Amax = theResidualNucleus->GetA_asInt();
160 
161  // Starts loop over evaporated particles, loop is limited by number
162  // of nucleons
163  for(G4int ia=0; ia<Amax; ++ia) {
164 
165  // g,n,p and light fragments - evaporation is finished
166  G4int Z = theResidualNucleus->GetZ_asInt();
167  G4int A = theResidualNucleus->GetA_asInt();
168  G4double Eex = theResidualNucleus->GetExcitationEnergy();
169 
170  // stop deecitation loop if residual can be deexcited by FBU
171  if(theFBU->IsApplicable(Z, A, Eex)) { break; }
172 
173  // check if it is stable, then finish evaporation
174  G4double abun = nist->GetIsotopeAbundance(Z, A);
175  if(fVerbose > 0) {
176  G4cout << "### G4Evaporation::BreakItUp step " << ia << " Z= " << Z
177  << " A= " << A << " Eex(MeV)= "
178  << theResidualNucleus->GetExcitationEnergy()
179  << " aban= " << abun << G4endl;
180  }
181  // stop deecitation loop in the case of a cold stable fragment
182  if(Eex <= minExcitation && abun > 0.0) { break; }
183 
184  totprob = 0.0;
185  maxchannel = nChannels;
186  if(fVerbose > 1) {
187  G4cout << "### Evaporation loop #" << ia
188  << " Fragment: " << theResidualNucleus << G4endl;
189  }
190  // loop over evaporation channels
191  for(i=0; i<nChannels; ++i) {
192  prob = (*theChannels)[i]->GetEmissionProbability(theResidualNucleus);
193  if(fVerbose > 0 && prob > 0.0) {
194  G4cout << " Channel# " << i << " prob= " << prob << G4endl;
195  }
196  totprob += prob;
197  probabilities[i] = totprob;
198 
199  // if two recent probabilities are near zero stop computations
200  if(i>=8 && prob > 0.0) {
201  if(prob <= totprob*1.e-8 && oldprob <= totprob*1.e-8) {
202  maxchannel = i+1;
203  break;
204  }
205  }
206  oldprob = prob;
207  }
208 
209  // photon evaporation in the case of no other channels available
210  // do evaporation chain and reset total probability
211  if(0.0 < totprob && probabilities[0] == totprob) {
212  if(fVerbose > 1) {
213  G4cout << "Start chain of gamma evaporation" << G4endl;
214  }
215  (*theChannels)[0]->BreakUpChain(theResult, theResidualNucleus);
216  totprob = 0.0;
217  }
218 
219  // stable fragment - evaporation is finished
220  if(0.0 == totprob) {
221 
222  // if fragment is exotic, then force its decay
223  if(0.0 == abun) {
224  if(!unstableBreakUp->BreakUpChain(theResult, theResidualNucleus))
225  { break; }
226  if(fVerbose > 1) { G4cout << "$$$ Decay exotic fragment" << G4endl; }
227  } else {
228  break;
229  }
230  }
231 
232  // select channel
233  totprob *= G4UniformRand();
234  // loop over evaporation channels
235  for(i=0; i<maxchannel; ++i) { if(probabilities[i] >= totprob) { break; } }
236 
237  if(fVerbose > 1) { G4cout << "Channel # " << i << G4endl; }
238  G4Fragment* frag = (*theChannels)[i]->EmittedFragment(theResidualNucleus);
239  if(fVerbose > 1 && frag) { G4cout << " " << *frag << G4endl; }
240 
241  // normaly a fragment should be created
242  if(frag) { theResult->push_back(frag); }
243  else { break; }
244  }
245  // loop is stopped, residual fragment is added to the results
246  theResult->push_back(theResidualNucleus);
247 }
void SetDefaultChannel()
virtual std::vector< G4VEvaporationChannel * > * GetChannel()=0
G4DeexChannelType GetDeexChannelsType() const
G4VFermiBreakUp * theFBU
virtual G4bool IsApplicable(G4int Z, G4int A, G4double mass) const =0
virtual G4bool BreakUpChain(G4FragmentVector *, G4Fragment *) final
G4VEvaporationChannel * thePhotonEvaporation
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4Evaporation(G4VEvaporationChannel *photoEvaporation=nullptr)
std::vector< G4VEvaporationChannel * > * theChannels
void SetGEMChannel()
G4IonTable * GetIonTable() const
virtual void InitialiseChannels() final
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:266
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
G4DeexPrecoParameters * GetParameters()
G4double GetIsotopeAbundance(G4int Z, G4int N) const
virtual void BreakFragment(G4FragmentVector *, G4Fragment *theNucleus) final
void SetCombinedChannel()
static G4ParticleTable * GetParticleTable()
virtual void SetPhotonEvaporation(G4VEvaporationChannel *ptr)
virtual ~G4Evaporation()
G4int GetZ_asInt() const
Definition: G4Fragment.hh:271
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4VEvaporationFactory * theChannelFactory
static constexpr double keV
Definition: G4SIunits.hh:216
static G4NuclearLevelData * GetInstance()
G4double GetMinExcitation() const
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:283