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