#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< G4double >	fCoefficients

G4bool	fChanged

G4double	fTolerance

G4int	fVerbose

Detailed Description

Definition at line 49 of file G4PolynomialPDF.hh.

Constructor & Destructor Documentation

◆ G4PolynomialPDF()

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

Definition at line 43 of file G4PolynomialPDF.cc.

                                           :
  fX1(x1), fX2(x2), fChanged(true), fTolerance(1.e-8), fVerbose(0)
{
  if(coeffs != nullptr) SetCoefficients(n, coeffs);
  else if(n > 0) SetNCoefficients(n);
}

◆ ~G4PolynomialPDF()

G4PolynomialPDF::~G4PolynomialPDF ( )

Definition at line 51 of file G4PolynomialPDF.cc.

52{}

Member Function Documentation

◆ Bisect()

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

Definition at line 392 of file G4PolynomialPDF.cc.

                                                                       {
  // Bisect to get 1% precision, then use Newton-Raphson
  G4double z = (x2 + x1)/2.0; // [x1 z x2]
  if((x2 - x1)/(fX2 - fX1) < 0.01) return GetX(p, fX1, fX2, -1, z, false);
  G4double fz = Evaluate(z, -1) - p;
  if(fz < 0) return Bisect(p, z, x2); // [z x2]
  return Bisect(p, x1, z); // [x1 z]
}

Referenced by Bisect(), and GetX().

◆ Dump()

void G4PolynomialPDF::Dump ( )

Definition at line 401 of file G4PolynomialPDF.cc.

{
  G4cout << "G4PolynomialPDF::Dump() - PDF(x) = ";
  for(size_t i=0; i<GetNCoefficients(); i++) {
    if(i>0) G4cout << " + ";
    G4cout << GetCoefficient(i);
    if(i>0) G4cout << "*x";
    if(i>1) G4cout << "^" << i;
  }
  G4cout << G4endl;
  G4cout << "G4PolynomialPDF::Dump() - Interval: " << fX1 << " <= x < " 
     << fX2 << G4endl;
}

Referenced by Normalize().

◆ EvalInverseCDF()

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)); }

Referenced by GetRandomX().

◆ Evaluate()

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.

{
  /// Evaluate f(x)
  /// ddxPower = -1: f = CDF
  /// ddxPower = 0: f = PDF
  /// ddxPower = 1: f = (d/dx) PDF
  /// ddxPower = 2: f = (d2/dx2) PDF
  if(ddxPower < -1 || ddxPower > 2) {
    if(fVerbose > 0) {
      G4cout << "G4PolynomialPDF::GetX() WARNING: ddxPower " << ddxPower 
         << " not implemented" << G4endl;
    }
    return 0.0;
  }
 
  double f = 0.; // return value
  double xN = 1.; // x to the power N
  double x1N = 1.; // endpoint x1 to the power N; only used by CDF
  for(size_t i=0; i<=GetNCoefficients(); ++i) {
    if(ddxPower == -1) { // CDF
      if(i>0) f += GetCoefficient(i-1)*(xN - x1N)/i;
      x1N *= fX1;
    }
    else if(ddxPower == 0 && i<GetNCoefficients()) f += GetCoefficient(i)*xN; // PDF
    else if(ddxPower == 1) { // (d/dx) PDF
      if(i<GetNCoefficients()-1) f += GetCoefficient(i+1)*xN*(i+1);
    }
    else if(ddxPower == 2) { // (d2/dx2) PDF
      if(i<GetNCoefficients()-2) f += GetCoefficient(i+2)*xN*((i+2)*(i+1));
    }
    xN *= x;
  }
  return f;
}

Referenced by Bisect(), and HasNegativeMinimum().

◆ GetCoefficient()

G4double G4PolynomialPDF::GetCoefficient ( size_t i ) const

inline

Definition at line 62 of file G4PolynomialPDF.hh.

62{ return fCoefficients[i]; }

G4PolynomialPDF::fCoefficients

std::vector< G4double > fCoefficients

Definition: G4PolynomialPDF.hh:113

