Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolynomialPDF.cc
Go to the documentation of this file.
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 // -------------------------------------------------------------------
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4PolynomialPDF
32 //
33 // Author: Jason Detwiler (jasondet@gmail.com)
34 //
35 // Creation date: Aug 2012
36 // -------------------------------------------------------------------
37 
38 #include "G4PolynomialPDF.hh"
39 #include "Randomize.hh"
40 
41 using namespace std;
42 
44  G4double x1, G4double x2) :
45  fX1(x1), fX2(x2), fChanged(true), fTolerance(1.e-8), fVerbose(0)
46 {
47  if(coeffs != nullptr) SetCoefficients(n, coeffs);
48  else if(n > 0) SetNCoefficients(n);
49 }
50 
52 {}
53 
54 void G4PolynomialPDF::SetCoefficient(size_t i, G4double value, bool doSimplify)
55 {
56  while(i >= fCoefficients.size()) fCoefficients.push_back(0);
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 nCoeffs,
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 
75 {
76  while(fCoefficients.size() && fCoefficients[fCoefficients.size()-1] == 0) {
77  if(fVerbose > 0) {
78  G4cout << "G4PolynomialPDF::Simplify() WARNING: had to pop coefficient "
79  << fCoefficients.size()-1 << G4endl;
80  }
81  fCoefficients.pop_back();
82  fChanged = true;
83  }
84 }
85 
87 {
88  if(x2 <= x1) {
89  if(fVerbose > 0) {
90  G4cout << "G4PolynomialPDF::SetDomain() WARNING: Invalide domain! "
91  << "(x1 = " << x1 << ", x2 = " << x2 << ")." << G4endl;
92  }
93  return;
94  }
95  fX1 = x1;
96  fX2 = x2;
97  fChanged = true;
98 }
99 
101 {
104  while(fCoefficients.size()) { /* Loop checking, 30-Oct-2015, G.Folger */
105  if(fCoefficients[fCoefficients.size()-1] == 0.0) fCoefficients.pop_back();
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)/G4double(i+1);
113  x1N*=fX1;
114  x2N*=fX2;
115  }
116  if(sum <= 0) {
117  if(fVerbose > 0) {
118  G4cout << "G4PolynomialPDF::Normalize() WARNING: PDF has non-positive area: "
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, false);
127  }
128  Simplify();
129 }
130 
132 {
138  if(ddxPower < -1 || ddxPower > 2) {
139  if(fVerbose > 0) {
140  G4cout << "G4PolynomialPDF::GetX() WARNING: ddxPower " << ddxPower
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 N; only used by CDF
149  for(size_t i=0; i<=GetNCoefficients(); ++i) {
150  if(ddxPower == -1) { // CDF
151  if(i>0) f += GetCoefficient(i-1)*(xN - x1N)/i;
152  x1N *= fX1;
153  }
154  else if(ddxPower == 0 && i<GetNCoefficients()) f += GetCoefficient(i)*xN; // PDF
155  else if(ddxPower == 1) { // (d/dx) PDF
156  if(i<GetNCoefficients()-1) f += GetCoefficient(i+1)*xN*(i+1);
157  }
158  else if(ddxPower == 2) { // (d2/dx2) PDF
159  if(i<GetNCoefficients()-2) f += GetCoefficient(i+2)*xN*((i+2)*(i+1));
160  }
161  xN *= x;
162  }
163  return f;
164 }
165 
167 {
168  // ax2 + bx + c = 0
169  // p': 2ax + b = 0 -> = 0 at min: x_extreme = -b/2a
170 
171  if(x1 < fX1 || x2 > fX2 || x2 < x1) {
172  if(fVerbose > 0) {
173  G4cout << "G4PolynomialPDF::HasNegativeMinimum() WARNING: Invalid range "
174  << x1 << " - " << x2 << G4endl;
175  }
176  return false;
177  }
178 
179  // If flat, then check anywhere.
180  if(GetNCoefficients() == 1) return (Evaluate(x1) < -fTolerance);
181 
182  // If linear, or if quadratic with negative second derivative,
183  // just check the endpoints
184  if(GetNCoefficients() == 2 ||
185  (GetNCoefficients() == 3 && GetCoefficient(2) <= 0)) {
186  return (Evaluate(x1) < -fTolerance) || (Evaluate(x2) < -fTolerance);
187  }
188 
189  // If quadratic and second dervative is positive, check at the mininum
190  if(GetNCoefficients() == 3) {
191  G4double xMin = -GetCoefficient(1)*0.5/GetCoefficient(2);
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 extremum between x1 and x2. If none
198  // are found, check the endpoints.
199  G4double extremum = GetX(0, x1, x2, 1);
200  if(Evaluate(extremum) < -fTolerance) return true;
201  else if(extremum <= x1+(x2-x1)*fTolerance ||
202  extremum >= x2-(x2-x1)*fTolerance) return false;
203  else return
204  HasNegativeMinimum(x1, extremum) || HasNegativeMinimum(extremum, x2);
205 }
206 
208 {
209  if(fChanged) {
210  Normalize();
211  if(HasNegativeMinimum(fX1, fX2)) {
212  if(fVerbose > 0) {
213  G4cout << "G4PolynomialPDF::GetRandomX() WARNING: PDF has negative values, returning 0..."
214  << G4endl;
215  }
216  return 0.0;
217  }
218  fChanged = false;
219  }
220  return EvalInverseCDF(G4UniformRand());
221 }
222 
224  G4int ddxPower, G4double guess, G4bool bisect)
225 {
232 
233  // input range checking
234  if(GetNCoefficients() == 0) {
235  if(fVerbose > 0) {
236  G4cout << "G4PolynomialPDF::GetX() WARNING: no PDF defined!" << G4endl;
237  }
238  return x2;
239  }
240  if(ddxPower < -1 || ddxPower > 1) {
241  if(fVerbose > 0) {
242  G4cout << "G4PolynomialPDF::GetX() WARNING: ddxPower " << ddxPower
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() WARNING: p is out of range" << G4endl;
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() WARNING: domain must have fX1 <= x1 < x2 <= fX2. "
258  << "You sent x1 = " << x1 << ", x2 = " << x2 << "." << G4endl;
259  }
260  return x2;
261  }
262 
263  // Return x2 for flat lines
264  if((ddxPower == 0 && GetNCoefficients() == 1) ||
265  (ddxPower == 1 && GetNCoefficients() == 2)) return x2;
266 
267  // Solve p = mx + b -> x = (p-b)/m for linear functions
268  if((ddxPower == -1 && GetNCoefficients() == 1) ||
269  (ddxPower == 0 && GetNCoefficients() == 2) ||
270  (ddxPower == 1 && GetNCoefficients() == 3)) {
271  G4double b = (ddxPower > -1) ? GetCoefficient(ddxPower) : -GetCoefficient(0)*fX1;
272  G4double slope = GetCoefficient(ddxPower+1); // the highest-order coefficient
273  if(slope == 0) { // the highest-order coefficient should never be zero if simplified
274  if(fVerbose > 0) {
275  G4cout << "G4PolynomialPDF::GetX() WARNING: Got slope = 0. "
276  << "Did you forget to Simplify()?" << G4endl;
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 is quadratic
294  if((ddxPower == -1 && GetNCoefficients() == 2) ||
295  (ddxPower == 0 && GetNCoefficients() == 3) ||
296  (ddxPower == 1 && GetNCoefficients() == 4)) {
297  G4double c = -p + ((ddxPower > -1) ? GetCoefficient(ddxPower) : 0);
298  if(ddxPower == -1) c -= (GetCoefficient(0) + GetCoefficient(1)/2.*fX1)*fX1;
299  G4double b = GetCoefficient(ddxPower+1);
300  if(ddxPower == 1) b *= 2.;
301  G4double a = GetCoefficient(ddxPower+2); // the highest-order coefficient
302  if(a == 0) { // the highest-order coefficient should never be 0 if simplified
303  if(fVerbose > 0) {
304  G4cout << "G4PolynomialPDF::GetX() WARNING: Got a = 0. "
305  << "Did you forget to Simplify()?" << G4endl;
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 equation has no solution (p not in range of f)
313  sqrtFactor = sqrt(sqrtFactor)/2./fabs(a);
314  G4double valueMinus = -b/2./a - sqrtFactor;
315  if(valueMinus >= x1 && valueMinus <= x2) return valueMinus;
316  else if(valueMinus > x2) return x2;
317  G4double valuePlus = -b/2./a + sqrtFactor;
318  if(valuePlus >= x1 && valuePlus <= x2) return valuePlus;
319  else if(valuePlus < x1) return x2;
320  return (x1-valueMinus <= valuePlus-x2) ? x1 : x2;
321  }
322 
323  // f is non-trivial, so use Newton-Raphson
324  // start in the middle if no good guess is provided
325  if(guess < x1 || guess > x2) guess = (x2+x1)*0.5;
326  G4double lastChange = 1;
327  size_t iterations = 0;
328  while(fabs(lastChange) > fTolerance) { /* Loop checking, 02.11.2015, A.Ribon */
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 - x1N)/G4double(i);
337  if(i<GetNCoefficients()) dfdx += GetCoefficient(i)*xN;
338  x1N *= fX1;
339  }
340  else if(ddxPower == 0) { // PDF
341  if(i<GetNCoefficients()) f += GetCoefficient(i)*xN;
342  if(i+1<GetNCoefficients()) dfdx += GetCoefficient(i+1)*xN*G4double(i+1);
343  }
344  else { // ddxPower == 1: (d/dx) PDF
345  if(i+1<GetNCoefficients()) f += GetCoefficient(i+1)*xN*G4double(i+1);
346  if(i+2<GetNCoefficients()) dfdx += GetCoefficient(i+2)*xN*G4double(i+2);
347  }
348  xN *= guess;
349  }
350  if(f == 0) return guess;
351  if(dfdx == 0) {
352  if(fVerbose > 0) {
353  G4cout << "G4PolynomialPDF::GetX() WARNING: got f != 0 but slope = 0 for ddxPower = "
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: got stuck searching for " << p
374  << " between " << x1 << " and " << x2 << " with ddxPower = "
375  << ddxPower
376  << ". Last guess was " << guess << "." << G4endl;
377  }
378  }
379  if(ddxPower==-1 && bisect) {
380  if(fVerbose > 0) {
381  G4cout << "G4PolynomialPDF::GetX() WARNING: Bisceting and trying again..."
382  << G4endl;
383  }
384  return Bisect(p, x1, x2);
385  }
386  else return guess;
387  }
388  }
389  return guess;
390 }
391 
393  // Bisect to get 1% precision, then use Newton-Raphson
394  G4double z = (x2 + x1)/2.0; // [x1 z x2]
395  if((x2 - x1)/(fX2 - fX1) < 0.01) return GetX(p, fX1, fX2, -1, z, false);
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 
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() - Interval: " << fX1 << " <= x < "
412  << fX2 << G4endl;
413 }
414 
415 
std::vector< G4double > fCoefficients
void SetCoefficients(const std::vector< G4double > &v)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4double GetRandomX()
const char * p
Definition: xmltok.h:285
void SetNCoefficients(size_t n)
G4PolynomialPDF(size_t n=0, const double *coeffs=nullptr, G4double x1=0, G4double x2=1)
tuple x
Definition: test.py:50
G4double GetX(G4double p, G4double x1, G4double x2, G4int ddxPower=0, G4double guess=1.e99, G4bool bisect=true)
int G4int
Definition: G4Types.hh:78
G4double Bisect(G4double p, G4double x1, G4double x2)
tuple b
Definition: test.py:12
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double GetCoefficient(size_t i) const
const XML_Char int const XML_Char * value
Definition: expat.h:331
void SetDomain(G4double x1, G4double x2)
bool G4bool
Definition: G4Types.hh:79
const G4int n
G4bool HasNegativeMinimum(G4double x1, G4double x2)
size_t GetNCoefficients() const
tuple z
Definition: test.py:28
#define G4endl
Definition: G4ios.hh:61
G4double EvalInverseCDF(G4double p)
G4double Evaluate(G4double x, G4int ddxPower=0)
double G4double
Definition: G4Types.hh:76
void SetCoefficient(size_t i, G4double value, bool doSimplify)
tuple c
Definition: test.py:13