ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups 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 <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 
58 VavilovFast::~VavilovFast()
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 
551 double VavilovFast::GetLambdaMin() const {
552  return fAC[0];
553 }
554 
555 double VavilovFast::GetLambdaMax() const {
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 
567 VavilovFast *VavilovFast::GetInstance() {
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
ROOT::Math::VavilovFast::GetBeta2
virtual double GetBeta2() const
Return the current value of .
Definition: VavilovFast.cxx:563
ROOT::Math::vavilov_fast_pdf
double vavilov_fast_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
Definition: VavilovFast.cxx:578
ROOT::Math::VavilovFast::Quantile_c
double Quantile_c(double z) const
Evaluate the inverse of the complementary Vavilov cummulative probability density function...
Definition: VavilovFast.cxx:541
ROOT::Math::VavilovFast::~VavilovFast
virtual ~VavilovFast()
Destructor.
Definition: VavilovFast.cxx:58
fu
static float fu(float x)
Definition: main.cpp:53
p3
static double p3(double t, double a, double b, double c, double d)
Definition: RooChebychev.cxx:98
ROOT::Math::VavilovFast::SetKappaBeta2
virtual void SetKappaBeta2(double kappa, double beta2)
Change and and recalculate coefficients if necessary.
Definition: VavilovFast.cxx:63
V1
Double_t * V1
Definition: memstatExample.C:59
h
TH1 * h
Definition: legend2.C:5
ROOT::Math::VavilovFast::GetLambdaMin
virtual double GetLambdaMin() const
Return the minimum value of for which is nonzero in the current approximation.
Definition: VavilovFast.cxx:551
ROOT::Math::VavilovFast::VavilovFast
VavilovFast(double kappa=1, double beta2=1)
Initialize an object to calculate the Vavilov distribution.
Definition: VavilovFast.cxx:52
ROOT::Math::VavilovFast::Cdf_c
double Cdf_c(double x) const
Evaluate the Vavilov complementary cummulative probability density function.
Definition: VavilovFast.cxx:440
ROOT::Math::VavilovFast::Cdf
double Cdf(double x) const
Evaluate the Vavilov cummulative probability density function.
Definition: VavilovFast.cxx:426
ROOT::Math::VavilovFast
Class describing a Vavilov distribution.
Definition: VavilovFast.h:116
sqrt
double sqrt(double)
x2
static const double x2[5]
Definition: RooGaussKronrodIntegrator1D.cxx:345
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Math::landau_pdf
double landau_pdf(double x, double xi=1, double x0=0)
Probability density function of the Landau distribution: with where .
Definition: PdfFuncMathCore.cxx:191
p2
static double p2(double t, double a, double b, double c)
Definition: RooChebychev.cxx:97
z
Float_t z[5]
Definition: Ifit.C:16
ROOT::Math::vavilov_fast_quantile_c
double vavilov_fast_quantile_c(double z, double kappa, double beta2)
The inverse of the complementary Vavilov cummulative probability density function.
Definition: VavilovFast.cxx:598
ROOT::Math::landau_cdf
double landau_cdf(double x, double xi=1, double x0=0)
Cumulative distribution function of the Landau distribution (lower tail).
Definition: ProbFuncMathCore.cxx:336
ROOT::Math::VavilovFast::GetInstance
static VavilovFast * GetInstance()
Returns a static instance of class VavilovFast.
Definition: VavilovFast.cxx:567
t
TThread * t[5]
Definition: threadsh1.C:13
V3
Double_t * V3
Definition: memstatExample.C:59
v
SVector< double, 2 > v
Definition: Dict.h:5
xy
XPoint xy[kMAXMK]
Definition: TGX11.cxx:122
ROOT::Math::vavilov_fast_cdf
double vavilov_fast_cdf(double x, double kappa, double beta2)
The Vavilov cummulative probability density function.
Definition: VavilovFast.cxx:583
y
Double_t y[n]
Definition: legend1.C:17
ROOT::Math::vavilov_fast_quantile
double vavilov_fast_quantile(double z, double kappa, double beta2)
The inverse of the Vavilov cummulative probability density function.
Definition: VavilovFast.cxx:593
ROOT::Math::VavilovFast::fgInstance
static VavilovFast * fgInstance
Definition: VavilovFast.h:279
s0
TSocket * s0
Definition: hserv2.C:36
ROOT::Math::VavilovFast::GetKappa
virtual double GetKappa() const
Return the current value of .
Definition: VavilovFast.cxx:559
ROOT::Math::VavilovFast::Quantile
double Quantile(double z) const
Evaluate the inverse of the Vavilov cummulative probability density function.
Definition: VavilovFast.cxx:480
exp
double exp(double)
V4
Double_t * V4
Definition: memstatExample.C:59
n
const Int_t n
Definition: legend1.C:16
log
double log(double)
ROOT::Math::vavilov_fast_cdf_c
double vavilov_fast_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cummulative probability density function.
Definition: VavilovFast.cxx:588
ROOT::Math::VavilovFast::GetLambdaMax
virtual double GetLambdaMax() const
Return the maximum value of for which is nonzero in the current approximation.
Definition: VavilovFast.cxx:555
V2
Double_t * V2
Definition: memstatExample.C:59
x3
static const double x3[11]
Definition: RooGaussKronrodIntegrator1D.cxx:373
ROOT::Math::VavilovFast::Pdf
double Pdf(double x) const
Evaluate the Vavilov probability density function.
Definition: VavilovFast.cxx:364