Referenced by Dump(), Evaluate(), GetX(), HasNegativeMinimum(), and Normalize().

◆ GetNCoefficients()

size_t G4PolynomialPDF::GetNCoefficients ( ) const

inline

Definition at line 58 of file G4PolynomialPDF.hh.

58{ return fCoefficients.size(); }

Referenced by Dump(), Evaluate(), GetX(), HasNegativeMinimum(), Normalize(), and SetCoefficients().

◆ GetRandomX()

G4double G4PolynomialPDF::GetRandomX ( )

Definition at line 207 of file G4PolynomialPDF.cc.

{
  if(fChanged) {
    Normalize(); 
    if(HasNegativeMinimum(fX1, fX2)) {
      if(fVerbose > 0) {
    G4cout << "G4PolynomialPDF::GetRandomX() WARNING: PDF has negative values, returning 0..." 
           << G4endl;
      }
      return 0.0;
    }
    fChanged = false;
  }
  return EvalInverseCDF(G4UniformRand());
}

◆ GetX()

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.

{
  /// 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
 
  // input range checking
  if(GetNCoefficients() == 0) {
    if(fVerbose > 0) {
      G4cout << "G4PolynomialPDF::GetX() WARNING: no PDF defined!" << G4endl;
    }
    return x2;
  }
  if(ddxPower < -1 || ddxPower > 1) {
    if(fVerbose > 0) {
      G4cout << "G4PolynomialPDF::GetX() WARNING: ddxPower " << ddxPower 
         << " not implemented" << G4endl;
    }
    return x2;
  }
  if(ddxPower == -1 && (p<0 || p>1)) {
    if(fVerbose > 0) {
      G4cout << "G4PolynomialPDF::GetX() WARNING: p is out of range" << G4endl;
    }
    return fX2;
  }
 
  // check limits
  if(x2 <= x1 || x1 < fX1 || x2 > fX2) {
    if(fVerbose > 0) {
      G4cout << "G4PolynomialPDF::GetX() WARNING: domain must have fX1 <= x1 < x2 <= fX2. "
         << "You sent x1 = " << x1 << ", x2 = " << x2 << "." << G4endl;
    }
    return x2;
  }
 
  // Return x2 for flat lines
  if((ddxPower == 0 && GetNCoefficients() == 1) ||
     (ddxPower == 1 && GetNCoefficients() == 2)) return x2;
 
  // Solve p = mx + b -> x = (p-b)/m for linear functions
  if((ddxPower == -1 && GetNCoefficients() == 1) ||
     (ddxPower ==  0 && GetNCoefficients() == 2) ||
     (ddxPower ==  1 && GetNCoefficients() == 3)) {
    G4double b = (ddxPower > -1) ? GetCoefficient(ddxPower) : -GetCoefficient(0)*fX1;
    G4double slope = GetCoefficient(ddxPower+1); // the highest-order coefficient
    if(slope == 0) { // the highest-order coefficient should never be zero if simplified
      if(fVerbose > 0) {
        G4cout << "G4PolynomialPDF::GetX() WARNING: Got slope = 0. "
               << "Did you forget to Simplify()?" << G4endl;
      }
      return x2;
    }
    if(ddxPower == 1) slope *= 2.;
    G4double value = (p-b)/slope;
    if(value < x1) {
      return x1;
    }
    else if(value > x2) {
      return x2;
    }
    else {
      return value;
    }
  }
 
  // Solve quadratic equation for f-p=0 when f is quadratic
  if((ddxPower == -1 && GetNCoefficients() == 2) ||
     (ddxPower ==  0 && GetNCoefficients() == 3) ||
     (ddxPower ==  1 && GetNCoefficients() == 4)) {
    G4double c = -p + ((ddxPower > -1) ? GetCoefficient(ddxPower) : 0);
    if(ddxPower == -1) c -= (GetCoefficient(0) + GetCoefficient(1)/2.*fX1)*fX1;
    G4double b = GetCoefficient(ddxPower+1);
    if(ddxPower == 1) b *= 2.;
    G4double a = GetCoefficient(ddxPower+2); // the highest-order coefficient
    if(a == 0) { // the highest-order coefficient should never be 0 if simplified
      if(fVerbose > 0) {
        G4cout << "G4PolynomialPDF::GetX() WARNING: Got a = 0. "
               << "Did you forget to Simplify()?" << G4endl;
      }
      return x2;
    }
    if(ddxPower == 1) a *= 3;
    else if(ddxPower == -1) a *= 0.5;
    double sqrtFactor = b*b - 4.*a*c;
    if(sqrtFactor < 0) return x2; // quadratic equation has no solution (p not in range of f)
    sqrtFactor = sqrt(sqrtFactor)/2./fabs(a);
    G4double valueMinus = -b/2./a - sqrtFactor;
    if(valueMinus >= x1 && valueMinus <= x2) return valueMinus;
    else if(valueMinus > x2) return x2;
    G4double valuePlus = -b/2./a + sqrtFactor;
    if(valuePlus >= x1 && valuePlus <= x2) return valuePlus;
    else if(valuePlus < x1) return x2;
    return (x1-valueMinus <= valuePlus-x2) ? x1 : x2;
  }
 
  // f is non-trivial, so use Newton-Raphson
  // start in the middle if no good guess is provided
  if(guess < x1 || guess > x2) guess = (x2+x1)*0.5; 
  G4double lastChange = 1;
  size_t iterations = 0;
  while(fabs(lastChange) > fTolerance) {  /* Loop checking, 02.11.2015, A.Ribon */
    // calculate f and f' simultaneously
    G4double f = -p;
    G4double dfdx = 0;
    G4double xN = 1;
    G4double x1N = 1; // only used by CDF
    for(size_t i=0; i<=GetNCoefficients(); ++i) {
      if(ddxPower == -1) { // CDF
        if(i>0) f += GetCoefficient(i-1)*(xN - x1N)/G4double(i);
        if(i<GetNCoefficients()) dfdx += GetCoefficient(i)*xN;
        x1N *= fX1;
      }
      else if(ddxPower == 0) { // PDF
        if(i<GetNCoefficients()) f += GetCoefficient(i)*xN;
        if(i+1<GetNCoefficients()) dfdx += GetCoefficient(i+1)*xN*G4double(i+1);
      }
      else { // ddxPower == 1: (d/dx) PDF
        if(i+1<GetNCoefficients()) f += GetCoefficient(i+1)*xN*G4double(i+1);
        if(i+2<GetNCoefficients()) dfdx += GetCoefficient(i+2)*xN*G4double(i+2);
      }
      xN *= guess;
    }
    if(f == 0) return guess;
    if(dfdx == 0) {
      if(fVerbose > 0) {
    G4cout << "G4PolynomialPDF::GetX() WARNING: got f != 0 but slope = 0 for ddxPower = " 
           << ddxPower << G4endl;
      }
      return x2;
    }
    lastChange = - f/dfdx;
 
    if(guess + lastChange < x1) {
      lastChange = x1 - guess;
    } else if(guess + lastChange > x2) {
      lastChange = x2 - guess;
    }
 
    guess += lastChange;
    lastChange /= (fX2-fX1);
     
    ++iterations;
    if(iterations > 50) {
      if(p!=0) {
    if(fVerbose > 0) {
      G4cout << "G4PolynomialPDF::GetX() WARNING: got stuck searching for " << p 
         << " between " << x1 << " and " << x2 << " with ddxPower = " 
         << ddxPower 
         << ". Last guess was " << guess << "." << G4endl;
    }
      }
      if(ddxPower==-1 && bisect) {
    if(fVerbose > 0) {
      G4cout << "G4PolynomialPDF::GetX() WARNING: Biseting and trying again..." 
         << G4endl;
    }
        return Bisect(p, x1, x2);
      }
      else return guess;
    }
  }
  return guess;
}

