| Geant4 Cross Reference |
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"
52 #include "G4EvaporationGEMFactory.hh"
53 #include "G4EvaporationGEMFactoryVI.hh"
54 #include "G4EvaporationDefaultGEMFactory.hh"
55 #include "G4NistManager.hh"
56 #include "G4VFermiBreakUp.hh"
57 #include "G4FermiBreakUpVI.hh"
58 #include "G4PhotonEvaporation.hh"
59 #include "G4VEvaporationChannel.hh"
60 #include "G4ParticleTable.hh"
61 #include "G4IonTable.hh"
62 #include "G4NuclearLevelData.hh"
63 #include "G4LevelManager.hh"
64 #include "G4UnstableFragmentBreakUp.hh"
65 #include "Randomize.hh"
66
67 G4Evaporation::G4Evaporation(G4VEvaporationChannel* photoEvaporation)
68 : fVerbose(0), minExcitation(0.1*CLHEP::keV)
69 {
70 if (nullptr != photoEvaporation) {
71 SetPhotonEvaporation(photoEvaporation);
72 } else {
73 SetPhotonEvaporation(new G4PhotonEvaporation());
74 }
75
76 channelType = fCombined;
77
78 fLevelData = G4NuclearLevelData::GetInstance();
79 theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
80 nist = G4NistManager::Instance();
81 unstableBreakUp = new G4UnstableFragmentBreakUp();
82 }
83
84 G4Evaporation::~G4Evaporation()
85 {
86 delete unstableBreakUp;
87 }
88
89 void G4Evaporation::InitialiseChannels()
90 {
91 if (isInitialised) { return; }
92
93 G4DeexPrecoParameters* param = fLevelData->GetParameters();
94 minExcitation = param->GetMinExcitation();
95 fVerbose = param->GetVerbose();
96 unstableBreakUp->SetVerbose(fVerbose);
97
98 if (nullptr == theChannelFactory) {
99 G4DeexChannelType type = param->GetDeexChannelsType();
100 if(type == fCombined) { SetCombinedChannel(); }
101 else if(type == fGEM) { SetGEMChannel(); }
102 else if(type == fEvaporation) { SetDefaultChannel(); }
103 else if(type == fGEMVI) { SetGEMVIChannel(); }
104 }
105 isInitialised = true;
106 }
107
108 void G4Evaporation::InitialiseChannelFactory()
109 {
110 if (nullptr == theFBU) {
111 theFBU = new G4FermiBreakUpVI();
112 theFBU->Initialise();
113 }
114 theChannels = theChannelFactory->GetChannel();
115 nChannels = theChannels->size();
116 probabilities.resize(nChannels, 0.0);
117
118 if(fVerbose > 1) {
119 G4cout << "### G4Evaporation::InitialiseChannelFactory for "
120 << nChannels << " channels " << this << G4endl;
121 }
122 for (std::size_t i=0; i<nChannels; ++i) {
123 (*theChannels)[i]->SetOPTxs(OPTxs);
124 (*theChannels)[i]->Initialise();
125 }
126 }
127
128 void G4Evaporation::SetDefaultChannel()
129 {
130 if (fEvaporation != channelType || nullptr == theChannelFactory) {
131 channelType = fEvaporation;
132 CleanChannels();
133 delete theChannelFactory;
134 theChannelFactory = new G4EvaporationFactory(thePhotonEvaporation);
135 InitialiseChannelFactory();
136 }
137 }
138
139 void G4Evaporation::SetGEMChannel()
140 {
141 if (fGEM != channelType || nullptr == theChannelFactory) {
142 channelType = fGEM;
143 CleanChannels();
144 delete theChannelFactory;
145 theChannelFactory = new G4EvaporationGEMFactory(thePhotonEvaporation);
146 InitialiseChannelFactory();
147 }
148 }
149
150 void G4Evaporation::SetGEMVIChannel()
151 {
152 if (fGEMVI != channelType || nullptr == theChannelFactory) {
153 channelType = fGEMVI;
154 CleanChannels();
155 delete theChannelFactory;
156 theChannelFactory = new G4EvaporationGEMFactoryVI(thePhotonEvaporation);
157 InitialiseChannelFactory();
158 }
159 }
160
161 void G4Evaporation::SetCombinedChannel()
162 {
163 if (fCombined != channelType || nullptr == theChannelFactory) {
164 channelType = fCombined;
165 CleanChannels();
166 delete theChannelFactory;
167 theChannelFactory =
168 new G4EvaporationDefaultGEMFactory(thePhotonEvaporation);
169 InitialiseChannelFactory();
170 }
171 }
172
173 void G4Evaporation::BreakFragment(G4FragmentVector* theResult,
174 G4Fragment* theResidualNucleus)
175 {
176 if (!isInitialised) { InitialiseChannels(); }
177
178 G4double totprob, prob, oldprob = 0.0;
179 const G4double limFact = 1.e-6;
180 std::size_t maxchannel, i;
181
182 G4int Amax = theResidualNucleus->GetA_asInt();
183 if(fVerbose > 1) {
184 G4cout << "### G4Evaporation::BreakItUp loop" << G4endl;
185 }
186 CLHEP::HepRandomEngine* rndm = G4Random::getTheEngine();
187
188 // Starts loop over evaporated particles, loop is limited by number
189 // of nucleons
190 for(G4int ia=0; ia<Amax; ++ia) {
191
192 // g,n,p and light fragments - evaporation is finished
193 G4int Z = theResidualNucleus->GetZ_asInt();
194 G4int A = theResidualNucleus->GetA_asInt();
195 if(A <= 1) { break; }
196 G4double Eex = theResidualNucleus->GetExcitationEnergy();
197
198 // stop deexcitation loop if residual can be deexcited by FBU
199 if(theFBU->IsApplicable(Z, A, Eex)) { break; }
200
201 // check if it is stable, then finish evaporation
202 G4double abun = nist->GetIsotopeAbundance(Z, A);
203 // stop deecitation loop in the case of a cold stable fragment
204 if(Eex <= minExcitation &&
205 (abun > 0.0 || (A == 3 && (Z == 1 || Z == 2)))) { break; }
206
207 totprob = 0.0;
208 maxchannel = nChannels;
209 if(fVerbose > 1) {
210 G4cout << "Evaporation# " << ia << " Z= " << Z << " A= " << A
211 << " Eex(MeV)= " << theResidualNucleus->GetExcitationEnergy()
212 << " aban= " << abun << G4endl;
213 }
214 // loop over evaporation channels
215 for(i=0; i<nChannels; ++i) {
216 prob = (*theChannels)[i]->GetEmissionProbability(theResidualNucleus);
217 if(fVerbose > 1 && prob > 0.0) {
218 G4cout << " Channel# " << i << " prob= " << prob << G4endl;
219 }
220 totprob += prob;
221 probabilities[i] = totprob;
222
223 // if two recent probabilities are near zero stop computations
224 if (i > 8) {
225 if (prob <= totprob*limFact && oldprob <= totprob*limFact) {
226 maxchannel = i + 1;
227 break;
228 }
229 }
230 oldprob = prob;
231 }
232
233 // photon evaporation in the case of no other channels available
234 // do evaporation chain and return back ground state fragment
235 if(0.0 < totprob && probabilities[0] == totprob) {
236 if(fVerbose > 1) {
237 G4cout << "$$$ Start chain of gamma evaporation" << G4endl;
238 }
239 (*theChannels)[0]->BreakUpChain(theResult, theResidualNucleus);
240
241 // release residual stable fragment
242 if(abun > 0.0) {
243 theResidualNucleus->SetLongLived(true);
244 break;
245 }
246 // release residual fragment known to FBU
247 Eex = theResidualNucleus->GetExcitationEnergy();
248 if(theFBU->IsApplicable(Z, A, Eex)) { break; }
249
250 // release residual fragment with non-zero life time
251 if(theResidualNucleus->IsLongLived()) { break; }
252 totprob = 0.0;
253 }
254
255 if(0.0 == totprob && A < 30) {
256 // if residual fragment is exotic, then it forced to decay
257 // if success, then decay product is added to results
258 if(fVerbose > 1) {
259 G4cout << "$$$ Decay exotic fragment" << G4endl;
260 }
261 if(unstableBreakUp->BreakUpChain(theResult, theResidualNucleus)) {
262 continue;
263 }
264 // release if it is not possible to decay
265 break;
266 }
267
268 // select channel
269 totprob *= rndm->flat();
270
271 // loop over evaporation channels
272 for (i=0; i<maxchannel; ++i) {
273 if (probabilities[i] >= totprob) { break; }
274 }
275
276 if(fVerbose > 1) { G4cout << "$$$ Channel # " << i << G4endl; }
277 G4Fragment* frag = (*theChannels)[i]->EmittedFragment(theResidualNucleus);
278 if(fVerbose > 2 && frag) { G4cout << " " << *frag << G4endl; }
279
280 // normaly a fragment should be created
281 if(nullptr != frag) { theResult->push_back(frag); }
282 else { break; }
283 }
284 }
285