include/Geant4/G4PolynomialPDF.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // -------------------------------------------------------------------
0028 //      GEANT4 Class file
0029 //
0030 //
0031 //      File name:     G4PolynomialPDF
0032 //
0033 //      Author:        Jason Detwiler (jasondet@gmail.com)
0034 //
0035 //      Creation date: Aug 2012
0036 //
0037 //      Description:   Evaluates, generates random numbers from, and evaluates
0038 //      the inverse of a polynomial PDF, its CDF, and its first and second
0039 //      derivative.
0040 //
0041 // -------------------------------------------------------------------
0042
0043 #ifndef G4POLYNOMIALPDF_HH
0044 #define G4POLYNOMIALPDF_HH
0045
0046 #include "globals.hh"
0047 #include <vector>
0048
0049 class G4PolynomialPDF
0050 {
0051   public:
0052     G4PolynomialPDF(size_t n = 0, const double* coeffs = nullptr,
0053             G4double x1=0, G4double x2=1);
0054
0055     ~G4PolynomialPDF();
0056     // Setters and Getters for coefficients
0057     inline void SetNCoefficients(size_t n) { fCoefficients.resize(n); fChanged = true; }
0058     inline size_t GetNCoefficients() const { return fCoefficients.size(); }
0059     inline void SetCoefficients(const std::vector<G4double>& v) {
0060       fCoefficients = v; fChanged = true; Simplify();
0061     }
0062     inline G4double GetCoefficient(size_t i) const { return fCoefficients[i]; }
0063     void SetCoefficient(size_t i, G4double value, bool doSimplify);
0064     void SetCoefficients(size_t n, const G4double* coeffs);
0065     void Simplify();
0066
0067     // Set the domain over which random numbers are generated and over which
0068     // the CDF is evaluated
0069     void SetDomain(G4double x1, G4double x2);
0070
0071     // Normalize PDF to 1 over domain fX1 to fX2. Used internally by
0072     // GetRandomX(), but the user may want to call this as well for evaluation
0073     // purposes.
0074     void Normalize();
0075
0076     // Evaluate (d/dx)^ddxPower f(x) (-1 <= ddxPower <= 2)
0077     // ddxPower = -1 -> CDF;
0078     // ddxPower = 0 -> PDF
0079     // ddxPower = 1 -> PDF'
0080     // ddxPower = 2 -> PDF''
0081     G4double Evaluate(G4double x, G4int ddxPower = 0);
0082
0083     // Generate a random number from this PDF
0084     G4double GetRandomX();
0085
0086     // Set the tolerance to within negative minima are checked
0087     inline void SetTolerance(G4double tolerance) { fTolerance = tolerance; }
0088
0089     // Find a value x between x1 and x2 at which ddxPower[PDF](x) = p.
0090     // ddxPower = -1 -> CDF;
0091     // ddxPower = 0 -> PDF
0092     // ddxPower = 1 -> PDF'
0093     // (ddxPower = 2 not implemented)
0094     // Solves analytically when possible, and otherwise uses the Newton-Raphson
0095     // method to find the zero of ddxPower[PDF](x) - p.
0096     // If not found in range, returns the nearest boundary.
0097     // Beware that if x1 and x2 are not set carefully there may be multiple
0098     // solutions, and care is not taken to select a particular one among them.
0099     // Returns x2 on error
0100     G4double GetX( G4double p, G4double x1, G4double x2, G4int ddxPower = 0,
0101                    G4double guess = 1.e99, G4bool bisect = true );
0102     inline G4double EvalInverseCDF(G4double p) { return GetX(p, fX1, fX2, -1, fX1 + p*(fX2-fX1)); }
0103     G4double Bisect( G4double p, G4double x1, G4double x2 );
0104
0105     void Dump();
0106
0107   protected:
0108     // Checks for negative values between x1 and x2. Used by GetRandomX()
0109     G4bool HasNegativeMinimum(G4double x1, G4double x2);
0110
0111     G4double fX1;
0112     G4double fX2;
0113     std::vector<G4double> fCoefficients;
0114     G4bool   fChanged;
0115     G4double fTolerance;
0116     G4int    fVerbose;
0117 };
0118
0119 #endif