| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4PolynomialPDF
32 //
33 // Author: Jason Detwiler (jasonde
34 //
35 // Creation date: Aug 2012
36 // -------------------------------------------
37
38 #include "G4PolynomialPDF.hh"
39 #include "Randomize.hh"
40
41 using namespace std;
42
43 G4PolynomialPDF::G4PolynomialPDF(size_t n, con
44 G4double x1, G4double x2) :
45 fX1(x1), fX2(x2), fChanged(true), fTolerance
46 {
47 if(coeffs != nullptr) SetCoefficients(n, coe
48 else if(n > 0) SetNCoefficients(n);
49 }
50
51 G4PolynomialPDF::~G4PolynomialPDF()
52 {}
53
54 void G4PolynomialPDF::SetCoefficient(size_t i,
55 {
56 while(i >= fCoefficients.size()) fCoefficien
57 /* Loop checking, 30-Oct-2015, G.Folger */
58 fCoefficients[i] = value;
59 fChanged = true;
60 if(doSimplify) Simplify();
61 }
62
63 void G4PolynomialPDF::SetCoefficients(size_t n
64 const G4double* coefficients)
65 {
66 SetNCoefficients(nCoeffs);
67 for(size_t i=0; i<GetNCoefficients(); ++i) {
68 SetCoefficient(i, coefficients[i], false);
69 }
70 fChanged = true;
71 Simplify();
72 }
73
74 void G4PolynomialPDF::Simplify()
75 {
76 while(fCoefficients.size() && fCoefficients[
77 if(fVerbose > 0) {
78 G4cout << "G4PolynomialPDF::Simplify() W
79 << fCoefficients.size()-1 << G4en
80 }
81 fCoefficients.pop_back();
82 fChanged = true;
83 }
84 }
85
86 void G4PolynomialPDF::SetDomain(G4double x1, G
87 {
88 if(x2 <= x1) {
89 if(fVerbose > 0) {
90 G4cout << "G4PolynomialPDF::SetDomain()
91 << "(x1 = " << x1 << ", x2 = " << x2 <<
92 }
93 return;
94 }
95 fX1 = x1;
96 fX2 = x2;
97 fChanged = true;
98 }
99
100 void G4PolynomialPDF::Normalize()
101 {
102 /// Normalize PDF to 1 over domain fX1 to fX
103 /// Double-check that the highest-order coef
104 while(fCoefficients.size()) { /* Loop check
105 if(fCoefficients[fCoefficients.size()-1] =
106 else break;
107 }
108
109 G4double x1N = fX1, x2N = fX2;
110 G4double sum = 0;
111 for(size_t i=0; i<GetNCoefficients(); ++i) {
112 sum += GetCoefficient(i)*(x2N - x1N)/G4dou
113 x1N*=fX1;
114 x2N*=fX2;
115 }
116 if(sum <= 0) {
117 if(fVerbose > 0) {
118 G4cout << "G4PolynomialPDF::Normalize()
119 << sum << G4endl;
120 Dump();
121 }
122 return;
123 }
124
125 for(size_t i=0; i<GetNCoefficients(); ++i) {
126 SetCoefficient(i, GetCoefficient(i)/sum, f
127 }
128 Simplify();
129 }
130
131 G4double G4PolynomialPDF::Evaluate(G4double x,
132 {
133 /// Evaluate f(x)
134 /// ddxPower = -1: f = CDF
135 /// ddxPower = 0: f = PDF
136 /// ddxPower = 1: f = (d/dx) PDF
137 /// ddxPower = 2: f = (d2/dx2) PDF
138 if(ddxPower < -1 || ddxPower > 2) {
139 if(fVerbose > 0) {
140 G4cout << "G4PolynomialPDF::GetX() WARNI
141 << " not implemented" << G4endl;
142 }
143 return 0.0;
144 }
145
146 double f = 0.; // return value
147 double xN = 1.; // x to the power N
148 double x1N = 1.; // endpoint x1 to the power
149 for(size_t i=0; i<=GetNCoefficients(); ++i)
150 if(ddxPower == -1) { // CDF
151 if(i>0) f += GetCoefficient(i-1)*(xN - x
152 x1N *= fX1;
153 }
154 else if(ddxPower == 0 && i<GetNCoefficient
155 else if(ddxPower == 1) { // (d/dx) PDF
156 if(i<GetNCoefficients()-1) f += GetCoeff
157 }
158 else if(ddxPower == 2) { // (d2/dx2) PDF
159 if(i<GetNCoefficients()-2) f += GetCoeff
160 }
161 xN *= x;
162 }
163 return f;
164 }
165
166 G4bool G4PolynomialPDF::HasNegativeMinimum(G4d
167 {
168 // ax2 + bx + c = 0
169 // p': 2ax + b = 0 -> = 0 at min: x_extreme
170
171 if(x1 < fX1 || x2 > fX2 || x2 < x1) {
172 if(fVerbose > 0) {
173 G4cout << "G4PolynomialPDF::HasNegativeM
174 << x1 << " - " << x2 << G4endl;
175 }
176 return false;
177 }
178
179 // If flat, then check anywhere.
180 if(GetNCoefficients() == 1) return (Evaluate
181
182 // If linear, or if quadratic with negative
183 // just check the endpoints
184 if(GetNCoefficients() == 2 ||
185 (GetNCoefficients() == 3 && GetCoefficien
186 return (Evaluate(x1) < -fTolerance) || (Ev
187 }
188
189 // If quadratic and second dervative is posi
190 if(GetNCoefficients() == 3) {
191 G4double xMin = -GetCoefficient(1)*0.5/Get
192 if(xMin < x1) xMin = x1;
193 if(xMin > x2) xMin = x2;
194 return Evaluate(xMin) < -fTolerance;
195 }
196
197 // Higher-order polynomials: consider any ex
198 // are found, check the endpoints.
199 G4double extremum = GetX(0, x1, x2, 1);
200 if(Evaluate(extremum) < -fTolerance) return
201 else if(extremum <= x1+(x2-x1)*fTolerance ||
202 extremum >= x2-(x2-x1)*fTolerance) return
203 else return
204 HasNegativeMinimum(x1, extremum) || HasNega
205 }
206
207 G4double G4PolynomialPDF::GetRandomX()
208 {
209 if(fChanged) {
210 Normalize();
211 if(HasNegativeMinimum(fX1, fX2)) {
212 if(fVerbose > 0) {
213 G4cout << "G4PolynomialPDF::GetRandomX() WAR
214 << G4endl;
215 }
216 return 0.0;
217 }
218 fChanged = false;
219 }
220 return EvalInverseCDF(G4UniformRand());
221 }
222
223 G4double G4PolynomialPDF::GetX(G4double p, G4d
224 G4int ddxPower, G4double guess, G
225 {
226 /// Find a value of X between x1 and x2 at w
227 /// ddxPower = -1: f = CDF
228 /// ddxPower = 0: f = PDF
229 /// ddxPower = 1: f = (d/dx) PDF
230 /// Uses the Newton-Raphson method to find t
231 /// If not found in range, returns the neare
232
233 // input range checking
234 if(GetNCoefficients() == 0) {
235 if(fVerbose > 0) {
236 G4cout << "G4PolynomialPDF::GetX() WARNI
237 }
238 return x2;
239 }
240 if(ddxPower < -1 || ddxPower > 1) {
241 if(fVerbose > 0) {
242 G4cout << "G4PolynomialPDF::GetX() WARNI
243 << " not implemented" << G4endl;
244 }
245 return x2;
246 }
247 if(ddxPower == -1 && (p<0 || p>1)) {
248 if(fVerbose > 0) {
249 G4cout << "G4PolynomialPDF::GetX() WARNI
250 }
251 return fX2;
252 }
253
254 // check limits
255 if(x2 <= x1 || x1 < fX1 || x2 > fX2) {
256 if(fVerbose > 0) {
257 G4cout << "G4PolynomialPDF::GetX() WARNI
258 << "You sent x1 = " << x1 << ", x2 = "
259 }
260 return x2;
261 }
262
263 // Return x2 for flat lines
264 if((ddxPower == 0 && GetNCoefficients() == 1
265 (ddxPower == 1 && GetNCoefficients() == 2
266
267 // Solve p = mx + b -> x = (p-b)/m for linea
268 if((ddxPower == -1 && GetNCoefficients() ==
269 (ddxPower == 0 && GetNCoefficients() ==
270 (ddxPower == 1 && GetNCoefficients() ==
271 G4double b = (ddxPower > -1) ? GetCoeffici
272 G4double slope = GetCoefficient(ddxPower+1
273 if(slope == 0) { // the highest-order coef
274 if(fVerbose > 0) {
275 G4cout << "G4PolynomialPDF::GetX() WAR
276 << "Did you forget to Simplify(
277 }
278 return x2;
279 }
280 if(ddxPower == 1) slope *= 2.;
281 G4double value = (p-b)/slope;
282 if(value < x1) {
283 return x1;
284 }
285 else if(value > x2) {
286 return x2;
287 }
288 else {
289 return value;
290 }
291 }
292
293 // Solve quadratic equation for f-p=0 when f
294 if((ddxPower == -1 && GetNCoefficients() ==
295 (ddxPower == 0 && GetNCoefficients() ==
296 (ddxPower == 1 && GetNCoefficients() ==
297 G4double c = -p + ((ddxPower > -1) ? GetCo
298 if(ddxPower == -1) c -= (GetCoefficient(0)
299 G4double b = GetCoefficient(ddxPower+1);
300 if(ddxPower == 1) b *= 2.;
301 G4double a = GetCoefficient(ddxPower+2); /
302 if(a == 0) { // the highest-order coeffici
303 if(fVerbose > 0) {
304 G4cout << "G4PolynomialPDF::GetX() WAR
305 << "Did you forget to Simplify(
306 }
307 return x2;
308 }
309 if(ddxPower == 1) a *= 3;
310 else if(ddxPower == -1) a *= 0.5;
311 double sqrtFactor = b*b - 4.*a*c;
312 if(sqrtFactor < 0) return x2; // quadratic
313 sqrtFactor = sqrt(sqrtFactor)/2./fabs(a);
314 G4double valueMinus = -b/2./a - sqrtFactor
315 if(valueMinus >= x1 && valueMinus <= x2) r
316 else if(valueMinus > x2) return x2;
317 G4double valuePlus = -b/2./a + sqrtFactor;
318 if(valuePlus >= x1 && valuePlus <= x2) ret
319 else if(valuePlus < x1) return x2;
320 return (x1-valueMinus <= valuePlus-x2) ? x
321 }
322
323 // f is non-trivial, so use Newton-Raphson
324 // start in the middle if no good guess is p
325 if(guess < x1 || guess > x2) guess = (x2+x1)
326 G4double lastChange = 1;
327 size_t iterations = 0;
328 while(fabs(lastChange) > fTolerance) { /* L
329 // calculate f and f' simultaneously
330 G4double f = -p;
331 G4double dfdx = 0;
332 G4double xN = 1;
333 G4double x1N = 1; // only used by CDF
334 for(size_t i=0; i<=GetNCoefficients(); ++i
335 if(ddxPower == -1) { // CDF
336 if(i>0) f += GetCoefficient(i-1)*(xN -
337 if(i<GetNCoefficients()) dfdx += GetCo
338 x1N *= fX1;
339 }
340 else if(ddxPower == 0) { // PDF
341 if(i<GetNCoefficients()) f += GetCoeff
342 if(i+1<GetNCoefficients()) dfdx += Get
343 }
344 else { // ddxPower == 1: (d/dx) PDF
345 if(i+1<GetNCoefficients()) f += GetCoe
346 if(i+2<GetNCoefficients()) dfdx += Get
347 }
348 xN *= guess;
349 }
350 if(f == 0) return guess;
351 if(dfdx == 0) {
352 if(fVerbose > 0) {
353 G4cout << "G4PolynomialPDF::GetX() WARNING:
354 << ddxPower << G4endl;
355 }
356 return x2;
357 }
358 lastChange = - f/dfdx;
359
360 if(guess + lastChange < x1) {
361 lastChange = x1 - guess;
362 } else if(guess + lastChange > x2) {
363 lastChange = x2 - guess;
364 }
365
366 guess += lastChange;
367 lastChange /= (fX2-fX1);
368
369 ++iterations;
370 if(iterations > 50) {
371 if(p!=0) {
372 if(fVerbose > 0) {
373 G4cout << "G4PolynomialPDF::GetX() WARNING
374 << " between " << x1 << " and " << x2 <<
375 << ddxPower
376 << ". Last guess was " << guess << "." <<
377 }
378 }
379 if(ddxPower==-1 && bisect) {
380 if(fVerbose > 0) {
381 G4cout << "G4PolynomialPDF::GetX() WARNING
382 << G4endl;
383 }
384 return Bisect(p, x1, x2);
385 }
386 else return guess;
387 }
388 }
389 return guess;
390 }
391
392 G4double G4PolynomialPDF::Bisect( G4double p,
393 // Bisect to get 1% precision, then use Newt
394 G4double z = (x2 + x1)/2.0; // [x1 z x2]
395 if((x2 - x1)/(fX2 - fX1) < 0.01) return GetX
396 G4double fz = Evaluate(z, -1) - p;
397 if(fz < 0) return Bisect(p, z, x2); // [z x2
398 return Bisect(p, x1, z); // [x1 z]
399 }
400
401 void G4PolynomialPDF::Dump()
402 {
403 G4cout << "G4PolynomialPDF::Dump() - PDF(x)
404 for(size_t i=0; i<GetNCoefficients(); i++) {
405 if(i>0) G4cout << " + ";
406 G4cout << GetCoefficient(i);
407 if(i>0) G4cout << "*x";
408 if(i>1) G4cout << "^" << i;
409 }
410 G4cout << G4endl;
411 G4cout << "G4PolynomialPDF::Dump() - Interva
412 << fX2 << G4endl;
413 }
414
415
416