Referenced by Bisect(), EvalInverseCDF(), and HasNegativeMinimum().

◆ HasNegativeMinimum()

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

protected

Definition at line 166 of file G4PolynomialPDF.cc.

{
  // ax2 + bx + c = 0
  // p': 2ax + b = 0 -> = 0 at min: x_extreme = -b/2a
 
  if(x1 < fX1 || x2 > fX2 || x2 < x1) {
    if(fVerbose > 0) {
      G4cout << "G4PolynomialPDF::HasNegativeMinimum() WARNING: Invalid range " 
         << x1 << " - " << x2 << G4endl;
    }
    return false;
  }
 
  // If flat, then check anywhere.
  if(GetNCoefficients() == 1) return (Evaluate(x1) < -fTolerance);
 
  // If linear, or if quadratic with negative second derivative, 
  // just check the endpoints
  if(GetNCoefficients() == 2 || 
     (GetNCoefficients() == 3 && GetCoefficient(2) <= 0)) {
    return (Evaluate(x1) < -fTolerance) || (Evaluate(x2) < -fTolerance);
  }
 
  // If quadratic and second dervative is positive, check at the mininum
  if(GetNCoefficients() == 3) {
    G4double xMin = -GetCoefficient(1)*0.5/GetCoefficient(2);
    if(xMin < x1) xMin = x1;
    if(xMin > x2) xMin = x2;
    return Evaluate(xMin) < -fTolerance;
  }
 
  // Higher-order polynomials: consider any extremum between x1 and x2. If none
  // are found, check the endpoints.
  G4double extremum = GetX(0, x1, x2, 1);
  if(Evaluate(extremum) < -fTolerance) return true;
  else if(extremum <= x1+(x2-x1)*fTolerance || 
      extremum >= x2-(x2-x1)*fTolerance) return false;
  else return 
     HasNegativeMinimum(x1, extremum) || HasNegativeMinimum(extremum, x2);
}

