Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/util/src/G4PolynomialPDF.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 
 43 G4PolynomialPDF::G4PolynomialPDF(size_t n, const G4double* coeffs, 
 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 
 51 G4PolynomialPDF::~G4PolynomialPDF() 
 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 
 74 void G4PolynomialPDF::Simplify()
 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 
 86 void G4PolynomialPDF::SetDomain(G4double x1, G4double x2) 
 87 { 
 88   if(x2 <= x1) {
 89     if(fVerbose > 0) {
 90       G4cout << "G4PolynomialPDF::SetDomain() WARNING: Invalid domain! "
 91        << "(x1 = " << x1 << ", x2 = " << x2 << ")." << G4endl;
 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 fX2.
103   /// Double-check that the highest-order coefficient is non-zero.
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 
131 G4double G4PolynomialPDF::Evaluate(G4double x, G4int ddxPower)
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() 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 
166 G4bool G4PolynomialPDF::HasNegativeMinimum(G4double x1, G4double x2)
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 
207 G4double G4PolynomialPDF::GetRandomX()
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 
223 G4double G4PolynomialPDF::GetX(G4double p, G4double x1, G4double x2, 
224              G4int ddxPower, G4double guess, G4bool bisect)
225 {
226   /// Find a value of X between x1 and x2 at which f(x) = p.
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 the zero of f(x) - p.
231   /// If not found in range, returns the nearest boundary
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: Biseting and trying again..." 
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, G4double x1, G4double x2 ) {
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 
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() - Interval: " << fX1 << " <= x < " 
412    << fX2 << G4endl;
413 }
414 
415 
416