Geant4 Cross Reference |
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 // Description: Evaluates, generates random numbers from, and evaluates 38 // the inverse of a polynomial PDF, its CDF, and its first and second 39 // derivative. 40 // 41 // ------------------------------------------------------------------- 42 43 #ifndef G4POLYNOMIALPDF_HH 44 #define G4POLYNOMIALPDF_HH 45 46 #include "globals.hh" 47 #include <vector> 48 49 class G4PolynomialPDF 50 { 51 public: 52 G4PolynomialPDF(size_t n = 0, const double* coeffs = nullptr, 53 G4double x1=0, G4double x2=1); 54 55 ~G4PolynomialPDF(); 56 // Setters and Getters for coefficients 57 inline void SetNCoefficients(size_t n) { fCoefficients.resize(n); fChanged = true; } 58 inline size_t GetNCoefficients() const { return fCoefficients.size(); } 59 inline void SetCoefficients(const std::vector<G4double>& v) { 60 fCoefficients = v; fChanged = true; Simplify(); 61 } 62 inline G4double GetCoefficient(size_t i) const { return fCoefficients[i]; } 63 void SetCoefficient(size_t i, G4double value, bool doSimplify); 64 void SetCoefficients(size_t n, const G4double* coeffs); 65 void Simplify(); 66 67 // Set the domain over which random numbers are generated and over which 68 // the CDF is evaluated 69 void SetDomain(G4double x1, G4double x2); 70 71 // Normalize PDF to 1 over domain fX1 to fX2. Used internally by 72 // GetRandomX(), but the user may want to call this as well for evaluation 73 // purposes. 74 void Normalize(); 75 76 // Evaluate (d/dx)^ddxPower f(x) (-1 <= ddxPower <= 2) 77 // ddxPower = -1 -> CDF; 78 // ddxPower = 0 -> PDF 79 // ddxPower = 1 -> PDF' 80 // ddxPower = 2 -> PDF'' 81 G4double Evaluate(G4double x, G4int ddxPower = 0); 82 83 // Generate a random number from this PDF 84 G4double GetRandomX(); 85 86 // Set the tolerance to within negative minima are checked 87 inline void SetTolerance(G4double tolerance) { fTolerance = tolerance; } 88 89 // Find a value x between x1 and x2 at which ddxPower[PDF](x) = p. 90 // ddxPower = -1 -> CDF; 91 // ddxPower = 0 -> PDF 92 // ddxPower = 1 -> PDF' 93 // (ddxPower = 2 not implemented) 94 // Solves analytically when possible, and otherwise uses the Newton-Raphson 95 // method to find the zero of ddxPower[PDF](x) - p. 96 // If not found in range, returns the nearest boundary. 97 // Beware that if x1 and x2 are not set carefully there may be multiple 98 // solutions, and care is not taken to select a particular one among them. 99 // Returns x2 on error 100 G4double GetX( G4double p, G4double x1, G4double x2, G4int ddxPower = 0, 101 G4double guess = 1.e99, G4bool bisect = true ); 102 inline G4double EvalInverseCDF(G4double p) { return GetX(p, fX1, fX2, -1, fX1 + p*(fX2-fX1)); } 103 G4double Bisect( G4double p, G4double x1, G4double x2 ); 104 105 void Dump(); 106 107 protected: 108 // Checks for negative values between x1 and x2. Used by GetRandomX() 109 G4bool HasNegativeMinimum(G4double x1, G4double x2); 110 111 G4double fX1; 112 G4double fX2; 113 std::vector<G4double> fCoefficients; 114 G4bool fChanged; 115 G4double fTolerance; 116 G4int fVerbose; 117 }; 118 119 #endif 120