Referenced by GetRandomX(), and HasNegativeMinimum().

◆ Normalize()

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.

{
  /// Normalize PDF to 1 over domain fX1 to fX2.
  /// Double-check that the highest-order coefficient is non-zero.
  while(fCoefficients.size()) {  /* Loop checking, 30-Oct-2015, G.Folger */
    if(fCoefficients[fCoefficients.size()-1] == 0.0) fCoefficients.pop_back();
    else break;
  }
 
  G4double x1N = fX1, x2N = fX2;
  G4double sum = 0;
  for(size_t i=0; i<GetNCoefficients(); ++i) {
    sum += GetCoefficient(i)*(x2N - x1N)/G4double(i+1);
    x1N*=fX1;
    x2N*=fX2;
  }
  if(sum <= 0) {
    if(fVerbose > 0) {
      G4cout << "G4PolynomialPDF::Normalize() WARNING: PDF has non-positive area: " 
         << sum << G4endl;
      Dump();
    }
    return;
  }
 
  for(size_t i=0; i<GetNCoefficients(); ++i) {
    SetCoefficient(i, GetCoefficient(i)/sum, false);
  }
  Simplify();
}

Referenced by GetRandomX().

◆ SetCoefficient()

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

Definition at line 54 of file G4PolynomialPDF.cc.

{
  while(i >= fCoefficients.size()) fCoefficients.push_back(0);  
  /* Loop checking, 30-Oct-2015, G.Folger */
  fCoefficients[i] = value;
  fChanged = true;
  if(doSimplify) Simplify();
}

Referenced by Normalize(), and SetCoefficients().

◆ SetCoefficients() [1/2]

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

inline

Definition at line 59 of file G4PolynomialPDF.hh.

                                                              { 
      fCoefficients = v; fChanged = true; Simplify(); 
    }

Referenced by G4PolynomialPDF().

◆ SetCoefficients() [2/2]

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

Definition at line 63 of file G4PolynomialPDF.cc.

{
  SetNCoefficients(nCoeffs);
  for(size_t i=0; i<GetNCoefficients(); ++i) {
    SetCoefficient(i, coefficients[i], false);
  }
  fChanged = true;
  Simplify();
}

