Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolynomialPDF Class Reference

#include <G4PolynomialPDF.hh>

Public Member Functions

 G4PolynomialPDF (size_t n=0, const double *coeffs=nullptr, G4double x1=0, G4double x2=1)
 
 ~G4PolynomialPDF ()
 
void SetNCoefficients (size_t n)
 
size_t GetNCoefficients () const
 
void SetCoefficients (const std::vector< G4double > &v)
 
G4double GetCoefficient (size_t i) const
 
void SetCoefficient (size_t i, G4double value, bool doSimplify)
 
void SetCoefficients (size_t n, const G4double *coeffs)
 
void Simplify ()
 
void SetDomain (G4double x1, G4double x2)
 
void Normalize ()
 
G4double Evaluate (G4double x, G4int ddxPower=0)
 
G4double GetRandomX ()
 
void SetTolerance (G4double tolerance)
 
G4double GetX (G4double p, G4double x1, G4double x2, G4int ddxPower=0, G4double guess=1.e99, G4bool bisect=true)
 
G4double EvalInverseCDF (G4double p)
 
G4double Bisect (G4double p, G4double x1, G4double x2)
 
void Dump ()
 

Protected Member Functions

G4bool HasNegativeMinimum (G4double x1, G4double x2)
 

Protected Attributes

G4double fX1
 
G4double fX2
 
std::vector< G4doublefCoefficients
 
G4bool fChanged
 
G4double fTolerance
 
G4int fVerbose
 

Detailed Description

Definition at line 49 of file G4PolynomialPDF.hh.

Constructor & Destructor Documentation

G4PolynomialPDF::G4PolynomialPDF ( size_t  n = 0,
const double *  coeffs = nullptr,
G4double  x1 = 0,
G4double  x2 = 1 
)

Definition at line 43 of file G4PolynomialPDF.cc.

44  :
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 }
void SetCoefficients(const std::vector< G4double > &v)
void SetNCoefficients(size_t n)
const G4int n

Here is the call graph for this function:

G4PolynomialPDF::~G4PolynomialPDF ( )

Definition at line 51 of file G4PolynomialPDF.cc.

52 {}

Member Function Documentation

G4double G4PolynomialPDF::Bisect ( G4double  p,
G4double  x1,
G4double  x2 
)

Definition at line 392 of file G4PolynomialPDF.cc.

392  {
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 }
const char * p
Definition: xmltok.h:285
G4double GetX(G4double p, G4double x1, G4double x2, G4int ddxPower=0, G4double guess=1.e99, G4bool bisect=true)
G4double Bisect(G4double p, G4double x1, G4double x2)
tuple z
Definition: test.py:28
G4double Evaluate(G4double x, G4int ddxPower=0)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4PolynomialPDF::Dump ( )

Definition at line 401 of file G4PolynomialPDF.cc.

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 }
G4GLOB_DLL std::ostream G4cout
G4double GetCoefficient(size_t i) const
size_t GetNCoefficients() const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PolynomialPDF::EvalInverseCDF ( G4double  p)
inline

Definition at line 102 of file G4PolynomialPDF.hh.

102 { return GetX(p, fX1, fX2, -1, fX1 + p*(fX2-fX1)); }
const char * p
Definition: xmltok.h:285
G4double GetX(G4double p, G4double x1, G4double x2, G4int ddxPower=0, G4double guess=1.e99, G4bool bisect=true)

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PolynomialPDF::Evaluate ( G4double  x,
G4int  ddxPower = 0 
)

Evaluate f(x) ddxPower = -1: f = CDF ddxPower = 0: f = PDF ddxPower = 1: f = (d/dx) PDF ddxPower = 2: f = (d2/dx2) PDF

Definition at line 131 of file G4PolynomialPDF.cc.

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 }
tuple x
Definition: test.py:50
G4GLOB_DLL std::ostream G4cout
G4double GetCoefficient(size_t i) const
size_t GetNCoefficients() const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PolynomialPDF::GetCoefficient ( size_t  i) const
inline

Definition at line 62 of file G4PolynomialPDF.hh.

62 { return fCoefficients[i]; }
std::vector< G4double > fCoefficients

Here is the caller graph for this function:

size_t G4PolynomialPDF::GetNCoefficients ( ) const
inline

Definition at line 58 of file G4PolynomialPDF.hh.

58 { return fCoefficients.size(); }
std::vector< G4double > fCoefficients

Here is the caller graph for this function:

G4double G4PolynomialPDF::GetRandomX ( )

Definition at line 207 of file G4PolynomialPDF.cc.

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 }
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4bool HasNegativeMinimum(G4double x1, G4double x2)
#define G4endl
Definition: G4ios.hh:61
G4double EvalInverseCDF(G4double p)

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PolynomialPDF::GetX ( G4double  p,
G4double  x1,
G4double  x2,
G4int  ddxPower = 0,
G4double  guess = 1.e99,
G4bool  bisect = true 
)

