Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
57 VavilovFast::~VavilovFast()
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]*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 
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]*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 
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 
550 double VavilovFast::GetLambdaMin() const {
551  return fAC[0];
552 }
553 
554 double VavilovFast::GetLambdaMax() const {
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 
566 VavilovFast *VavilovFast::GetInstance() {
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