◆ SetDomain()

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

Definition at line 86 of file G4PolynomialPDF.cc.

{ 
  if(x2 <= x1) {
    if(fVerbose > 0) {
      G4cout << "G4PolynomialPDF::SetDomain() WARNING: Invalid domain! "
         << "(x1 = " << x1 << ", x2 = " << x2 << ")." << G4endl;
    }
    return;
  }
  fX1 = x1; 
  fX2 = x2; 
  fChanged = true; 
}

◆ SetNCoefficients()

void G4PolynomialPDF::SetNCoefficients ( size_t n )

inline

Definition at line 57 of file G4PolynomialPDF.hh.

57{ fCoefficients.resize(n); fChanged = true; }

Referenced by G4PolynomialPDF(), and SetCoefficients().

◆ SetTolerance()

void G4PolynomialPDF::SetTolerance ( G4double tolerance )

inline

Definition at line 87 of file G4PolynomialPDF.hh.

87{ fTolerance = tolerance; }

◆ Simplify()

void G4PolynomialPDF::Simplify ( )

Definition at line 74 of file G4PolynomialPDF.cc.

{
  while(fCoefficients.size() && fCoefficients[fCoefficients.size()-1] == 0) {
    if(fVerbose > 0) {
      G4cout << "G4PolynomialPDF::Simplify() WARNING: had to pop coefficient "
             << fCoefficients.size()-1 << G4endl;
    }
    fCoefficients.pop_back();
    fChanged = true;
  }
}

Referenced by Normalize(), SetCoefficient(), and SetCoefficients().

Member Data Documentation

◆ fChanged

G4bool G4PolynomialPDF::fChanged

protected

Definition at line 114 of file G4PolynomialPDF.hh.

Referenced by GetRandomX(), SetCoefficient(), SetCoefficients(), SetDomain(), SetNCoefficients(), and Simplify().

◆ fCoefficients

std::vector<G4double> G4PolynomialPDF::fCoefficients

protected

Definition at line 113 of file G4PolynomialPDF.hh.

Referenced by GetCoefficient(), GetNCoefficients(), Normalize(), SetCoefficient(), SetCoefficients(), SetNCoefficients(), and Simplify().

◆ fTolerance

G4double G4PolynomialPDF::fTolerance

protected

Definition at line 115 of file G4PolynomialPDF.hh.

Referenced by GetX(), HasNegativeMinimum(), and SetTolerance().

◆ fVerbose

G4int G4PolynomialPDF::fVerbose

protected

Definition at line 116 of file G4PolynomialPDF.hh.

Referenced by Evaluate(), GetRandomX(), GetX(), HasNegativeMinimum(), Normalize(), SetDomain(), and Simplify().

◆ fX1

G4double G4PolynomialPDF::fX1

protected

Definition at line 111 of file G4PolynomialPDF.hh.

Referenced by Bisect(), Dump(), EvalInverseCDF(), Evaluate(), GetRandomX(), GetX(), Normalize(), and SetDomain().

◆ fX2

G4double G4PolynomialPDF::fX2

protected

Definition at line 112 of file G4PolynomialPDF.hh.

Referenced by Bisect(), Dump(), EvalInverseCDF(), GetRandomX(), GetX(), HasNegativeMinimum(), Normalize(), and SetDomain().

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

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ G4PolynomialPDF()

◆ ~G4PolynomialPDF()

Member Function Documentation

◆ Bisect()

◆ Dump()

◆ EvalInverseCDF()

◆ Evaluate()

◆ GetCoefficient()

◆ GetNCoefficients()

◆ GetRandomX()

◆ GetX()

◆ HasNegativeMinimum()

◆ Normalize()

◆ SetCoefficient()

◆ SetCoefficients() [1/2]

◆ SetCoefficients() [2/2]

◆ SetDomain()

◆ SetNCoefficients()

◆ SetTolerance()

◆ Simplify()

Member Data Documentation

◆ fChanged

◆ fCoefficients

◆ fTolerance

◆ fVerbose

◆ fX1

◆ fX2