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"
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
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) {
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
551 return fAC[0];
552}
553
555 return fAC[8];
556}
557
558double VavilovFast::GetKappa() const {
559 return fKappa;
560}
561
562double VavilovFast::GetBeta2() const {
563 return fBeta2;
564}
565
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
#define s0(x)
Definition: RSha256.hxx:90
#define h(i)
Definition: RSha256.hxx:106
static const double x2[5]
static const double x3[11]
XPoint xy[kMAXMK]
Definition: TGX11.cxx:122
double sqrt(double)
double exp(double)
double log(double)
Class describing a Vavilov distribution.
Definition: VavilovFast.h:116
double Cdf(double x) const
Evaluate the Vavilov cumulative probability density function.
static VavilovFast * fgInstance
Definition: VavilovFast.h:279
virtual double GetKappa() const
Return the current value of .
virtual double GetLambdaMax() const
Return the maximum value of for which is nonzero in the current approximation.
virtual ~VavilovFast()
Destructor.
Definition: VavilovFast.cxx:57
double Cdf_c(double x) const
Evaluate the Vavilov complementary cumulative probability density function.
virtual double GetBeta2() const
Return the current value of .
double Quantile_c(double z) const
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
double Pdf(double x) const
Evaluate the Vavilov probability density function.
double Quantile(double z) const
Evaluate the inverse of the Vavilov cumulative probability density function.
VavilovFast(double kappa=1, double beta2=1)
Initialize an object to calculate the Vavilov distribution.
Definition: VavilovFast.cxx:51
virtual void SetKappaBeta2(double kappa, double beta2)
Change and and recalculate coefficients if necessary.
Definition: VavilovFast.cxx:62
virtual double GetLambdaMin() const
Return the minimum value of for which is nonzero in the current approximation.
static VavilovFast * GetInstance()
Returns a static instance of class VavilovFast.
double vavilov_fast_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
double landau_pdf(double x, double xi=1, double x0=0)
Probability density function of the Landau distribution:
double vavilov_fast_cdf(double x, double kappa, double beta2)
The Vavilov cumulative probability density function.
double vavilov_fast_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cumulative probability density function.
double landau_cdf(double x, double xi=1, double x0=0)
Cumulative distribution function of the Landau distribution (lower tail).
double vavilov_fast_quantile(double z, double kappa, double beta2)
The inverse of the Vavilov cumulative probability density function.
double vavilov_fast_quantile_c(double z, double kappa, double beta2)
The inverse of the complementary Vavilov cumulative probability density function.
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
VSD Structures.
Definition: StringConv.hxx:21
static constexpr double s