| 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 * G4ErrorFunction.cc
28 *
29 * Created on: Jul 23, 2019
30 * Author: W. G. Shin
31 * J. Ramos-Mendez and B. Faddego
32 */
33
34 /*
35 Extracted from http://ab-initio.mit.edu/Fadde
36 Steven G. Johnson, October 2012.
37
38 Copyright © 2012 Massachusetts Institute of
39
40 Permission is hereby granted, free of charge,
41 obtaining a copy of this software and associa
42 files (the "Software"), to deal in the Softwa
43 including without limitation the rights to us
44 publish, distribute, sublicense, and/or sell
45 and to permit persons to whom the Software is
46 subject to the following conditions:
47
48 The above copyright notice and this permissio
49 included in all copies or substantial portion
50 */
51
52 #include "G4ErrorFunction.hh"
53 #include "globals.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "Randomize.hh"
56 #include <cmath>
57
58
59 G4ErrorFunction::G4ErrorFunction() = default;
60
61 G4ErrorFunction::~G4ErrorFunction() = default;
62
63
64
65 G4double G4ErrorFunction::erfcx_y100(G4double
66
67 {
68 switch ((int) y100) {
69 case 0: {
70 G4double t = 2*y100 - 1;
71 return 0.70878032454106438663e-3 +
72 }
73 case 1: {
74 G4double t = 2*y100 - 3;
75 return 0.21479143208285144230e-2 +
76 }
77 case 2: {
78 G4double t = 2*y100 - 5;
79 return 0.36165255935630175090e-2 +
80 }
81 case 3: {
82 G4double t = 2*y100 - 7;
83 return 0.51154983860031979264e-2 +
84 }
85 case 4: {
86 G4double t = 2*y100 - 9;
87 return 0.66457513172673049824e-2 +
88 }
89 case 5: {
90 G4double t = 2*y100 - 11;
91 return 0.82082389970241207883e-2 +
92 }
93 case 6: {
94 G4double t = 2*y100 - 13;
95 return 0.98039537275352193165e-2 +
96 }
97 case 7: {
98 G4double t = 2*y100 - 15;
99 return 0.11433927298290302370e-1 +
100 }
101 case 8: {
102 G4double t = 2*y100 - 17;
103 return 0.13099232878814653979e-1 +
104 }
105 case 9: {
106 G4double t = 2*y100 - 19;
107 return 0.14800987015587535621e-1 +
108 }
109 case 10: {
110 G4double t = 2*y100 - 21;
111 return 0.16540351739394069380e-1 +
112 }
113 case 11: {
114 G4double t = 2*y100 - 23;
115 return 0.18318536789842392647e-1 +
116 }
117 case 12: {
118 G4double t = 2*y100 - 25;
119 return 0.20136801964214276775e-1 +
120 }
121 case 13: {
122 G4double t = 2*y100 - 27;
123 return 0.21996459598282740954e-1 +
124 }
125 case 14: {
126 G4double t = 2*y100 - 29;
127 return 0.23898877187226319502e-1 +
128 }
129 case 15: {
130 G4double t = 2*y100 - 31;
131 return 0.25845480155298518485e-1 +
132 }
133 case 16: {
134 G4double t = 2*y100 - 33;
135 return 0.27837754783474696598e-1 +
136 }
137 case 17: {
138 G4double t = 2*y100 - 35;
139 return 0.29877251304899307550e-1 +
140 }
141 case 18: {
142 G4double t = 2*y100 - 37;
143 return 0.31965587178596443475e-1 +
144 }
145 case 19: {
146 G4double t = 2*y100 - 39;
147 return 0.34104450552588334840e-1 +
148 }
149 case 20: {
150 G4double t = 2*y100 - 41;
151 return 0.36295603928292425716e-1 +
152 }
153 case 21: {
154 G4double t = 2*y100 - 43;
155 return 0.38540888038840509795e-1 +
156 }
157 case 22: {
158 G4double t = 2*y100 - 45;
159 return 0.40842225954785960651e-1 +
160 }
161 case 23: {
162 G4double t = 2*y100 - 47;
163 return 0.43201627431540222422e-1 +
164 }
165 case 24: {
166 G4double t = 2*y100 - 49;
167 return 0.45621193513810471438e-1 +
168 }
169 case 25: {
170 G4double t = 2*y100 - 51;
171 return 0.48103121413299865517e-1 +
172 }
173 case 26: {
174 G4double t = 2*y100 - 53;
175 return 0.50649709676983338501e-1 +
176 }
177 case 27: {
178 G4double t = 2*y100 - 55;
179 return 0.53263363664388864181e-1 +
180 }
181 case 28: {
182 G4double t = 2*y100 - 57;
183 return 0.55946601353500013794e-1 +
184 }
185 case 29: {
186 G4double t = 2*y100 - 59;
187 return 0.58702059496154081813e-1 +
188 }
189 case 30: {
190 G4double t = 2*y100 - 61;
191 return 0.61532500145144778048e-1 +
192 }
193 case 31: {
194 G4double t = 2*y100 - 63;
195 return 0.64440817576653297993e-1 +
196 }
197 case 32: {
198 G4double t = 2*y100 - 65;
199 return 0.67430045633130393282e-1 +
200 }
201 case 33: {
202 G4double t = 2*y100 - 67;
203 return 0.70503365513338850709e-1 +
204 }
205 case 34: {
206 G4double t = 2*y100 - 69;
207 return 0.73664114037944596353e-1 +
208 }
209 case 35: {
210 G4double t = 2*y100 - 71;
211 return 0.76915792420819562379e-1 +
212 }
213 case 36: {
214 G4double t = 2*y100 - 73;
215 return 0.80262075578094612819e-1 +
216 }
217 case 37: {
218 G4double t = 2*y100 - 75;
219 return 0.83706822008980357446e-1 +
220 }
221 case 38: {
222 G4double t = 2*y100 - 77;
223 return 0.87254084284461718231e-1 +
224 }
225 case 39: {
226 G4double t = 2*y100 - 79;
227 return 0.90908120182172748487e-1 +
228 }
229 case 40: {
230 G4double t = 2*y100 - 81;
231 return 0.94673404508075481121e-1 +
232 }
233 case 41: {
234 G4double t = 2*y100 - 83;
235 return 0.98554641648004456555e-1 +
236 }
237 case 42: {
238 G4double t = 2*y100 - 85;
239 return 0.10255677889470089531e0 +
240 }
241 case 43: {
242 G4double t = 2*y100 - 87;
243 return 0.10668502059865093318e0 +
244 }
245 case 44: {
246 G4double t = 2*y100 - 89;
247 return 0.11094484319386444474e0 +
248 }
249 case 45: {
250 G4double t = 2*y100 - 91;
251 return 0.11534201115268804714e0 +
252 }
253 case 46: {
254 G4double t = 2*y100 - 93;
255 return 0.11988259392684094740e0 +
256 }
257 case 47: {
258 G4double t = 2*y100 - 95;
259 return 0.12457298393509812907e0 +
260 }
261 case 48: {
262 G4double t = 2*y100 - 97;
263 return 0.12941991566142438816e0 +
264 }
265 case 49: {
266 G4double t = 2*y100 - 99;
267 return 0.13443048593088696613e0 +
268 }
269 case 50: {
270 G4double t = 2*y100 - 101;
271 return 0.13961217543434561353e0 +
272 }
273 case 51: {
274 G4double t = 2*y100 - 103;
275 return 0.14497287157673800690e0 +
276 }
277 case 52: {
278 G4double t = 2*y100 - 105;
279 return 0.15052089272774618151e0 +
280 }
281 case 53: {
282 G4double t = 2*y100 - 107;
283 return 0.15626501395774612325e0 +
284 }
285 case 54: {
286 G4double t = 2*y100 - 109;
287 return 0.16221449434620737567e0 +
288 }
289 case 55: {
290 G4double t = 2*y100 - 111;
291 return 0.16837910595412130659e0 +
292 }
293 case 56: {
294 G4double t = 2*y100 - 113;
295 return 0.17476916455659369953e0 +
296 }
297 case 57: {
298 G4double t = 2*y100 - 115;
299 return 0.18139556223643701364e0 +
300 }
301 case 58: {
302 G4double t = 2*y100 - 117;
303 return 0.18826980194443664549e0 +
304 }
305 case 59: {
306 G4double t = 2*y100 - 119;
307 return 0.19540403413693967350e0 +
308 }
309 case 60: {
310 G4double t = 2*y100 - 121;
311 return 0.20281109560651886959e0 +
312 }
313 case 61: {
314 G4double t = 2*y100 - 123;
315 return 0.21050455062669334978e0 +
316 }
317 case 62: {
318 G4double t = 2*y100 - 125;
319 return 0.21849873453703332479e0 +
320 }
321 case 63: {
322 G4double t = 2*y100 - 127;
323 return 0.22680879990043229327e0 +
324 }
325 case 64: {
326 G4double t = 2*y100 - 129;
327 return 0.23545076536988703937e0 +
328 }
329 case 65: {
330 G4double t = 2*y100 - 131;
331 return 0.24444156740777432838e0 +
332 }
333 case 66: {
334 G4double t = 2*y100 - 133;
335 return 0.25379911500634264643e0 +
336 }
337 case 67: {
338 G4double t = 2*y100 - 135;
339 return 0.26354234756393613032e0 +
340 }
341 case 68: {
342 G4double t = 2*y100 - 137;
343 return 0.27369129607732343398e0 +
344 }
345 case 69: {
346 G4double t = 2*y100 - 139;
347 return 0.28426714781640316172e0 +
348 }
349 case 70: {
350 G4double t = 2*y100 - 141;
351 return 0.29529231465348519920e0 +
352 }
353 case 71: {
354 G4double t = 2*y100 - 143;
355 return 0.30679050522528838613e0 +
356 }
357 case 72: {
358 G4double t = 2*y100 - 145;
359 return 0.31878680111173319425e0 +
360 }
361 case 73: {
362 G4double t = 2*y100 - 147;
363 return 0.33130773722152622027e0 +
364 }
365 case 74: {
366 G4double t = 2*y100 - 149;
367 return 0.34438138658041336523e0 +
368 }
369 case 75: {
370 G4double t = 2*y100 - 151;
371 return 0.35803744972380175583e0 +
372 }
373 case 76: {
374 G4double t = 2*y100 - 153;
375 return 0.37230734890119724188e0 +
376 }
377 case 77: {
378 G4double t = 2*y100 - 155;
379 return 0.38722432730555448223e0 +
380 }
381 case 78: {
382 G4double t = 2*y100 - 157;
383 return 0.40282355354616940667e0 +
384 }
385 case 79: {
386 G4double t = 2*y100 - 159;
387 return 0.41914223158913787649e0 +
388 }
389 case 80: {
390 G4double t = 2*y100 - 161;
391 return 0.43621971639463786896e0 +
392 }
393 case 81: {
394 G4double t = 2*y100 - 163;
395 return 0.45409763548534330981e0 +
396 }
397 case 82: {
398 G4double t = 2*y100 - 165;
399 return 0.47282001668512331468e0 +
400 }
401 case 83: {
402 G4double t = 2*y100 - 167;
403 return 0.49243342227179841649e0 +
404 }
405 case 84: {
406 G4double t = 2*y100 - 169;
407 return 0.51298708979209258326e0 +
408 }
409 case 85: {
410 G4double t = 2*y100 - 171;
411 return 0.53453307979101369843e0 +
412 }
413 case 86: {
414 G4double t = 2*y100 - 173;
415 return 0.55712643071169299478e0 +
416 }
417 case 87: {
418 G4double t = 2*y100 - 175;
419 return 0.58082532122519320968e0 +
420 }
421 case 88: {
422 G4double t = 2*y100 - 177;
423 return 0.60569124025293375554e0 +
424 }
425 case 89: {
426 G4double t = 2*y100 - 179;
427 return 0.63178916494715716894e0 +
428 }
429 case 90: {
430 G4double t = 2*y100 - 181;
431 return 0.65918774689725319200e0 +
432 }
433 case 91: {
434 G4double t = 2*y100 - 183;
435 return 0.68795950683174433822e0 +
436 }
437 case 92: {
438 G4double t = 2*y100 - 185;
439 return 0.71818103808729967036e0 +
440 }
441 case 93: {
442 G4double t = 2*y100 - 187;
443 return 0.74993321911726254661e0 +
444 }
445 case 94: {
446 G4double t = 2*y100 - 189;
447 return 0.78330143531283492729e0 +
448 }
449 case 95: {
450 G4double t = 2*y100 - 191;
451 return 0.81837581041023811832e0 +
452 }
453 case 96: {
454 G4double t = 2*y100 - 193;
455 return 0.85525144775685126237e0 +
456 }
457 case 97: {
458 G4double t = 2*y100 - 195;
459 return 0.89402868170849933734e0 +
460 }
461 case 98: {
462 G4double t = 2*y100 - 197;
463 return 0.93481333942870796363e0 +
464 }
465 case 99: {
466 G4double t = 2*y100 - 199;
467 return 0.97771701335885035464e0 +
468 }
469 }
470 // we only get here if y = 1, i.e. |x| < 4
471 // erfcx is within 1e-15 of 1..
472 return 1.0;
473 }
474
475
476 G4double G4ErrorFunction::NormQuantile(G4doubl
477 {
478 G4double a0 = 3.3871328727963666080e0;
479 G4double a1 = 1.3314166789178437745e+2;
480 G4double a2 = 1.9715909503065514427e+3;
481 G4double a3 = 1.3731693765509461125e+4;
482 G4double a4 = 4.5921953931549871457e+4;
483 G4double a5 = 6.7265770927008700853e+4;
484 G4double a6 = 3.3430575583588128105e+4;
485 G4double a7 = 2.5090809287301226727e+3;
486 G4double b1 = 4.2313330701600911252e+1;
487 G4double b2 = 6.8718700749205790830e+2;
488 G4double b3 = 5.3941960214247511077e+3;
489 G4double b4 = 2.1213794301586595867e+4;
490 G4double b5 = 3.9307895800092710610e+4;
491 G4double b6 = 2.8729085735721942674e+4;
492 G4double b7 = 5.2264952788528545610e+3;
493 G4double c0 = 1.42343711074968357734e0;
494 G4double c1 = 4.63033784615654529590e0;
495 G4double c2 = 5.76949722146069140550e0;
496 G4double c3 = 3.64784832476320460504e0;
497 G4double c4 = 1.27045825245236838258e0;
498 G4double c5 = 2.41780725177450611770e-1;
499 G4double c6 = 2.27238449892691845833e-2;
500 G4double c7 = 7.74545014278341407640e-4;
501 G4double d1 = 2.05319162663775882187e0;
502 G4double d2 = 1.67638483018380384940e0;
503 G4double d3 = 6.89767334985100004550e-1;
504 G4double d4 = 1.48103976427480074590e-1;
505 G4double d5 = 1.51986665636164571966e-2;
506 G4double d6 = 5.47593808499534494600e-4;
507 G4double d7 = 1.05075007164441684324e-9;
508 G4double e0 = 6.65790464350110377720e0;
509 G4double e1 = 5.46378491116411436990e0;
510 G4double e2 = 1.78482653991729133580e0;
511 G4double e3 = 2.96560571828504891230e-1;
512 G4double e4 = 2.65321895265761230930e-2;
513 G4double e5 = 1.24266094738807843860e-3;
514 G4double e6 = 2.71155556874348757815e-5;
515 G4double e7 = 2.01033439929228813265e-7;
516 G4double f1 = 5.99832206555887937690e-1;
517 G4double f2 = 1.36929880922735805310e-1;
518 G4double f3 = 1.48753612908506148525e-2;
519 G4double f4 = 7.86869131145613259100e-4;
520 G4double f5 = 1.84631831751005468180e-5;
521 G4double f6 = 1.42151175831644588870e-7;
522 G4double f7 = 2.04426310338993978564e-15;
523
524 G4double split1 = 0.425;
525 G4double split2=5.;
526 G4double konst1=0.180625;
527 G4double konst2=1.6;
528
529 G4double q, r, quantile;
530 q=p-0.5;
531 if (std::abs(q)<split1) {
532 r=konst1-q*q;
533 quantile = q* (((((((a7 * r + a6) * r
534 * r + a2) * r + a1) *
535 (((((((b7 * r + b6) * r + b5) * r + b4
536 * r + b2) * r + b1) * r + 1.);
537 } else {
538 if(q<0) r=p;
539 else r=1-p;
540 //error case
541 if (r<=0)
542 quantile=0;
543 else {
544 r=std::sqrt(-std::log(r));
545 if (r<=split2) {
546 r=r-konst2;
547 quantile=(((((((c7 * r + c6) *
548 * r + c2) * r + c1
549 (((((((d7 * r + d6) * r + d5)
550 * r + d2) * r + d1) * r + 1
551 } else{
552 r=r-split2;
553 quantile=(((((((e7 * r + e6) *
554 * r + e2) * r + e1
555 (((((((f7 * r + f6) * r + f5)
556 * r + f2) * r + f1) * r + 1
557 }
558 if (q<0) quantile=-quantile;
559 }
560 }
561 return quantile;
562 }
563
564
565
566
567 G4double G4ErrorFunction::erfcx(G4double x)
568 {
569 if (x >= 0) {
570 if (x > 50) { // continued-fraction ex
571 const G4double ispi = 1./std::sqrt
572 if (x > 5e7) // 1-term expansion,
573 return ispi / x;
574 /* 5-term expansion (rely on compi
575 ispi / (x+0.5/(x+1/(x+1.5/(x+2/x)
576 return ispi*((x*x) * (x*x+4.5) + 2
577 }
578 return erfcx_y100(400/(4+x));
579 }
580
581 return x < -26.7 ? HUGE_VAL : (x < -6.1 ?
582 : 2*std::ex
583 }
584
585
586 G4double G4ErrorFunction::erfc(G4double x) {
587 return std::erfc(x);
588 }
589
590
591 G4double G4ErrorFunction::erfcWxy(G4double c,
592 return c * ( erfc(x) - std::exp(-x*x) * er
593 }
594
595
596 G4double G4ErrorFunction::Lambda(G4double x, G
597 return std::exp(-beta*beta/x) * ( 1.0 - al
598 }
599
600
601 G4double G4ErrorFunction::erfcInv(G4double x)
602 return - 0.70710678118654752440 * NormQuan
603 }
604
605
606