| 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 // Hadronic Process: Nuclear De-excitations 26 // Hadronic Process: Nuclear De-excitations
27 // by V. Lara (Oct 1998) 27 // by V. Lara (Oct 1998)
28 // 28 //
29 // Modifications: 29 // Modifications:
30 // 28.10.2010 V.Ivanchenko defined members in 30 // 28.10.2010 V.Ivanchenko defined members in constructor and cleaned up
31 31
32 #include "G4VEmissionProbability.hh" 32 #include "G4VEmissionProbability.hh"
33 #include "G4NuclearLevelData.hh" 33 #include "G4NuclearLevelData.hh"
34 #include "G4LevelManager.hh" 34 #include "G4LevelManager.hh"
35 #include "G4DeexPrecoParameters.hh" 35 #include "G4DeexPrecoParameters.hh"
36 #include "Randomize.hh" 36 #include "Randomize.hh"
37 #include "G4Pow.hh" 37 #include "G4Pow.hh"
38 #include "G4Log.hh" 38 #include "G4Log.hh"
39 #include "G4Exp.hh" 39 #include "G4Exp.hh"
40 40
41 G4VEmissionProbability::G4VEmissionProbability 41 G4VEmissionProbability::G4VEmissionProbability(G4int Z, G4int A)
42 : pVerbose(1), theZ(Z), theA(A), elimit(CLHE << 42 : OPTxs(3), pVerbose(1), theZ(Z), theA(A), elimit(CLHEP::MeV)
43 { 43 {
44 pNuclearLevelData = G4NuclearLevelData::GetI 44 pNuclearLevelData = G4NuclearLevelData::GetInstance();
45 pG4pow = G4Pow::GetInstance(); 45 pG4pow = G4Pow::GetInstance();
46 if(A > 0) { pEvapMass = G4NucleiProperties:: 46 if(A > 0) { pEvapMass = G4NucleiProperties::GetNuclearMass(theA, theZ); }
47 G4DeexPrecoParameters* param = pNuclearLevel <<
48 OPTxs = param->GetDeexModelType(); <<
49 } 47 }
50 48
51 void G4VEmissionProbability::Initialise() 49 void G4VEmissionProbability::Initialise()
52 { 50 {
53 G4DeexPrecoParameters* param = pNuclearLevel 51 G4DeexPrecoParameters* param = pNuclearLevelData->GetParameters();
>> 52 OPTxs = param->GetDeexModelType();
54 pVerbose = param->GetVerbose(); 53 pVerbose = param->GetVerbose();
55 fFD = param->GetDiscreteExcitationFlag(); 54 fFD = param->GetDiscreteExcitationFlag();
56 pTolerance = param->GetMinExcitation(); 55 pTolerance = param->GetMinExcitation();
57 pWidth = param->GetNuclearLevelWidth(); <<
58 } 56 }
59 57
60 void G4VEmissionProbability::ResetIntegrator(s 58 void G4VEmissionProbability::ResetIntegrator(size_t, G4double de, G4double eps)
61 { 59 {
62 if(de > 0.0) { elimit = de; } 60 if(de > 0.0) { elimit = de; }
63 if(eps > 0.0) { accuracy = eps; } 61 if(eps > 0.0) { accuracy = eps; }
64 } 62 }
65 63
66 G4double G4VEmissionProbability::EmissionProba 64 G4double G4VEmissionProbability::EmissionProbability(const G4Fragment&, G4double)
67 { 65 {
68 return 0.0; 66 return 0.0;
69 } 67 }
70 68
71 G4double G4VEmissionProbability::ComputeProbab 69 G4double G4VEmissionProbability::ComputeProbability(G4double, G4double)
72 { 70 {
73 return 0.0; 71 return 0.0;
74 } 72 }
75 73
76 G4double G4VEmissionProbability::IntegrateProb 74 G4double G4VEmissionProbability::IntegrateProbability(G4double elow,
77 75 G4double ehigh,
78 76 G4double cb)
79 { 77 {
80 pProbability = 0.0; 78 pProbability = 0.0;
81 if(elow >= ehigh) { return pProbability; } 79 if(elow >= ehigh) { return pProbability; }
82 80
83 emin = elow; 81 emin = elow;
84 emax = ehigh; 82 emax = ehigh;
85 eCoulomb = cb; 83 eCoulomb = cb;
86 84
87 const G4double edeltamin = 0.1*CLHEP::MeV; 85 const G4double edeltamin = 0.1*CLHEP::MeV;
88 const G4double edeltamax = 2*CLHEP::MeV; 86 const G4double edeltamax = 2*CLHEP::MeV;
89 G4double edelta = std::min(std::min(elimit, 87 G4double edelta = std::min(std::min(elimit, edeltamax), edeltamin);
90 G4double xbin = (emax - emin)/edelta + 1.0; 88 G4double xbin = (emax - emin)/edelta + 1.0;
91 G4int ibin = std::max((G4int)xbin, 4); 89 G4int ibin = std::max((G4int)xbin, 4);
92 90
93 // providing smart binning 91 // providing smart binning
94 G4int nbin = ibin*5; 92 G4int nbin = ibin*5;
95 edelta = (emax - emin)/ibin; 93 edelta = (emax - emin)/ibin;
96 94
97 G4double x(emin), y(0.0); 95 G4double x(emin), y(0.0);
98 G4double edelmicro = edelta*0.02; 96 G4double edelmicro = edelta*0.02;
99 probmax = ComputeProbability(x + edelmicro, 97 probmax = ComputeProbability(x + edelmicro, eCoulomb);
100 G4double problast = probmax; 98 G4double problast = probmax;
101 if(pVerbose > 1) { 99 if(pVerbose > 1) {
102 G4cout << "### G4VEmissionProbability::Int 100 G4cout << "### G4VEmissionProbability::IntegrateProbability: "
103 << "probmax=" << probmax << " Emin=" << e 101 << "probmax=" << probmax << " Emin=" << emin
104 << " Emax=" << emax << " QB=" << cb << " 102 << " Emax=" << emax << " QB=" << cb << " nbin=" << nbin
105 << G4endl; 103 << G4endl;
106 } 104 }
107 fE1 = fE2 = fP2 = 0.0; 105 fE1 = fE2 = fP2 = 0.0;
108 G4double emax0 = emax - edelmicro; 106 G4double emax0 = emax - edelmicro;
109 G4bool endpoint = false; 107 G4bool endpoint = false;
110 for(G4int i=0; i<nbin; ++i) { 108 for(G4int i=0; i<nbin; ++i) {
111 x += edelta; 109 x += edelta;
112 if(x >= emax0) { 110 if(x >= emax0) {
113 x = emax0; 111 x = emax0;
114 endpoint = true; 112 endpoint = true;
115 } 113 }
116 y = ComputeProbability(x, eCoulomb); 114 y = ComputeProbability(x, eCoulomb);
117 if(pVerbose > 2) { 115 if(pVerbose > 2) {
118 G4cout << " " << i << ". E= " << x < 116 G4cout << " " << i << ". E= " << x << " prob= " << y
119 << " Edel= " << edelta << G4endl; 117 << " Edel= " << edelta << G4endl;
120 } 118 }
121 if(y >= probmax) { 119 if(y >= probmax) {
122 probmax = y; 120 probmax = y;
123 } else if(0.0 == fE1 && 2*y < probmax) { 121 } else if(0.0 == fE1 && 2*y < probmax) {
124 fE1 = x; 122 fE1 = x;
125 } 123 }
126 124
127 G4double del = (y + problast)*edelta*0.5; 125 G4double del = (y + problast)*edelta*0.5;
128 pProbability += del; 126 pProbability += del;
129 // end of the loop 127 // end of the loop
130 if(del < accuracy*pProbability || endpoint 128 if(del < accuracy*pProbability || endpoint) { break; }
131 problast = y; 129 problast = y;
132 130
133 // smart step definition 131 // smart step definition
134 if(del != pProbability && del > 0.8*pProba 132 if(del != pProbability && del > 0.8*pProbability &&
135 0.7*edelta > edeltamin) { 133 0.7*edelta > edeltamin) {
136 edelta *= 0.7; 134 edelta *= 0.7;
137 } else if(del < 0.1*pProbability && 1.5*ed 135 } else if(del < 0.1*pProbability && 1.5*edelta < edeltamax) {
138 edelta *= 1.5; 136 edelta *= 1.5;
139 } 137 }
140 } 138 }
141 if(fE1 > emin && fE1 < emax) { 139 if(fE1 > emin && fE1 < emax) {
142 fE2 = std::max(0.5*(fE1 + emax), emax - ed 140 fE2 = std::max(0.5*(fE1 + emax), emax - edelta);
143 fP2 = 2*ComputeProbability(fE2, eCoulomb); 141 fP2 = 2*ComputeProbability(fE2, eCoulomb);
144 } 142 }
145 143
146 if(pVerbose > 1) { 144 if(pVerbose > 1) {
147 G4cout << " Probability= " << pProbability 145 G4cout << " Probability= " << pProbability << " probmax= "
148 << probmax << " emin=" << emin << " 146 << probmax << " emin=" << emin << " emax=" << emax
149 << " E1=" << fE1 << " E2=" << fE2 << G4en 147 << " E1=" << fE1 << " E2=" << fE2 << G4endl;
150 } 148 }
151 return pProbability; 149 return pProbability;
152 } 150 }
153 151
154 G4double G4VEmissionProbability::SampleEnergy( 152 G4double G4VEmissionProbability::SampleEnergy()
155 { 153 {
156 static const G4double fact = 1.05; 154 static const G4double fact = 1.05;
157 static const G4double alim = 0.05; 155 static const G4double alim = 0.05;
158 static const G4double blim = 20.; 156 static const G4double blim = 20.;
159 probmax *= fact; 157 probmax *= fact;
160 158
161 // two regions with flat and exponential maj 159 // two regions with flat and exponential majorant
162 G4double del = emax - emin; 160 G4double del = emax - emin;
163 G4double p1 = 1.0; 161 G4double p1 = 1.0;
164 G4double p2 = 0.0; 162 G4double p2 = 0.0;
165 G4double a0 = 0.0; 163 G4double a0 = 0.0;
166 G4double a1 = 1.0; 164 G4double a1 = 1.0;
167 G4double x; 165 G4double x;
168 if(fE1 > 0.0 && fP2 > 0.0 && fP2 < 0.5*probm 166 if(fE1 > 0.0 && fP2 > 0.0 && fP2 < 0.5*probmax) {
169 a0 = G4Log(probmax/fP2)/(fE2 - fE1); 167 a0 = G4Log(probmax/fP2)/(fE2 - fE1);
170 del= fE1 - emin; 168 del= fE1 - emin;
171 p1 = del; 169 p1 = del;
172 x = a0*(emax - fE1); 170 x = a0*(emax - fE1);
173 if(x < blim) { 171 if(x < blim) {
174 a1 = (x > alim) ? 1.0 - G4Exp(-x) : x*(1 172 a1 = (x > alim) ? 1.0 - G4Exp(-x) : x*(1.0 - 0.5*x);
175 } 173 }
176 p2 = a1/a0; 174 p2 = a1/a0;
177 p1 /= (p1 + p2); 175 p1 /= (p1 + p2);
178 p2 = 1.0 - p1; 176 p2 = 1.0 - p1;
179 } 177 }
180 178
181 if(pVerbose > 1) { 179 if(pVerbose > 1) {
182 G4cout << "### G4VEmissionProbability::Sam 180 G4cout << "### G4VEmissionProbability::SampleEnergy: "
183 << " Emin= " << emin << " Emax= " << emax 181 << " Emin= " << emin << " Emax= " << emax
184 << "/n E1=" << fE1 << " p1=" << 182 << "/n E1=" << fE1 << " p1=" << p1
185 << " probmax=" << probmax << " P2=" << fP 183 << " probmax=" << probmax << " P2=" << fP2 << G4endl;
186 } 184 }
187 185
188 CLHEP::HepRandomEngine* rndm = G4Random::get 186 CLHEP::HepRandomEngine* rndm = G4Random::getTheEngine();
189 const G4int nmax = 1000; 187 const G4int nmax = 1000;
190 G4double ekin, gg, gmax; << 188 G4double ekin, g, gmax;
191 G4int n = 0; 189 G4int n = 0;
192 do { 190 do {
193 ++n; 191 ++n;
194 G4double q = rndm->flat(); 192 G4double q = rndm->flat();
195 if (p2 == 0.0) { << 193 if(q <= p1) {
196 gmax = probmax; <<
197 ekin = del*q + emin; <<
198 } else if (q <= p1) { <<
199 gmax = probmax; 194 gmax = probmax;
200 ekin = del*q/p1 + emin; 195 ekin = del*q/p1 + emin;
201 } else { 196 } else {
202 ekin = fE1 - G4Log(1.0 - (q - p1)*a1/p2) 197 ekin = fE1 - G4Log(1.0 - (q - p1)*a1/p2)/a0;
203 x = a0*(ekin - fE1); 198 x = a0*(ekin - fE1);
204 gmax = fP2; 199 gmax = fP2;
205 if(x < blim) { 200 if(x < blim) {
206 gmax = probmax*((x > alim) ? G4Exp(-x) : 1.0 201 gmax = probmax*((x > alim) ? G4Exp(-x) : 1.0 - x*(1.0 - 0.5*x));
207 } 202 }
208 } 203 }
209 gg = ComputeProbability(ekin, eCoulomb); << 204 g = ComputeProbability(ekin, eCoulomb);
210 if(pVerbose > 2) { 205 if(pVerbose > 2) {
211 G4cout << " " << n 206 G4cout << " " << n
212 << ". prob= " << gg << " probmax= " << << 207 << ". prob= " << g << " probmax= " << probmax
213 << " Ekin= " << ekin << G4endl; 208 << " Ekin= " << ekin << G4endl;
214 } 209 }
215 if((gg > gmax || n > nmax) && pVerbose > 1 << 210 if((g > gmax || n > nmax) && pVerbose > 1) {
216 G4cout << "### G4VEmissionProbability::S 211 G4cout << "### G4VEmissionProbability::SampleEnergy for Z= " << theZ
217 << " A= " << theA << " Eex(MeV)=" 212 << " A= " << theA << " Eex(MeV)=" << fExc << " p1=" << p1
218 << "\n Warning n= " << n 213 << "\n Warning n= " << n
219 << " prob/gmax=" << gg/gmax << 214 << " prob/gmax=" << g/gmax
220 << " prob=" << gg << " gmax=" << gmax < << 215 << " prob=" << g << " gmax=" << gmax << " probmax=" << probmax
221 << "\n Ekin= " << ekin << " Emin= " 216 << "\n Ekin= " << ekin << " Emin= " << emin
222 << " Emax= " << emax << G4endl; 217 << " Emax= " << emax << G4endl;
223 } 218 }
224 } while(gmax*rndm->flat() > gg && n < nmax); << 219 } while(gmax*rndm->flat() > g && n < nmax);
225 G4double enew = FindRecoilExcitation(ekin); 220 G4double enew = FindRecoilExcitation(ekin);
226 if(pVerbose > 1) { 221 if(pVerbose > 1) {
227 G4cout << "### SampleEnergy: Efinal= " 222 G4cout << "### SampleEnergy: Efinal= "
228 << enew << " E=" << ekin << " Eexc=" << 223 << enew << " E=" << ekin << " Eexc=" << fExcRes << G4endl;
229 } 224 }
230 return enew; 225 return enew;
231 } 226 }
232 227
233 G4double G4VEmissionProbability::FindRecoilExc 228 G4double G4VEmissionProbability::FindRecoilExcitation(const G4double e)
234 { 229 {
235 G4double mass = pEvapMass + fExc; 230 G4double mass = pEvapMass + fExc;
236 231
237 G4double m02 = pMass*pMass; 232 G4double m02 = pMass*pMass;
238 G4double m12 = mass*mass; 233 G4double m12 = mass*mass;
239 G4double m22 = pResMass*pResMass; 234 G4double m22 = pResMass*pResMass;
240 G4double mres = std::sqrt(m02 + m12 - 2.*pMa 235 G4double mres = std::sqrt(m02 + m12 - 2.*pMass*(mass + e));
241 236
242 fExcRes = mres - pResMass; 237 fExcRes = mres - pResMass;
243 238
244 if(pVerbose > 1) { 239 if(pVerbose > 1) {
245 G4cout << "### FindRecoilExcitation for re 240 G4cout << "### FindRecoilExcitation for resZ= "
246 << resZ << " resA= " << resA 241 << resZ << " resA= " << resA
247 << " evaporated Z= " << theZ << " A 242 << " evaporated Z= " << theZ << " A= " << theA
248 << " Ekin= " << e << " Eexc= " << fExcRes 243 << " Ekin= " << e << " Eexc= " << fExcRes << G4endl;
249 } 244 }
250 245
251 // residual nucleus is in the ground state 246 // residual nucleus is in the ground state
252 if(fExcRes < pTolerance) { 247 if(fExcRes < pTolerance) {
253 fExcRes = 0.0; 248 fExcRes = 0.0;
254 return std::max(0.5*(m02 + m12 - m22)/pMas 249 return std::max(0.5*(m02 + m12 - m22)/pMass - mass, 0.0);
255 } 250 }
256 if(!fFD) { return e; } 251 if(!fFD) { return e; }
257 252
258 // select final state excitation 253 // select final state excitation
259 auto lManager = pNuclearLevelData->GetLevelM 254 auto lManager = pNuclearLevelData->GetLevelManager(resZ, resA);
260 if(nullptr == lManager) { return e; } 255 if(nullptr == lManager) { return e; }
261 256
262 // levels are not known 257 // levels are not known
263 if(fExcRes > lManager->MaxLevelEnergy() + pT 258 if(fExcRes > lManager->MaxLevelEnergy() + pTolerance) { return e; }
264 259
265 // find level 260 // find level
266 std::size_t idx = lManager->NearestLevelInde << 261 G4double elevel = lManager->NearestLevelEnergy(fExcRes);
267 auto level = lManager->GetLevel(idx); <<
268 <<
269 // unstable level <<
270 if (level->GetTimeGamma() == 0.0) { return e <<
271 262
272 // is possible to use level energy? << 263 // excited level
273 G4double elevel = lManager->LevelEnergy(idx) << 264 if(pMass > mass + pResMass + elevel &&
274 if (std::abs(elevel - fExcRes) > pWidth || p << 265 std::abs(elevel - fExcRes) <= pTolerance) {
275 return e; << 266 G4double massR = pResMass + elevel;
276 } << 267 G4double mr2 = massR*massR;
277 << 268 fExcRes = elevel;
278 // long-lived level << 269 return std::max(0.5*(m02 + m12 - mr2)/pMass - mass, 0.0);
279 G4double massR = pResMass + elevel; << 270 }
280 G4double mr2 = massR*massR; << 271 return e;
281 fExcRes = elevel; <<
282 return std::max(0.5*(m02 + m12 - mr2)/pMass <<
283 } 272 }
284 273