| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 //
>> 27 // $Id: G4VEmissionProbability.cc 66241 2012-12-13 18:34:42Z gunter $
>> 28 //
26 // Hadronic Process: Nuclear De-excitations 29 // Hadronic Process: Nuclear De-excitations
27 // by V. Lara (Oct 1998) 30 // by V. Lara (Oct 1998)
28 // 31 //
29 // Modifications: 32 // Modifications:
30 // 28.10.2010 V.Ivanchenko defined members in 33 // 28.10.2010 V.Ivanchenko defined members in constructor and cleaned up
31 34
32 #include "G4VEmissionProbability.hh" 35 #include "G4VEmissionProbability.hh"
33 #include "G4NuclearLevelData.hh" 36 #include "G4NuclearLevelData.hh"
34 #include "G4LevelManager.hh" <<
35 #include "G4DeexPrecoParameters.hh" 37 #include "G4DeexPrecoParameters.hh"
36 #include "Randomize.hh" <<
37 #include "G4Pow.hh" <<
38 #include "G4Log.hh" <<
39 #include "G4Exp.hh" <<
40 <<
41 G4VEmissionProbability::G4VEmissionProbability <<
42 : pVerbose(1), theZ(Z), theA(A), elimit(CLHE <<
43 { <<
44 pNuclearLevelData = G4NuclearLevelData::GetI <<
45 pG4pow = G4Pow::GetInstance(); <<
46 if(A > 0) { pEvapMass = G4NucleiProperties:: <<
47 G4DeexPrecoParameters* param = pNuclearLevel <<
48 OPTxs = param->GetDeexModelType(); <<
49 } <<
50 <<
51 void G4VEmissionProbability::Initialise() <<
52 { <<
53 G4DeexPrecoParameters* param = pNuclearLevel <<
54 pVerbose = param->GetVerbose(); <<
55 fFD = param->GetDiscreteExcitationFlag(); <<
56 pTolerance = param->GetMinExcitation(); <<
57 pWidth = param->GetNuclearLevelWidth(); <<
58 } <<
59 38
60 void G4VEmissionProbability::ResetIntegrator(s << 39 G4VEmissionProbability::G4VEmissionProbability()
>> 40 :OPTxs(3),useSICB(false),LevelDensity(0.1)
61 { 41 {
62 if(de > 0.0) { elimit = de; } << 42 fG4pow = G4Pow::GetInstance();
63 if(eps > 0.0) { accuracy = eps; } << 43 fPairCorr = G4PairingCorrection::GetInstance();
64 } <<
65 <<
66 G4double G4VEmissionProbability::EmissionProba <<
67 { <<
68 return 0.0; <<
69 } <<
70 <<
71 G4double G4VEmissionProbability::ComputeProbab <<
72 { <<
73 return 0.0; <<
74 } <<
75 <<
76 G4double G4VEmissionProbability::IntegrateProb <<
77 <<
78 <<
79 { <<
80 pProbability = 0.0; <<
81 if(elow >= ehigh) { return pProbability; } <<
82 <<
83 emin = elow; <<
84 emax = ehigh; <<
85 eCoulomb = cb; <<
86 <<
87 const G4double edeltamin = 0.1*CLHEP::MeV; <<
88 const G4double edeltamax = 2*CLHEP::MeV; <<
89 G4double edelta = std::min(std::min(elimit, <<
90 G4double xbin = (emax - emin)/edelta + 1.0; <<
91 G4int ibin = std::max((G4int)xbin, 4); <<
92 <<
93 // providing smart binning <<
94 G4int nbin = ibin*5; <<
95 edelta = (emax - emin)/ibin; <<
96 <<
97 G4double x(emin), y(0.0); <<
98 G4double edelmicro = edelta*0.02; <<
99 probmax = ComputeProbability(x + edelmicro, <<
100 G4double problast = probmax; <<
101 if(pVerbose > 1) { <<
102 G4cout << "### G4VEmissionProbability::Int <<
103 << "probmax=" << probmax << " Emin=" << e <<
104 << " Emax=" << emax << " QB=" << cb << " <<
105 << G4endl; <<
106 } <<
107 fE1 = fE2 = fP2 = 0.0; <<
108 G4double emax0 = emax - edelmicro; <<
109 G4bool endpoint = false; <<
110 for(G4int i=0; i<nbin; ++i) { <<
111 x += edelta; <<
112 if(x >= emax0) { <<
113 x = emax0; <<
114 endpoint = true; <<
115 } <<
116 y = ComputeProbability(x, eCoulomb); <<
117 if(pVerbose > 2) { <<
118 G4cout << " " << i << ". E= " << x < <<
119 << " Edel= " << edelta << G4endl; <<
120 } <<
121 if(y >= probmax) { <<
122 probmax = y; <<
123 } else if(0.0 == fE1 && 2*y < probmax) { <<
124 fE1 = x; <<
125 } <<
126 <<
127 G4double del = (y + problast)*edelta*0.5; <<
128 pProbability += del; <<
129 // end of the loop <<
130 if(del < accuracy*pProbability || endpoint <<
131 problast = y; <<
132 <<
133 // smart step definition <<
134 if(del != pProbability && del > 0.8*pProba <<
135 0.7*edelta > edeltamin) { <<
136 edelta *= 0.7; <<
137 } else if(del < 0.1*pProbability && 1.5*ed <<
138 edelta *= 1.5; <<
139 } <<
140 } <<
141 if(fE1 > emin && fE1 < emax) { <<
142 fE2 = std::max(0.5*(fE1 + emax), emax - ed <<
143 fP2 = 2*ComputeProbability(fE2, eCoulomb); <<
144 } <<
145 <<
146 if(pVerbose > 1) { <<
147 G4cout << " Probability= " << pProbability <<
148 << probmax << " emin=" << emin << " <<
149 << " E1=" << fE1 << " E2=" << fE2 << G4en <<
150 } <<
151 return pProbability; <<
152 } 44 }
153 45
154 G4double G4VEmissionProbability::SampleEnergy( << 46 G4VEmissionProbability::~G4VEmissionProbability()
155 { << 47 {}
156 static const G4double fact = 1.05; <<
157 static const G4double alim = 0.05; <<
158 static const G4double blim = 20.; <<
159 probmax *= fact; <<
160 <<
161 // two regions with flat and exponential maj <<
162 G4double del = emax - emin; <<
163 G4double p1 = 1.0; <<
164 G4double p2 = 0.0; <<
165 G4double a0 = 0.0; <<
166 G4double a1 = 1.0; <<
167 G4double x; <<
168 if(fE1 > 0.0 && fP2 > 0.0 && fP2 < 0.5*probm <<
169 a0 = G4Log(probmax/fP2)/(fE2 - fE1); <<
170 del= fE1 - emin; <<
171 p1 = del; <<
172 x = a0*(emax - fE1); <<
173 if(x < blim) { <<
174 a1 = (x > alim) ? 1.0 - G4Exp(-x) : x*(1 <<
175 } <<
176 p2 = a1/a0; <<
177 p1 /= (p1 + p2); <<
178 p2 = 1.0 - p1; <<
179 } <<
180 <<
181 if(pVerbose > 1) { <<
182 G4cout << "### G4VEmissionProbability::Sam <<
183 << " Emin= " << emin << " Emax= " << emax <<
184 << "/n E1=" << fE1 << " p1=" << <<
185 << " probmax=" << probmax << " P2=" << fP <<
186 } <<
187 <<
188 CLHEP::HepRandomEngine* rndm = G4Random::get <<
189 const G4int nmax = 1000; <<
190 G4double ekin, gg, gmax; <<
191 G4int n = 0; <<
192 do { <<
193 ++n; <<
194 G4double q = rndm->flat(); <<
195 if (p2 == 0.0) { <<
196 gmax = probmax; <<
197 ekin = del*q + emin; <<
198 } else if (q <= p1) { <<
199 gmax = probmax; <<
200 ekin = del*q/p1 + emin; <<
201 } else { <<
202 ekin = fE1 - G4Log(1.0 - (q - p1)*a1/p2) <<
203 x = a0*(ekin - fE1); <<
204 gmax = fP2; <<
205 if(x < blim) { <<
206 gmax = probmax*((x > alim) ? G4Exp(-x) : 1.0 <<
207 } <<
208 } <<
209 gg = ComputeProbability(ekin, eCoulomb); <<
210 if(pVerbose > 2) { <<
211 G4cout << " " << n <<
212 << ". prob= " << gg << " probmax= " << <<
213 << " Ekin= " << ekin << G4endl; <<
214 } <<
215 if((gg > gmax || n > nmax) && pVerbose > 1 <<
216 G4cout << "### G4VEmissionProbability::S <<
217 << " A= " << theA << " Eex(MeV)=" <<
218 << "\n Warning n= " << n <<
219 << " prob/gmax=" << gg/gmax <<
220 << " prob=" << gg << " gmax=" << gmax < <<
221 << "\n Ekin= " << ekin << " Emin= " <<
222 << " Emax= " << emax << G4endl; <<
223 } <<
224 } while(gmax*rndm->flat() > gg && n < nmax); <<
225 G4double enew = FindRecoilExcitation(ekin); <<
226 if(pVerbose > 1) { <<
227 G4cout << "### SampleEnergy: Efinal= " <<
228 << enew << " E=" << ekin << " Eexc=" << <<
229 } <<
230 return enew; <<
231 } <<
232 48
233 G4double G4VEmissionProbability::FindRecoilExc << 49 void G4VEmissionProbability::Initialise()
234 { 50 {
235 G4double mass = pEvapMass + fExc; << 51 G4DeexPrecoParameters* param = G4NuclearLevelData::GetInstance()->GetParameters();
236 << 52 OPTxs = param->GetDeexModelType();
237 G4double m02 = pMass*pMass; << 53 LevelDensity = param->GetLevelDensity();
238 G4double m12 = mass*mass; <<
239 G4double m22 = pResMass*pResMass; <<
240 G4double mres = std::sqrt(m02 + m12 - 2.*pMa <<
241 <<
242 fExcRes = mres - pResMass; <<
243 <<
244 if(pVerbose > 1) { <<
245 G4cout << "### FindRecoilExcitation for re <<
246 << resZ << " resA= " << resA <<
247 << " evaporated Z= " << theZ << " A <<
248 << " Ekin= " << e << " Eexc= " << fExcRes <<
249 } <<
250 <<
251 // residual nucleus is in the ground state <<
252 if(fExcRes < pTolerance) { <<
253 fExcRes = 0.0; <<
254 return std::max(0.5*(m02 + m12 - m22)/pMas <<
255 } <<
256 if(!fFD) { return e; } <<
257 <<
258 // select final state excitation <<
259 auto lManager = pNuclearLevelData->GetLevelM <<
260 if(nullptr == lManager) { return e; } <<
261 <<
262 // levels are not known <<
263 if(fExcRes > lManager->MaxLevelEnergy() + pT <<
264 <<
265 // find level <<
266 std::size_t idx = lManager->NearestLevelInde <<
267 auto level = lManager->GetLevel(idx); <<
268 <<
269 // unstable level <<
270 if (level->GetTimeGamma() == 0.0) { return e <<
271 <<
272 // is possible to use level energy? <<
273 G4double elevel = lManager->LevelEnergy(idx) <<
274 if (std::abs(elevel - fExcRes) > pWidth || p <<
275 return e; <<
276 } <<
277 <<
278 // long-lived level <<
279 G4double massR = pResMass + elevel; <<
280 G4double mr2 = massR*massR; <<
281 fExcRes = elevel; <<
282 return std::max(0.5*(m02 + m12 - mr2)/pMass <<
283 } 54 }
284 55