Logo ROOT  

Reference Guide
VavilovFast.cxx
Go to the documentation of this file.
1// @(#)root/mathmore:$Id$
2// Authors: B. List 29.4.2010
3
4
5 /**********************************************************************
6 * *
7 * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
8 * *
9 * This library is free software; you can redistribute it and/or *
10 * modify it under the terms of the GNU General Public License *
11 * as published by the Free Software Foundation; either version 2 *
12 * of the License, or (at your option) any later version. *
13 * *
14 * This library is distributed in the hope that it will be useful, *
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
17 * General Public License for more details. *
18 * *
19 * You should have received a copy of the GNU General Public License *
20 * along with this library (see file COPYING); if not, write *
21 * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
22 * 330, Boston, MA 02111-1307 USA, or contact the author. *
23 * *
24 **********************************************************************/
25
26// Implementation file for class VavilovFast
27//
28// Created by: blist at Thu Apr 29 11:19:00 2010
29//
30// Last update: Thu Apr 29 11:19:00 2010
31//
32
33
34#include "Math/VavilovFast.h"
35#include "Math/PdfFuncMathCore.h"
36#include "Math/ProbFuncMathCore.h"
37#include "Math/SpecFuncMathCore.h"
38#include "Math/SpecFuncMathMore.h"
39
40#include <iostream>
41#include <cmath>
42#include <limits>
43
44
45namespace ROOT {
46namespace Math {
47
48VavilovFast *VavilovFast::fgInstance = 0;
49
50
51VavilovFast::VavilovFast(double kappa, double beta2)
52{
53 SetKappaBeta2 (kappa, beta2);
54}
55
56
57VavilovFast::~VavilovFast()
58{
59 // desctructor (clean up resources)
60}
61
62void VavilovFast::SetKappaBeta2 (double kappa, double beta2)
63{
64 // Modified version of void TMath::VavilovSet(Double_t rkappa, Double_t beta2, Bool_t mode, Double_t *WCM, Double_t *AC, Double_t *HC, Int_t &itype, Int_t &npt)
65 fKappa = kappa;
66 fBeta2 = beta2;
67
68 double BKMNX1 = 0.02, BKMNY1 = 0.05, BKMNX2 = 0.12, BKMNY2 = 0.05,
69 BKMNX3 = 0.22, BKMNY3 = 0.05, BKMXX1 = 0.1 , BKMXY1 = 1,
70 BKMXX2 = 0.2 , BKMXY2 = 1 , BKMXX3 = 0.3 , BKMXY3 = 1;
71
72 double FBKX1 = 2/(BKMXX1-BKMNX1), FBKX2 = 2/(BKMXX2-BKMNX2),
73 FBKX3 = 2/(BKMXX3-BKMNX3), FBKY1 = 2/(BKMXY1-BKMNY1),
74 FBKY2 = 2/(BKMXY2-BKMNY2), FBKY3 = 2/(BKMXY3-BKMNY3);
75
76 double FNINV[] = {0, 1, 0.5, 0.33333333, 0.25, 0.2};
77
78 double EDGEC[]= {0, 0, 0.16666667e+0, 0.41666667e-1, 0.83333333e-2,
79 0.13888889e-1, 0.69444444e-2, 0.77160493e-3};
80
81 double U1[] = {0, 0.25850868e+0, 0.32477982e-1, -0.59020496e-2,
82 0. , 0.24880692e-1, 0.47404356e-2,
83 -0.74445130e-3, 0.73225731e-2, 0. ,
84 0.11668284e-2, 0. , -0.15727318e-2,-0.11210142e-2};
85
86 double U2[] = {0, 0.43142611e+0, 0.40797543e-1, -0.91490215e-2,
87 0. , 0.42127077e-1, 0.73167928e-2,
88 -0.14026047e-2, 0.16195241e-1, 0.24714789e-2,
89 0.20751278e-2, 0. , -0.25141668e-2,-0.14064022e-2};
90
91 double U3[] = {0, 0.25225955e+0, 0.64820468e-1, -0.23615759e-1,
92 0. , 0.23834176e-1, 0.21624675e-2,
93 -0.26865597e-2, -0.54891384e-2, 0.39800522e-2,
94 0.48447456e-2, -0.89439554e-2, -0.62756944e-2,-0.24655436e-2};
95
96 double U4[] = {0, 0.12593231e+1, -0.20374501e+0, 0.95055662e-1,
97 -0.20771531e-1, -0.46865180e-1, -0.77222986e-2,
98 0.32241039e-2, 0.89882920e-2, -0.67167236e-2,
99 -0.13049241e-1, 0.18786468e-1, 0.14484097e-1};
100
101 double U5[] = {0, -0.24864376e-1, -0.10368495e-2, 0.14330117e-2,
102 0.20052730e-3, 0.18751903e-2, 0.12668869e-2,
103 0.48736023e-3, 0.34850854e-2, 0. ,
104 -0.36597173e-3, 0.19372124e-2, 0.70761825e-3, 0.46898375e-3};
105
106 double U6[] = {0, 0.35855696e-1, -0.27542114e-1, 0.12631023e-1,
107 -0.30188807e-2, -0.84479939e-3, 0. ,
108 0.45675843e-3, -0.69836141e-2, 0.39876546e-2,
109 -0.36055679e-2, 0. , 0.15298434e-2, 0.19247256e-2};
110
111 double U7[] = {0, 0.10234691e+2, -0.35619655e+1, 0.69387764e+0,
112 -0.14047599e+0, -0.19952390e+1, -0.45679694e+0,
113 0. , 0.50505298e+0};
114 double U8[] = {0, 0.21487518e+2, -0.11825253e+2, 0.43133087e+1,
115 -0.14500543e+1, -0.34343169e+1, -0.11063164e+1,
116 -0.21000819e+0, 0.17891643e+1, -0.89601916e+0,
117 0.39120793e+0, 0.73410606e+0, 0. ,-0.32454506e+0};
118
119 double V1[] = {0, 0.27827257e+0, -0.14227603e-2, 0.24848327e-2,
120 0. , 0.45091424e-1, 0.80559636e-2,
121 -0.38974523e-2, 0. , -0.30634124e-2,
122 0.75633702e-3, 0.54730726e-2, 0.19792507e-2};
123
124 double V2[] = {0, 0.41421789e+0, -0.30061649e-1, 0.52249697e-2,
125 0. , 0.12693873e+0, 0.22999801e-1,
126 -0.86792801e-2, 0.31875584e-1, -0.61757928e-2,
127 0. , 0.19716857e-1, 0.32596742e-2};
128
129 double V3[] = {0, 0.20191056e+0, -0.46831422e-1, 0.96777473e-2,
130 -0.17995317e-2, 0.53921588e-1, 0.35068740e-2,
131 -0.12621494e-1, -0.54996531e-2, -0.90029985e-2,
132 0.34958743e-2, 0.18513506e-1, 0.68332334e-2,-0.12940502e-2};
133
134 double V4[] = {0, 0.13206081e+1, 0.10036618e+0, -0.22015201e-1,
135 0.61667091e-2, -0.14986093e+0, -0.12720568e-1,
136 0.24972042e-1, -0.97751962e-2, 0.26087455e-1,
137 -0.11399062e-1, -0.48282515e-1, -0.98552378e-2};
138
139 double V5[] = {0, 0.16435243e-1, 0.36051400e-1, 0.23036520e-2,
140 -0.61666343e-3, -0.10775802e-1, 0.51476061e-2,
141 0.56856517e-2, -0.13438433e-1, 0. ,
142 0. , -0.25421507e-2, 0.20169108e-2,-0.15144931e-2};
143
144 double V6[] = {0, 0.33432405e-1, 0.60583916e-2, -0.23381379e-2,
145 0.83846081e-3, -0.13346861e-1, -0.17402116e-2,
146 0.21052496e-2, 0.15528195e-2, 0.21900670e-2,
147 -0.13202847e-2, -0.45124157e-2, -0.15629454e-2, 0.22499176e-3};
148
149 double V7[] = {0, 0.54529572e+1, -0.90906096e+0, 0.86122438e-1,
150 0. , -0.12218009e+1, -0.32324120e+0,
151 -0.27373591e-1, 0.12173464e+0, 0. ,
152 0. , 0.40917471e-1};
153
154 double V8[] = {0, 0.93841352e+1, -0.16276904e+1, 0.16571423e+0,
155 0. , -0.18160479e+1, -0.50919193e+0,
156 -0.51384654e-1, 0.21413992e+0, 0. ,
157 0. , 0.66596366e-1};
158
159 double W1[] = {0, 0.29712951e+0, 0.97572934e-2, 0. ,
160 -0.15291686e-2, 0.35707399e-1, 0.96221631e-2,
161 -0.18402821e-2, -0.49821585e-2, 0.18831112e-2,
162 0.43541673e-2, 0.20301312e-2, -0.18723311e-2,-0.73403108e-3};
163
164 double W2[] = {0, 0.40882635e+0, 0.14474912e-1, 0.25023704e-2,
165 -0.37707379e-2, 0.18719727e+0, 0.56954987e-1,
166 0. , 0.23020158e-1, 0.50574313e-2,
167 0.94550140e-2, 0.19300232e-1};
168
169 double W3[] = {0, 0.16861629e+0, 0. , 0.36317285e-2,
170 -0.43657818e-2, 0.30144338e-1, 0.13891826e-1,
171 -0.58030495e-2, -0.38717547e-2, 0.85359607e-2,
172 0.14507659e-1, 0.82387775e-2, -0.10116105e-1,-0.55135670e-2};
173
174 double W4[] = {0, 0.13493891e+1, -0.26863185e-2, -0.35216040e-2,
175 0.24434909e-1, -0.83447911e-1, -0.48061360e-1,
176 0.76473951e-2, 0.24494430e-1, -0.16209200e-1,
177 -0.37768479e-1, -0.47890063e-1, 0.17778596e-1, 0.13179324e-1};
178
179 double W5[] = {0, 0.10264945e+0, 0.32738857e-1, 0. ,
180 0.43608779e-2, -0.43097757e-1, -0.22647176e-2,
181 0.94531290e-2, -0.12442571e-1, -0.32283517e-2,
182 -0.75640352e-2, -0.88293329e-2, 0.52537299e-2, 0.13340546e-2};
183
184 double W6[] = {0, 0.29568177e-1, -0.16300060e-2, -0.21119745e-3,
185 0.23599053e-2, -0.48515387e-2, -0.40797531e-2,
186 0.40403265e-3, 0.18200105e-2, -0.14346306e-2,
187 -0.39165276e-2, -0.37432073e-2, 0.19950380e-2, 0.12222675e-2};
188
189 double W8[] = {0, 0.66184645e+1, -0.73866379e+0, 0.44693973e-1,
190 0. , -0.14540925e+1, -0.39529833e+0,
191 -0.44293243e-1, 0.88741049e-1};
192
193 fItype = 0;
194 if (fKappa <0.01 || fKappa >12) {
195 std::cerr << "VavilovFast::set: illegal value of kappa=" << kappa << std::endl;
196 if (fKappa < 0.01) fKappa = 0.01;
197 else if (fKappa > 12) fKappa = 12;
198 }
199
200 double DRK[6];
201 double DSIGM[6];
202 double ALFA[8];
203 int j;
204 double x, y, xx, yy, x2, x3, y2, y3, xy, p2, p3, q2, q3, pq;
205 if (fKappa >= 0.29) {
206 fItype = 1;
207 fNpt = 100;
208 double wk = 1./std::sqrt(fKappa);
209
210 fAC[0] = (-0.032227*fBeta2-0.074275)*fKappa + (0.24533*fBeta2+0.070152)*wk + (-0.55610*fBeta2-3.1579);
211 fAC[8] = (-0.013483*fBeta2-0.048801)*fKappa + (-1.6921*fBeta2+8.3656)*wk + (-0.73275*fBeta2-3.5226);
212 DRK[1] = wk*wk;
213 DSIGM[1] = std::sqrt(fKappa/(1-0.5*fBeta2));
214 for (j=1; j<=4; j++) {
215 DRK[j+1] = DRK[1]*DRK[j];
216 DSIGM[j+1] = DSIGM[1]*DSIGM[j];
217 ALFA[j+1] = (FNINV[j]-fBeta2*FNINV[j+1])*DRK[j];
218 }
219 fHC[0]=std::log(fKappa)+fBeta2+0.42278434;
220 fHC[1]=DSIGM[1];
221 fHC[2]=ALFA[3]*DSIGM[3];
222 fHC[3]=(3*ALFA[2]*ALFA[2] + ALFA[4])*DSIGM[4]-3;
223 fHC[4]=(10*ALFA[2]*ALFA[3]+ALFA[5])*DSIGM[5]-10*fHC[2];
224 fHC[5]=fHC[2]*fHC[2];
225 fHC[6]=fHC[2]*fHC[3];
226 fHC[7]=fHC[2]*fHC[5];
227 for (j=2; j<=7; j++)
228 fHC[j]*=EDGEC[j];
229 fHC[8]=0.39894228*fHC[1];
230 }
231 else if (fKappa >=0.22) {
232 fItype = 2;
233 fNpt = 150;
234 x = 1+(fKappa-BKMXX3)*FBKX3;
235 y = 1+(std::sqrt(fBeta2)-BKMXY3)*FBKY3;
236 xx = 2*x;
237 yy = 2*y;
238 x2 = xx*x-1;
239 x3 = xx*x2-x;
240 y2 = yy*y-1;
241 y3 = yy*y2-y;
242 xy = x*y;
243 p2 = x2*y;
244 p3 = x3*y;
245 q2 = y2*x;
246 q3 = y3*x;
247 pq = x2*y2;
248 fAC[1] = W1[1] + W1[2]*x + W1[4]*x3 + W1[5]*y + W1[6]*y2 + W1[7]*y3 +
249 W1[8]*xy + W1[9]*p2 + W1[10]*p3 + W1[11]*q2 + W1[12]*q3 + W1[13]*pq;
250 fAC[2] = W2[1] + W2[2]*x + W2[3]*x2 + W2[4]*x3 + W2[5]*y + W2[6]*y2 +
251 W2[8]*xy + W2[9]*p2 + W2[10]*p3 + W2[11]*q2;
252 fAC[3] = W3[1] + W3[3]*x2 + W3[4]*x3 + W3[5]*y + W3[6]*y2 + W3[7]*y3 +
253 W3[8]*xy + W3[9]*p2 + W3[10]*p3 + W3[11]*q2 + W3[12]*q3 + W3[13]*pq;
254 fAC[4] = W4[1] + W4[2]*x + W4[3]*x2 + W4[4]*x3 + W4[5]*y + W4[6]*y2 + W4[7]*y3 +
255 W4[8]*xy + W4[9]*p2 + W4[10]*p3 + W4[11]*q2 + W4[12]*q3 + W4[13]*pq;
256 fAC[5] = W5[1] + W5[2]*x + W5[4]*x3 + W5[5]*y + W5[6]*y2 + W5[7]*y3 +
257 W5[8]*xy + W5[9]*p2 + W5[10]*p3 + W5[11]*q2 + W5[12]*q3 + W5[13]*pq;
258 fAC[6] = W6[1] + W6[2]*x + W6[3]*x2 + W6[4]*x3 + W6[5]*y + W6[6]*y2 + W6[7]*y3 +
259 W6[8]*xy + W6[9]*p2 + W6[10]*p3 + W6[11]*q2 + W6[12]*q3 + W6[13]*pq;
260 fAC[8] = W8[1] + W8[2]*x + W8[3]*x2 + W8[5]*y + W8[6]*y2 + W8[7]*y3 + W8[8]*xy;
261 fAC[0] = -3.05;
262 } else if (fKappa >= 0.12) {
263 fItype = 3;
264 fNpt = 200;
265 x = 1 + (fKappa-BKMXX2)*FBKX2;
266 y = 1 + (std::sqrt(fBeta2)-BKMXY2)*FBKY2;
267 xx = 2*x;
268 yy = 2*y;
269 x2 = xx*x-1;
270 x3 = xx*x2-x;
271 y2 = yy*y-1;
272 y3 = yy*y2-y;
273 xy = x*y;
274 p2 = x2*y;
275 p3 = x3*y;
276 q2 = y2*x;
277 q3 = y3*x;
278 pq = x2*y2;
279 fAC[1] = V1[1] + V1[2]*x + V1[3]*x2 + V1[5]*y + V1[6]*y2 + V1[7]*y3 +
280 V1[9]*p2 + V1[10]*p3 + V1[11]*q2 + V1[12]*q3;
281 fAC[2] = V2[1] + V2[2]*x + V2[3]*x2 + V2[5]*y + V2[6]*y2 + V2[7]*y3 +
282 V2[8]*xy + V2[9]*p2 + V2[11]*q2 + V2[12]*q3;
283 fAC[3] = V3[1] + V3[2]*x + V3[3]*x2 + V3[4]*x3 + V3[5]*y + V3[6]*y2 + V3[7]*y3 +
284 V3[8]*xy + V3[9]*p2 + V3[10]*p3 + V3[11]*q2 + V3[12]*q3 + V3[13]*pq;
285 fAC[4] = V4[1] + V4[2]*x + V4[3]*x2 + V4[4]*x3 + V4[5]*y + V4[6]*y2 + V4[7]*y3 +
286 V4[8]*xy + V4[9]*p2 + V4[10]*p3 + V4[11]*q2 + V4[12]*q3;
287 fAC[5] = V5[1] + V5[2]*x + V5[3]*x2 + V5[4]*x3 + V5[5]*y + V5[6]*y2 + V5[7]*y3 +
288 V5[8]*xy + V5[11]*q2 + V5[12]*q3 + V5[13]*pq;
289 fAC[6] = V6[1] + V6[2]*x + V6[3]*x2 + V6[4]*x3 + V6[5]*y + V6[6]*y2 + V6[7]*y3 +
290 V6[8]*xy + V6[9]*p2 + V6[10]*p3 + V6[11]*q2 + V6[12]*q3 + V6[13]*pq;
291 fAC[7] = V7[1] + V7[2]*x + V7[3]*x2 + V7[5]*y + V7[6]*y2 + V7[7]*y3 +
292 V7[8]*xy + V7[11]*q2;
293 fAC[8] = V8[1] + V8[2]*x + V8[3]*x2 + V8[5]*y + V8[6]*y2 + V8[7]*y3 +
294 V8[8]*xy + V8[11]*q2;
295 fAC[0] = -3.04;
296 } else {
297 fItype = 4;
298 if (fKappa >=0.02) fItype = 3;
299 fNpt = 200;
300 x = 1+(fKappa-BKMXX1)*FBKX1;
301 y = 1+(std::sqrt(fBeta2)-BKMXY1)*FBKY1;
302 xx = 2*x;
303 yy = 2*y;
304 x2 = xx*x-1;
305 x3 = xx*x2-x;
306 y2 = yy*y-1;
307 y3 = yy*y2-y;
308 xy = x*y;
309 p2 = x2*y;
310 p3 = x3*y;
311 q2 = y2*x;
312 q3 = y3*x;
313 pq = x2*y2;
314 if (fItype==3){
315 fAC[1] = U1[1] + U1[2]*x + U1[3]*x2 + U1[5]*y + U1[6]*y2 + U1[7]*y3 +
316 U1[8]*xy + U1[10]*p3 + U1[12]*q3 + U1[13]*pq;
317 fAC[2] = U2[1] + U2[2]*x + U2[3]*x2 + U2[5]*y + U2[6]*y2 + U2[7]*y3 +
318 U2[8]*xy + U2[9]*p2 + U2[10]*p3 + U2[12]*q3 + U2[13]*pq;
319 fAC[3] = U3[1] + U3[2]*x + U3[3]*x2 + U3[5]*y + U3[6]*y2 + U3[7]*y3 +
320 U3[8]*xy + U3[9]*p2 + U3[10]*p3 + U3[11]*q2 + U3[12]*q3 + U3[13]*pq;
321 fAC[4] = U4[1] + U4[2]*x + U4[3]*x2 + U4[4]*x3 + U4[5]*y + U4[6]*y2 + U4[7]*y3 +
322 U4[8]*xy + U4[9]*p2 + U4[10]*p3 + U4[11]*q2 + U4[12]*q3;
323 fAC[5] = U5[1] + U5[2]*x + U5[3]*x2 + U5[4]*x3 + U5[5]*y + U5[6]*y2 + U5[7]*y3 +
324 U5[8]*xy + U5[10]*p3 + U5[11]*q2 + U5[12]*q3 + U5[13]*pq;
325 fAC[6] = U6[1] + U6[2]*x + U6[3]*x2 + U6[4]*x3 + U6[5]*y + U6[7]*y3 +
326 U6[8]*xy + U6[9]*p2 + U6[10]*p3 + U6[12]*q3 + U6[13]*pq;
327 fAC[7] = U7[1] + U7[2]*x + U7[3]*x2 + U7[4]*x3 + U7[5]*y + U7[6]*y2 + U7[8]*xy;
328 }
329 fAC[8] = U8[1] + U8[2]*x + U8[3]*x2 + U8[4]*x3 + U8[5]*y + U8[6]*y2 + U8[7]*y3 +
330 U8[8]*xy + U8[9]*p2 + U8[10]*p3 + U8[11]*q2 + U8[13]*pq;
331 fAC[0] = -3.03;
332 }
333
334 fAC[9] = (fAC[8] - fAC[0])/fNpt;
335 fAC[10] = 1./fAC[9];
336 if (fItype == 3) {
337 x = (fAC[7]-fAC[8])/(fAC[7]*fAC[8]);
338 y = 1./std::log (fAC[8]/fAC[7]);
339 p2 = fAC[7]*fAC[7];
340 fAC[11] = p2*(fAC[1]*std::exp(-fAC[2]*(fAC[7]+fAC[5]*p2)-
341 fAC[3]*std::exp(-fAC[4]*(fAC[7]+fAC[6]*p2)))-0.045*y/fAC[7])/(1+x*y*fAC[7]);
342 fAC[12] = (0.045+x*fAC[11])*y;
343 }
344 if (fItype == 4) fAC[13] = 0.995/ROOT::Math::landau_cdf(fAC[8]);
345
346 //
347 x = fAC[0];
348 fWCM[0] = 0;
349 double fl, fu;
350 int k;
351 fl = Pdf (x);
352 for (k=1; k<=fNpt; k++) {
353 x += fAC[9];
354 fu = Pdf (x);
355 fWCM[k] = fWCM[k-1] + fl + fu;
356 fl = fu;
357 }
358 x = 0.5*fAC[9];
359 for (k=1; k<=fNpt; k++)
360 fWCM[k]*=x;
361}
362
363double VavilovFast::Pdf (double x) const
364{
365 // Modified version of TMath::double VavilovDenEval(Double_t rlam, Double_t *AC, Double_t *HC, Int_t itype);
366 //Internal function, called by Vavilov and VavilovSet
367
368 double v = 0;
369 if (x < fAC[0] || x > fAC[8])
370 return 0;
371 int k;
372 double h[10];
373 if (fItype ==1 ) {
374 double fn = 1;
375 double xx = (x + fHC[0])*fHC[1];
376 h[1] = xx;
377 h[2] = xx*xx -1;
378 for (k=2; k<=8; k++) {
379 fn++;
380 h[k+1] = xx*h[k]-fn*h[k-1];
381 }
382 double s = 1 + fHC[7]*h[9];
383 for (k=2; k<=6; k++)
384 s += fHC[k]*h[k+1];
385 if (s>0) v = fHC[8]*s*std::exp(-0.5*xx*xx);
386 }
387 else if (fItype == 2) {
388 double xx = x*x;
389 v = fAC[1]*std::exp(-fAC[2]*(x+fAC[5]*xx) - fAC[3]*std::exp(-fAC[4]*(x+fAC[6]*xx)));
390 }
391 else if (fItype == 3) {
392 if (x < fAC[7]) {
393 double xx = x*x;
394 v = fAC[1]*std::exp(-fAC[2]*(x+fAC[5]*xx)-fAC[3]*std::exp(-fAC[4]*(x+fAC[6]*xx)));
395 } else {
396 double xx = 1./x;
397 v = (fAC[11]*xx + fAC[12])*xx;
398 }
399 }
400 else if (fItype == 4) {
401 v = fAC[13]*ROOT::Math::landau_pdf(x);
402 }
403 return v;
404}
405
406
407double VavilovFast::Pdf (double x, double kappa, double beta2) {
408 //Returns the value of the Vavilov density function
409 //Parameters: 1st - the point were the density function is evaluated
410 // 2nd - value of kappa (distribution parameter)
411 // 3rd - value of beta2 (distribution parameter)
412 //The algorithm was taken from the CernLib function vavden(G115)
413 //Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
414 //Nucl.Instr. and Meth. B47(1990), 215-224
415 //Accuracy: quote from the reference above:
416 //"The resuls of our code have been compared with the values of the Vavilov
417 //density function computed numerically in an accurate way: our approximation
418 //shows a difference of less than 3% around the peak of the density function, slowly
419 //increasing going towards the extreme tails to the right and to the left"
420
421 if (kappa != fKappa || beta2 != fBeta2) SetKappaBeta2 (kappa, beta2);
422 return Pdf (x);
423}
424
425double VavilovFast::Cdf (double x) const {
426 // Modified version of Double_t TMath::VavilovI(Double_t x, Double_t kappa, Double_t beta2)
427 double xx, v;
428 if (x < fAC[0]) v = 0;
429 else if (x >= fAC[8]) v = 1;
430 else {
431 xx = x - fAC[0];
432 int k = int (xx*fAC[10]);
433 v = fWCM[k] + (xx - k*fAC[9])*(fWCM[k+1]-fWCM[k])*fAC[10];
434 if (v > 1) v = 1;
435 }
436 return v;
437}
438
439double VavilovFast::Cdf_c (double x) const {
440 return 1-Cdf(x);
441}
442
443double VavilovFast::Cdf (double x, double kappa, double beta2) {
444 //Returns the value of the Vavilov distribution function
445 //Parameters: 1st - the point were the density function is evaluated
446 // 2nd - value of kappa (distribution parameter)
447 // 3rd - value of beta2 (distribution parameter)
448 //The algorithm was taken from the CernLib function vavden(G115)
449 //Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
450 //Nucl.Instr. and Meth. B47(1990), 215-224
451 //Accuracy: quote from the reference above:
452 //"The resuls of our code have been compared with the values of the Vavilov
453 //density function computed numerically in an accurate way: our approximation
454 //shows a difference of less than 3% around the peak of the density function, slowly
455 //increasing going towards the extreme tails to the right and to the left"
456
457 if (kappa != fKappa || beta2 != fBeta2) SetKappaBeta2 (kappa, beta2);
458 return Cdf (x);
459}
460
461double VavilovFast::Cdf_c (double x, double kappa, double beta2) {
462 //Returns the value of the Vavilov distribution function
463 //Parameters: 1st - the point were the density function is evaluated
464 // 2nd - value of kappa (distribution parameter)
465 // 3rd - value of beta2 (distribution parameter)
466 //The algorithm was taken from the CernLib function vavden(G115)
467 //Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
468 //Nucl.Instr. and Meth. B47(1990), 215-224
469 //Accuracy: quote from the reference above:
470 //"The resuls of our code have been compared with the values of the Vavilov
471 //density function computed numerically in an accurate way: our approximation
472 //shows a difference of less than 3% around the peak of the density function, slowly
473 //increasing going towards the extreme tails to the right and to the left"
474
475 if (kappa != fKappa || beta2 != fBeta2) SetKappaBeta2 (kappa, beta2);
476 return Cdf_c (x);
477}
478
479double VavilovFast::Quantile (double z) const {
480 if (z < 0 || z > 1) return std::numeric_limits<double>::signaling_NaN();
481
482 // translated from CERNLIB routine VAVRAN by B. List 14.5.2010
483
484 double t = 2*z/fAC[9];
485 double rlam = fAC[0];
486 double fl = 0;
487 double fu = 0;
488 double s = 0;
489 double h[10];
490 for (int n = 1; n <= fNpt; ++n) {
491 rlam += fAC[9];
492 if (fItype == 1) {
493 double fn = 1;
494 double x = (rlam+fHC[0])*fHC[1];
495 h[1] = x;
496 h[2] = x*x-1;
497 for (int k = 2; k <= 8; ++k) {
498 ++fn;
499 h[k+1] = x*h[k]-fn*h[k-1];
500 }
501 double y = 1+fHC[7]*h[9];
502 for (int k = 2; k <= 6; ++k) {
503 y += fHC[k]*h[k+1];
504 }
505 if (y > 0) fu = fHC[8]*y*std::exp(-0.5*x*x);
506 }
507 else if (fItype == 2) {
508 double x = rlam*rlam;
509 fu = fAC[1]*std::exp(-fAC[2]*(rlam+fAC[5]*x)-
510 fAC[3]*std::exp(-fAC[4]*(rlam+fAC[6]*x)));
511 }
512 else if (fItype == 3) {
513 if (rlam < fAC[7]) {
514 double x = rlam*rlam;
515 fu = fAC[1]*std::exp(-fAC[2]*(rlam+fAC[5]*x)-
516 fAC[3]*std::exp(-fAC[4]*(rlam+fAC[6]*x)));
517 } else {
518 double x = 1/rlam;
519 fu = (fAC[11]*x+fAC[12])*x;
520 }
521 }
522 else {
523 fu = fAC[13]*Pdf(rlam); // in VAVRAN: AC(10) -> difference between VAVRAN and VAVSET
524 }
525 s += fl+fu;
526 if (s > t) break;
527 fl = fu;
528 }
529 double s0 = s-fl-fu;
530 double v = rlam-fAC[9];
531 if (s > s0) v += fAC[9]*(t-s0)/(s-s0);
532 return v;
533}
534
535double VavilovFast::Quantile (double z, double kappa, double beta2) {
536 if (kappa != fKappa || beta2 != fBeta2) SetKappaBeta2 (kappa, beta2);
537 return Quantile (z);
538}
539
540double VavilovFast::Quantile_c (double z) const {
541 if (z < 0 || z > 1) return std::numeric_limits<double>::signaling_NaN();
542 return Quantile (1-z);
543}
544
545double VavilovFast::Quantile_c (double z, double kappa, double beta2) {
546 if (kappa != fKappa || beta2 != fBeta2) SetKappaBeta2 (kappa, beta2);
547 return Quantile_c (z);
548}
549
550double VavilovFast::GetLambdaMin() const {
551 return fAC[0];
552}
553
554double VavilovFast::GetLambdaMax() const {
555 return fAC[8];
556}
557
558double VavilovFast::GetKappa() const {
559 return fKappa;
560}
561
562double VavilovFast::GetBeta2() const {
563 return fBeta2;
564}
565
566VavilovFast *VavilovFast::GetInstance() {
567 if (!fgInstance) fgInstance = new VavilovFast (1, 1);
568 return fgInstance;
569}
570
571VavilovFast *VavilovFast::GetInstance(double kappa, double beta2) {
572 if (!fgInstance) fgInstance = new VavilovFast (kappa, beta2);
573 else if (kappa != fgInstance->fKappa || beta2 != fgInstance->fBeta2) fgInstance->SetKappaBeta2 (kappa, beta2);
574 return fgInstance;
575}
576
577double vavilov_fast_pdf (double x, double kappa, double beta2) {
578 VavilovFast *vavilov = VavilovFast::GetInstance (kappa, beta2);
579 return vavilov->Pdf (x);
580}
581
582double vavilov_fast_cdf (double x, double kappa, double beta2) {
583 VavilovFast *vavilov = VavilovFast::GetInstance (kappa, beta2);
584 return vavilov->Cdf (x);
585}
586
587double vavilov_fast_cdf_c (double x, double kappa, double beta2) {
588 VavilovFast *vavilov = VavilovFast::GetInstance (kappa, beta2);
589 return vavilov->Cdf_c (x);
590}
591
592double vavilov_fast_quantile (double z, double kappa, double beta2) {
593 VavilovFast *vavilov = VavilovFast::GetInstance (kappa, beta2);
594 return vavilov->Quantile (z);
595}
596
597double vavilov_fast_quantile_c (double z, double kappa, double beta2) {
598 VavilovFast *vavilov = VavilovFast::GetInstance (kappa, beta2);
599 return vavilov->Quantile_c (z);
600}
601
602
603} // namespace Math
604} // namespace ROOT
s0
#define s0(x)
Definition: RSha256.hxx:90
h
#define h(i)
Definition: RSha256.hxx:106
x3
static const double x3[11]
Definition: RooGaussKronrodIntegrator1D.cxx:366
x2
Option_t Option_t TPoint TPoint const char x2
Definition: TGWin32VirtualXProxy.cxx:70
xy
Option_t Option_t TPoint xy
Definition: TGWin32VirtualXProxy.cxx:62
y2
Option_t Option_t TPoint TPoint const char y2
Definition: TGWin32VirtualXProxy.cxx:70
ROOT::Math::VavilovFast
Class describing a Vavilov distribution.
Definition: VavilovFast.h:116
ROOT::Math::VavilovFast::GetBeta2
double GetBeta2() const override
Return the current value of .
Definition: VavilovFast.cxx:562
ROOT::Math::VavilovFast::~VavilovFast
~VavilovFast() override
Destructor.
Definition: VavilovFast.cxx:57
ROOT::Math::VavilovFast::GetLambdaMax
double GetLambdaMax() const override
Return the maximum value of for which is nonzero in the current approximation.
Definition: VavilovFast.cxx:554
ROOT::Math::VavilovFast::fgInstance
static VavilovFast * fgInstance
Definition: VavilovFast.h:279
ROOT::Math::VavilovFast::GetKappa
double GetKappa() const override
Return the current value of .
Definition: VavilovFast.cxx:558
ROOT::Math::VavilovFast::Pdf
double Pdf(double x) const override
Evaluate the Vavilov probability density function.
Definition: VavilovFast.cxx:363
ROOT::Math::VavilovFast::SetKappaBeta2
void SetKappaBeta2(double kappa, double beta2) override
Change and and recalculate coefficients if necessary.
Definition: VavilovFast.cxx:62
ROOT::Math::VavilovFast::Quantile_c
double Quantile_c(double z) const override
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
Definition: VavilovFast.cxx:540
ROOT::Math::VavilovFast::Cdf_c
double Cdf_c(double x) const override
Evaluate the Vavilov complementary cumulative probability density function.
Definition: VavilovFast.cxx:439
ROOT::Math::VavilovFast::Quantile
double Quantile(double z) const override
Evaluate the inverse of the Vavilov cumulative probability density function.
Definition: VavilovFast.cxx:479
ROOT::Math::VavilovFast::Cdf
double Cdf(double x) const override
Evaluate the Vavilov cumulative probability density function.
Definition: VavilovFast.cxx:425
ROOT::Math::VavilovFast::GetLambdaMin
double GetLambdaMin() const override
Return the minimum value of for which is nonzero in the current approximation.
Definition: VavilovFast.cxx:550
ROOT::Math::VavilovFast::VavilovFast
VavilovFast(double kappa=1, double beta2=1)
Initialize an object to calculate the Vavilov distribution.
Definition: VavilovFast.cxx:51
ROOT::Math::VavilovFast::GetInstance
static VavilovFast * GetInstance()
Returns a static instance of class VavilovFast.
Definition: VavilovFast.cxx:566
ROOT::Math::vavilov_fast_pdf
double vavilov_fast_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
Definition: VavilovFast.cxx:577
ROOT::Math::landau_pdf
double landau_pdf(double x, double xi=1, double x0=0)
Probability density function of the Landau distribution:
Definition: PdfFuncMathCore.cxx:21
ROOT::Math::vavilov_fast_cdf
double vavilov_fast_cdf(double x, double kappa, double beta2)
The Vavilov cumulative probability density function.
Definition: VavilovFast.cxx:582
ROOT::Math::vavilov_fast_cdf_c
double vavilov_fast_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cumulative probability density function.
Definition: VavilovFast.cxx:587
ROOT::Math::landau_cdf
double landau_cdf(double x, double xi=1, double x0=0)
Cumulative distribution function of the Landau distribution (lower tail).
Definition: ProbFuncMathCore.cxx:336
ROOT::Math::vavilov_fast_quantile
double vavilov_fast_quantile(double z, double kappa, double beta2)
The inverse of the Vavilov cumulative probability density function.
Definition: VavilovFast.cxx:592
ROOT::Math::vavilov_fast_quantile_c
double vavilov_fast_quantile_c(double z, double kappa, double beta2)
The inverse of the complementary Vavilov cumulative probability density function.
Definition: VavilovFast.cxx:597
ROOT::VecOps::log
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1765
ROOT::VecOps::exp
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1761
y
Double_t y[n]
Definition: legend1.C:17
x
Double_t x[n]
Definition: legend1.C:17
n
const Int_t n
Definition: legend1.C:16
Math
Namespace for new Math classes and functions.
ROOT::Math::sqrt
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:281
ROOT
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
Definition: EExecutionPolicy.hxx:4
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
v
@ v
Definition: rootcling_impl.cxx:3667