Find a value of X between x1 and x2 at which f(x) = p. ddxPower = -1: f = CDF ddxPower = 0: f = PDF ddxPower = 1: f = (d/dx) PDF Uses the Newton-Raphson method to find the zero of f(x) - p. If not found in range, returns the nearest boundary

Definition at line 223 of file G4PolynomialPDF.cc.

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 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const char * p
Definition: xmltok.h:285
G4double Bisect(G4double p, G4double x1, G4double x2)
tuple b
Definition: test.py:12
G4GLOB_DLL std::ostream G4cout
G4double GetCoefficient(size_t i) const
const XML_Char int const XML_Char * value
Definition: expat.h:331
size_t GetNCoefficients() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13

Here is the call graph for this function:

Here is the caller graph for this function:

G4bool G4PolynomialPDF::HasNegativeMinimum ( G4double  x1,
G4double  x2 
)
protected

Definition at line 166 of file G4PolynomialPDF.cc.

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 }
G4double GetX(G4double p, G4double x1, G4double x2, G4int ddxPower=0, G4double guess=1.e99, G4bool bisect=true)
G4GLOB_DLL std::ostream G4cout
G4double GetCoefficient(size_t i) const
G4bool HasNegativeMinimum(G4double x1, G4double x2)
size_t GetNCoefficients() const
#define G4endl
Definition: G4ios.hh:61
G4double Evaluate(G4double x, G4int ddxPower=0)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4PolynomialPDF::Normalize ( )

Normalize PDF to 1 over domain fX1 to fX2. Double-check that the highest-order coefficient is non-zero.

Definition at line 100 of file G4PolynomialPDF.cc.

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 }
std::vector< G4double > fCoefficients
G4GLOB_DLL std::ostream G4cout
G4double GetCoefficient(size_t i) const
size_t GetNCoefficients() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void SetCoefficient(size_t i, G4double value, bool doSimplify)

Here is the call graph for this function:

Here is the caller graph for this function:

void G4PolynomialPDF::SetCoefficient ( size_t  i,
G4double  value,
bool  doSimplify 
)

Definition at line 54 of file G4PolynomialPDF.cc.

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 }
std::vector< G4double > fCoefficients
const XML_Char int const XML_Char * value
Definition: expat.h:331

Here is the call graph for this function:

Here is the caller graph for this function:

void G4PolynomialPDF::SetCoefficients ( const std::vector< G4double > &  v)
inline

Definition at line 59 of file G4PolynomialPDF.hh.

59  {
60  fCoefficients = v; fChanged = true; Simplify();
61  }
std::vector< G4double > fCoefficients
tuple v
Definition: test.py:18

Here is the call graph for this function:

Here is the caller graph for this function:

void G4PolynomialPDF::SetCoefficients ( size_t  n,
const G4double coeffs 
)

Definition at line 63 of file G4PolynomialPDF.cc.

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 }
void SetNCoefficients(size_t n)
size_t GetNCoefficients() const
void SetCoefficient(size_t i, G4double value, bool doSimplify)

Here is the call graph for this function:

void G4PolynomialPDF::SetDomain ( G4double  x1,
G4double  x2 
)

Definition at line 86 of file G4PolynomialPDF.cc.

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 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
void G4PolynomialPDF::SetNCoefficients ( size_t  n)
inline

Definition at line 57 of file G4PolynomialPDF.hh.

57 { fCoefficients.resize(n); fChanged = true; }
std::vector< G4double > fCoefficients
const G4int n

Here is the caller graph for this function:

void G4PolynomialPDF::SetTolerance ( G4double  tolerance)
inline

Definition at line 87 of file G4PolynomialPDF.hh.

87 { fTolerance = tolerance; }
void G4PolynomialPDF::Simplify ( )

Definition at line 74 of file G4PolynomialPDF.cc.

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 }
std::vector< G4double > fCoefficients
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

Here is the caller graph for this function:

Member Data Documentation

G4bool G4PolynomialPDF::fChanged
protected

Definition at line 114 of file G4PolynomialPDF.hh.

std::vector<G4double> G4PolynomialPDF::fCoefficients
protected

Definition at line 113 of file G4PolynomialPDF.hh.

G4double G4PolynomialPDF::fTolerance
protected

Definition at line 115 of file G4PolynomialPDF.hh.

G4int G4PolynomialPDF::fVerbose
protected

Definition at line 116 of file G4PolynomialPDF.hh.

G4double G4PolynomialPDF::fX1
protected

Definition at line 111 of file G4PolynomialPDF.hh.

G4double G4PolynomialPDF::fX2
protected

Definition at line 112 of file G4PolynomialPDF.hh.


The documentation for this class was generated from the following files: