ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Evaporation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Evaporation.cc
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 // Hadronic Process: Nuclear De-excitations
28 // by V. Lara (Oct 1998)
29 //
30 // Alex Howard - added protection for negative probabilities in the sum, 14/2/07
31 //
32 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
33 // cross section option
34 // JMQ (06 September 2008) Also external choices have been added for
35 // superimposed Coulomb barrier (if useSICBis set true, by default is false)
36 //
37 // V.Ivanchenko (27 July 2009) added G4EvaporationDefaultGEMFactory option
38 // V.Ivanchenko (10 May 2010) rewrited BreakItUp method: do not make new/delete
39 // photon channel first, fission second,
40 // added G4UnstableFragmentBreakUp to decay
41 // unphysical fragments (like 2n or 2p),
42 // use Z and A integer
43 // V.Ivanchenko (22 April 2011) added check if a fragment can be deexcited
44 // by the FermiBreakUp model
45 // V.Ivanchenko (23 January 2012) added pointer of G4VPhotonEvaporation
46 // V.Ivanchenko (6 May 2013) added check of existence of residual ion
47 // in the ion table
48 
49 #include "G4Evaporation.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4EvaporationFactory.hh"
55 #include "G4NistManager.hh"
56 #include "G4VFermiBreakUp.hh"
57 #include "G4PhotonEvaporation.hh"
58 #include "G4VEvaporationChannel.hh"
59 #include "G4ParticleTable.hh"
60 #include "G4IonTable.hh"
61 #include "G4NuclearLevelData.hh"
62 #include "G4LevelManager.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 
74  theChannelFactory = nullptr;
75 
80  /*
81  G4cout << "G4Evaporation() " << this << " thePhotonEvaporation: "
82  << photoEvaporation << " UnstableFragmentBreakUp: "
83  << unstableBreakUp << G4endl;
84  */
85 }
86 
88 {
89  delete unstableBreakUp;
90 }
91 
93 {
94  if(isInitialised) { return; }
95 
98  fVerbose = param->GetVerbose();
100 
101  G4DeexChannelType type = param->GetDeexChannelsType();
102  if(type == fCombined) { SetCombinedChannel(); }
103  else if(type == fGEM) { SetGEMChannel(); }
104  else if(type == fEvaporation) { SetDefaultChannel(); }
105  else if(type == fGEMVI) { SetGEMVIChannel(); }
106 
107  isInitialised = true;
108 }
109 
111 {
113  nChannels = theChannels->size();
114  probabilities.resize(nChannels, 0.0);
115 
116  if(fVerbose > 1) {
117  G4cout << "### G4Evaporation::InitialiseChannelFactory for "
118  << nChannels << " channels " << this << G4endl;
119  }
120  for(size_t i=0; i<nChannels; ++i) {
121  (*theChannels)[i]->SetOPTxs(OPTxs);
122  (*theChannels)[i]->Initialise();
123  }
124 }
125 
127 {
128  if(fEvaporation != channelType) {
130  if(theChannelFactory) {
131  CleanChannels();
132  delete theChannelFactory;
133  }
136  }
137 }
138 
140 {
141  if(fGEM != channelType) {
143  if(theChannelFactory) {
144  CleanChannels();
145  delete theChannelFactory;
146  }
149  }
150 }
151 
153 {
154  if(fGEMVI != channelType) {
156  if(theChannelFactory) {
157  CleanChannels();
158  delete theChannelFactory;
159  }
162  }
163 }
164 
166 {
167  if(fCombined != channelType) {
169  if(theChannelFactory) {
170  CleanChannels();
171  delete theChannelFactory;
172  }
176  }
177 }
178 
180  G4Fragment* theResidualNucleus)
181 {
183 
184  G4double totprob, prob, oldprob = 0.0;
185  size_t maxchannel, i;
186 
187  G4int Amax = theResidualNucleus->GetA_asInt();
188  if(fVerbose > 1) {
189  G4cout << "### G4Evaporation::BreakItUp loop" << G4endl;
190  }
191 
192  // Starts loop over evaporated particles, loop is limited by number
193  // of nucleons
194  for(G4int ia=0; ia<Amax; ++ia) {
195 
196  // g,n,p and light fragments - evaporation is finished
197  G4int Z = theResidualNucleus->GetZ_asInt();
198  G4int A = theResidualNucleus->GetA_asInt();
199  if(A <= 1) { break; }
200  G4double Eex = theResidualNucleus->GetExcitationEnergy();
201 
202  // stop deecitation loop if residual can be deexcited by FBU
203  if(theFBU->IsApplicable(Z, A, Eex)) { break; }
204 
205  // check if it is stable, then finish evaporation
206  G4double abun = nist->GetIsotopeAbundance(Z, A);
207  // stop deecitation loop in the case of a cold stable fragment
208  if(Eex <= minExcitation &&
209  (abun > 0.0 || (A == 3 && (Z == 1 || Z == 2)))) { break; }
210 
211  totprob = 0.0;
212  maxchannel = nChannels;
213  if(fVerbose > 1) {
214  G4cout << "Evaporation# " << ia << " Z= " << Z << " A= " << A
215  << " Eex(MeV)= " << theResidualNucleus->GetExcitationEnergy()
216  << " aban= " << abun << G4endl;
217  }
218  // loop over evaporation channels
219  for(i=0; i<nChannels; ++i) {
220  prob = (*theChannels)[i]->GetEmissionProbability(theResidualNucleus);
221  if(fVerbose > 2 && prob > 0.0) {
222  G4cout << " Channel# " << i << " prob= " << prob << G4endl;
223  }
224  totprob += prob;
225  probabilities[i] = totprob;
226 
227  // if two recent probabilities are near zero stop computations
228  if(i>=8 && prob > 0.0) {
229  if(prob <= totprob*1.e-8 && oldprob <= totprob*1.e-8) {
230  maxchannel = i+1;
231  break;
232  }
233  }
234  oldprob = prob;
235  }
236 
237  // photon evaporation in the case of no other channels available
238  // do evaporation chain and return back ground state fragment
239  if(0.0 < totprob && probabilities[0] == totprob) {
240  if(fVerbose > 1) {
241  G4cout << "$$$ Start chain of gamma evaporation" << G4endl;
242  }
243  (*theChannels)[0]->BreakUpChain(theResult, theResidualNucleus);
244  totprob = 0.0;
245  }
246 
247  // stable fragment - evaporation is finished
248  if(0.0 == totprob) {
249 
250  // release fragment known to DB
251  if(fLevelData->GetLevelManager(Z, A)) { break; }
252 
253  // if fragment is exotic, then it forced to decay
254  // if success, then decay product is added to results
255  if(fVerbose > 1) {
256  G4cout << "$$$ Decay exotic fragment" << G4endl;
257  }
258  if(unstableBreakUp->BreakUpChain(theResult, theResidualNucleus)) {
259  continue;
260  }
261  // release if it is not possible to decay
262  break;
263  }
264 
265  // select channel
266  totprob *= G4UniformRand();
267  // loop over evaporation channels
268  for(i=0; i<maxchannel; ++i) { if(probabilities[i] >= totprob) { break; } }
269 
270  if(fVerbose > 1) { G4cout << "$$$ Channel # " << i << G4endl; }
271  G4Fragment* frag = (*theChannels)[i]->EmittedFragment(theResidualNucleus);
272  if(fVerbose > 2 && frag) { G4cout << " " << *frag << G4endl; }
273 
274  // normaly a fragment should be created
275  if(frag) { theResult->push_back(frag); }
276  else { break; }
277  }
278 }