Logo ROOT   6.12/07
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 
45 namespace ROOT {
46 namespace Math {
47 
48 VavilovFast *VavilovFast::fgInstance = 0;
49 
50 
51 VavilovFast::VavilovFast(double kappa, double beta2)
52 {
53  SetKappaBeta2 (kappa, beta2);
54 }
55 
56 
58 {
59  // desctructor (clean up resources)
60 }
61 
62 void 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 
363 double 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]*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 
407 double 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 
425 double 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 
439 double VavilovFast::Cdf_c (double x) const {
440  return 1-Cdf(x);
441 }
442 
443 double 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 
461 double 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 
479 double 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]*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 
535 double VavilovFast::Quantile (double z, double kappa, double beta2) {
536  if (kappa != fKappa || beta2 != fBeta2) SetKappaBeta2 (kappa, beta2);
537  return Quantile (z);
538 }
539 
540 double 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 
545 double 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 
558 double VavilovFast::GetKappa() const {
559  return fKappa;
560 }
561 
562 double VavilovFast::GetBeta2() const {
563  return fBeta2;
564 }
565 
567  if (!fgInstance) fgInstance = new VavilovFast (1, 1);
568  return fgInstance;
569 }
570 
571 VavilovFast *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 
577 double vavilov_fast_pdf (double x, double kappa, double beta2) {
578  VavilovFast *vavilov = VavilovFast::GetInstance (kappa, beta2);
579  return vavilov->Pdf (x);
580 }
581 
582 double vavilov_fast_cdf (double x, double kappa, double beta2) {
583  VavilovFast *vavilov = VavilovFast::GetInstance (kappa, beta2);
584  return vavilov->Cdf (x);
585 }
586 
587 double 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 
592 double vavilov_fast_quantile (double z, double kappa, double beta2) {
593  VavilovFast *vavilov = VavilovFast::GetInstance (kappa, beta2);
594  return vavilov->Quantile (z);
595 }
596 
597 double 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
double vavilov_fast_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
virtual double GetLambdaMax() const
Return the maximum value of for which is nonzero in the current approximation.
virtual ~VavilovFast()
Destructor.
Definition: VavilovFast.cxx:57
static double p3(double t, double a, double b, double c, double d)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
virtual double GetKappa() const
Return the current value of .
virtual void SetKappaBeta2(double kappa, double beta2)
Change and and recalculate coefficients if necessary.
Definition: VavilovFast.cxx:62
double Quantile(double z) const
Evaluate the inverse of the Vavilov cumulative probability density function.
double Cdf(double x) const
Evaluate the Vavilov cumulative probability density function.
virtual double GetLambdaMin() const
Return the minimum value of for which is nonzero in the current approximation.
TH1 * h
Definition: legend2.C:5
VavilovFast(double kappa=1, double beta2=1)
Initialize an object to calculate the Vavilov distribution.
Definition: VavilovFast.cxx:51
Class describing a Vavilov distribution.
Definition: VavilovFast.h:116
double sqrt(double)
double Cdf_c(double x) const
Evaluate the Vavilov complementary cumulative probability density function.
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
double landau_pdf(double x, double xi=1, double x0=0)
Probability density function of the Landau distribution: with where .
static double p2(double t, double a, double b, double c)
double vavilov_fast_quantile_c(double z, double kappa, double beta2)
The inverse of the complementary Vavilov cumulative probability density function. ...
double landau_cdf(double x, double xi=1, double x0=0)
Cumulative distribution function of the Landau distribution (lower tail).
static VavilovFast * GetInstance()
Returns a static instance of class VavilovFast.
SVector< double, 2 > v
Definition: Dict.h:5
XPoint xy[kMAXMK]
Definition: TGX11.cxx:122
double vavilov_fast_cdf(double x, double kappa, double beta2)
The Vavilov cumulative probability density function.
Double_t y[n]
Definition: legend1.C:17
double vavilov_fast_quantile(double z, double kappa, double beta2)
The inverse of the Vavilov cumulative probability density function.
static constexpr double s
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
virtual double GetBeta2() const
Return the current value of .
static VavilovFast * fgInstance
Definition: VavilovFast.h:279
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 exp(double)
const Int_t n
Definition: legend1.C:16
double log(double)
double vavilov_fast_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cumulative probability density function.
static const double x3[11]