Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // ------------------------------------------- 28 // GEANT4 Class file 29 // 30 // 31 // File name: G4PolynomialPDF 32 // 33 // Author: Jason Detwiler (jasonde 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, con 44 G4double x1, G4double x2) : 45 fX1(x1), fX2(x2), fChanged(true), fTolerance 46 { 47 if(coeffs != nullptr) SetCoefficients(n, coe 48 else if(n > 0) SetNCoefficients(n); 49 } 50 51 G4PolynomialPDF::~G4PolynomialPDF() 52 {} 53 54 void G4PolynomialPDF::SetCoefficient(size_t i, 55 { 56 while(i >= fCoefficients.size()) fCoefficien 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 n 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[ 77 if(fVerbose > 0) { 78 G4cout << "G4PolynomialPDF::Simplify() W 79 << fCoefficients.size()-1 << G4en 80 } 81 fCoefficients.pop_back(); 82 fChanged = true; 83 } 84 } 85 86 void G4PolynomialPDF::SetDomain(G4double x1, G 87 { 88 if(x2 <= x1) { 89 if(fVerbose > 0) { 90 G4cout << "G4PolynomialPDF::SetDomain() 91 << "(x1 = " << x1 << ", x2 = " << x2 << 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 fX 103 /// Double-check that the highest-order coef 104 while(fCoefficients.size()) { /* Loop check 105 if(fCoefficients[fCoefficients.size()-1] = 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)/G4dou 113 x1N*=fX1; 114 x2N*=fX2; 115 } 116 if(sum <= 0) { 117 if(fVerbose > 0) { 118 G4cout << "G4PolynomialPDF::Normalize() 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, f 127 } 128 Simplify(); 129 } 130 131 G4double G4PolynomialPDF::Evaluate(G4double x, 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() WARNI 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 149 for(size_t i=0; i<=GetNCoefficients(); ++i) 150 if(ddxPower == -1) { // CDF 151 if(i>0) f += GetCoefficient(i-1)*(xN - x 152 x1N *= fX1; 153 } 154 else if(ddxPower == 0 && i<GetNCoefficient 155 else if(ddxPower == 1) { // (d/dx) PDF 156 if(i<GetNCoefficients()-1) f += GetCoeff 157 } 158 else if(ddxPower == 2) { // (d2/dx2) PDF 159 if(i<GetNCoefficients()-2) f += GetCoeff 160 } 161 xN *= x; 162 } 163 return f; 164 } 165 166 G4bool G4PolynomialPDF::HasNegativeMinimum(G4d 167 { 168 // ax2 + bx + c = 0 169 // p': 2ax + b = 0 -> = 0 at min: x_extreme 170 171 if(x1 < fX1 || x2 > fX2 || x2 < x1) { 172 if(fVerbose > 0) { 173 G4cout << "G4PolynomialPDF::HasNegativeM 174 << x1 << " - " << x2 << G4endl; 175 } 176 return false; 177 } 178 179 // If flat, then check anywhere. 180 if(GetNCoefficients() == 1) return (Evaluate 181 182 // If linear, or if quadratic with negative 183 // just check the endpoints 184 if(GetNCoefficients() == 2 || 185 (GetNCoefficients() == 3 && GetCoefficien 186 return (Evaluate(x1) < -fTolerance) || (Ev 187 } 188 189 // If quadratic and second dervative is posi 190 if(GetNCoefficients() == 3) { 191 G4double xMin = -GetCoefficient(1)*0.5/Get 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 ex 198 // are found, check the endpoints. 199 G4double extremum = GetX(0, x1, x2, 1); 200 if(Evaluate(extremum) < -fTolerance) return 201 else if(extremum <= x1+(x2-x1)*fTolerance || 202 extremum >= x2-(x2-x1)*fTolerance) return 203 else return 204 HasNegativeMinimum(x1, extremum) || HasNega 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() WAR 214 << G4endl; 215 } 216 return 0.0; 217 } 218 fChanged = false; 219 } 220 return EvalInverseCDF(G4UniformRand()); 221 } 222 223 G4double G4PolynomialPDF::GetX(G4double p, G4d 224 G4int ddxPower, G4double guess, G 225 { 226 /// Find a value of X between x1 and x2 at w 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 t 231 /// If not found in range, returns the neare 232 233 // input range checking 234 if(GetNCoefficients() == 0) { 235 if(fVerbose > 0) { 236 G4cout << "G4PolynomialPDF::GetX() WARNI 237 } 238 return x2; 239 } 240 if(ddxPower < -1 || ddxPower > 1) { 241 if(fVerbose > 0) { 242 G4cout << "G4PolynomialPDF::GetX() WARNI 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() WARNI 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() WARNI 258 << "You sent x1 = " << x1 << ", x2 = " 259 } 260 return x2; 261 } 262 263 // Return x2 for flat lines 264 if((ddxPower == 0 && GetNCoefficients() == 1 265 (ddxPower == 1 && GetNCoefficients() == 2 266 267 // Solve p = mx + b -> x = (p-b)/m for linea 268 if((ddxPower == -1 && GetNCoefficients() == 269 (ddxPower == 0 && GetNCoefficients() == 270 (ddxPower == 1 && GetNCoefficients() == 271 G4double b = (ddxPower > -1) ? GetCoeffici 272 G4double slope = GetCoefficient(ddxPower+1 273 if(slope == 0) { // the highest-order coef 274 if(fVerbose > 0) { 275 G4cout << "G4PolynomialPDF::GetX() WAR 276 << "Did you forget to Simplify( 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 294 if((ddxPower == -1 && GetNCoefficients() == 295 (ddxPower == 0 && GetNCoefficients() == 296 (ddxPower == 1 && GetNCoefficients() == 297 G4double c = -p + ((ddxPower > -1) ? GetCo 298 if(ddxPower == -1) c -= (GetCoefficient(0) 299 G4double b = GetCoefficient(ddxPower+1); 300 if(ddxPower == 1) b *= 2.; 301 G4double a = GetCoefficient(ddxPower+2); / 302 if(a == 0) { // the highest-order coeffici 303 if(fVerbose > 0) { 304 G4cout << "G4PolynomialPDF::GetX() WAR 305 << "Did you forget to Simplify( 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 313 sqrtFactor = sqrt(sqrtFactor)/2./fabs(a); 314 G4double valueMinus = -b/2./a - sqrtFactor 315 if(valueMinus >= x1 && valueMinus <= x2) r 316 else if(valueMinus > x2) return x2; 317 G4double valuePlus = -b/2./a + sqrtFactor; 318 if(valuePlus >= x1 && valuePlus <= x2) ret 319 else if(valuePlus < x1) return x2; 320 return (x1-valueMinus <= valuePlus-x2) ? x 321 } 322 323 // f is non-trivial, so use Newton-Raphson 324 // start in the middle if no good guess is p 325 if(guess < x1 || guess > x2) guess = (x2+x1) 326 G4double lastChange = 1; 327 size_t iterations = 0; 328 while(fabs(lastChange) > fTolerance) { /* L 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 - 337 if(i<GetNCoefficients()) dfdx += GetCo 338 x1N *= fX1; 339 } 340 else if(ddxPower == 0) { // PDF 341 if(i<GetNCoefficients()) f += GetCoeff 342 if(i+1<GetNCoefficients()) dfdx += Get 343 } 344 else { // ddxPower == 1: (d/dx) PDF 345 if(i+1<GetNCoefficients()) f += GetCoe 346 if(i+2<GetNCoefficients()) dfdx += Get 347 } 348 xN *= guess; 349 } 350 if(f == 0) return guess; 351 if(dfdx == 0) { 352 if(fVerbose > 0) { 353 G4cout << "G4PolynomialPDF::GetX() WARNING: 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 374 << " between " << x1 << " and " << x2 << 375 << ddxPower 376 << ". Last guess was " << guess << "." << 377 } 378 } 379 if(ddxPower==-1 && bisect) { 380 if(fVerbose > 0) { 381 G4cout << "G4PolynomialPDF::GetX() WARNING 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, 393 // Bisect to get 1% precision, then use Newt 394 G4double z = (x2 + x1)/2.0; // [x1 z x2] 395 if((x2 - x1)/(fX2 - fX1) < 0.01) return GetX 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() - Interva 412 << fX2 << G4endl; 413 } 414 415 416