Logo ROOT   6.10/09
Reference Guide
VavilovTest.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 VavilovTest
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 "VavilovTest.h"
35 #include "Math/Vavilov.h"
36 #include "Math/VavilovAccurate.h"
37 #include "Math/VavilovFast.h"
38 //#include "Math/SpecFuncMathCore.h"
39 //#include "Math/SpecFuncMathMore.h"
40 
41 #include <cassert>
42 #include <iostream>
43 #include <cmath>
44 #include <iomanip>
45 #include <cmath>
46 #include <cstdlib>
47 #include <string>
48 #include <sstream>
49 
50 #ifndef roundl
51 #define roundl(x) (long double)((long long)((x == 0) ? 0.0 : ( (x) + ( ((x) > 0) ? 0.5 : -0.5) )))
52 #endif
53 
54 namespace ROOT {
55 namespace Math {
56 
57 static const double eu = 0.577215664901532860606; // Euler's constant
58 
59 
60 //kappa = 0.01
61 static double vavilovPdfValues0[45][12]={
62  {-3.00, 6.80E-4, 6.85E-4, 6.90E-4, 6.95E-4, 7.00E-4, 7.05E-4, 7.10E-4, 7.15E-4, 7.21E-4, 7.26E-4, 7.31E-4},
63  {-2.50, 9.73E-3, 9.80E-3, 9.87E-3, 9.93E-3, 1.00E-2, 1.01E-2, 1.01E-2, 1.02E-2, 1.03E-2, 1.03E-2, 1.04E-2},
64  {-2.00, 4.44E-2, 4.47E-2, 4.50E-2, 4.53E-2, 4.55E-2, 4.58E-2, 4.61E-2, 4.64E-2, 4.67E-2, 4.70E-2, 4.73E-2},
65  {-1.50, 1.02E-1, 1.02E-1, 1.03E-1, 1.03E-1, 1.04E-1, 1.04E-1, 1.05E-1, 1.06E-1, 1.06E-1, 1.07E-1, 1.07E-1},
66  {-1.00, 1.53E-1, 1.54E-1, 1.55E-1, 1.55E-1, 1.56E-1, 1.57E-1, 1.58E-1, 1.59E-1, 1.59E-1, 1.60E-1, 1.61E-1},
67  {-0.50, 1.79E-1, 1.80E-1, 1.81E-1, 1.82E-1, 1.83E-1, 1.83E-1, 1.84E-1, 1.85E-1, 1.86E-1, 1.87E-1, 1.88E-1},
68  { 0.00, 1.81E-1, 1.81E-1, 1.82E-1, 1.83E-1, 1.84E-1, 1.84E-1, 1.85E-1, 1.86E-1, 1.87E-1, 1.88E-1, 1.88E-1},
69  { 0.50, 1.67E-1, 1.68E-1, 1.68E-1, 1.69E-1, 1.69E-1, 1.70E-1, 1.71E-1, 1.71E-1, 1.72E-1, 1.73E-1, 1.73E-1},
70  { 1.00, 1.47E-1, 1.47E-1, 1.48E-1, 1.48E-1, 1.49E-1, 1.49E-1, 1.49E-1, 1.50E-1, 1.50E-1, 1.51E-1, 1.51E-1},
71  { 1.50, 1.25E-1, 1.26E-1, 1.26E-1, 1.26E-1, 1.27E-1, 1.27E-1, 1.27E-1, 1.28E-1, 1.28E-1, 1.29E-1, 1.29E-1},
72  { 2.00, 1.06E-1, 1.06E-1, 1.06E-1, 1.07E-1, 1.07E-1, 1.07E-1, 1.07E-1, 1.08E-1, 1.08E-1, 1.08E-1, 1.08E-1},
73  { 2.50, 8.91E-2, 8.93E-2, 8.94E-2, 8.96E-2, 8.97E-2, 8.99E-2, 9.00E-2, 9.02E-2, 9.03E-2, 9.04E-2, 9.06E-2},
74  { 3.00, 7.50E-2, 7.51E-2, 7.52E-2, 7.53E-2, 7.53E-2, 7.54E-2, 7.55E-2, 7.56E-2, 7.57E-2, 7.58E-2, 7.59E-2},
75  { 3.50, 6.34E-2, 6.34E-2, 6.34E-2, 6.35E-2, 6.35E-2, 6.36E-2, 6.36E-2, 6.36E-2, 6.37E-2, 6.37E-2, 6.37E-2},
76  { 4.00, 5.38E-2, 5.38E-2, 5.38E-2, 5.38E-2, 5.38E-2, 5.38E-2, 5.38E-2, 5.39E-2, 5.39E-2, 5.39E-2, 5.39E-2},
77  { 4.50, 4.60E-2, 4.60E-2, 4.60E-2, 4.59E-2, 4.59E-2, 4.59E-2, 4.59E-2, 4.59E-2, 4.58E-2, 4.58E-2, 4.58E-2},
78  { 5.00, 3.96E-2, 3.95E-2, 3.95E-2, 3.95E-2, 3.94E-2, 3.94E-2, 3.93E-2, 3.93E-2, 3.93E-2, 3.92E-2, 3.92E-2},
79  { 5.50, 3.43E-2, 3.42E-2, 3.42E-2, 3.41E-2, 3.41E-2, 3.40E-2, 3.40E-2, 3.39E-2, 3.39E-2, 3.38E-2, 3.38E-2},
80  { 6.00, 2.99E-2, 2.98E-2, 2.97E-2, 2.97E-2, 2.96E-2, 2.96E-2, 2.95E-2, 2.95E-2, 2.94E-2, 2.93E-2, 2.93E-2},
81  { 6.50, 2.62E-2, 2.61E-2, 2.61E-2, 2.60E-2, 2.59E-2, 2.59E-2, 2.58E-2, 2.57E-2, 2.57E-2, 2.56E-2, 2.55E-2},
82  { 7.00, 2.31E-2, 2.30E-2, 2.30E-2, 2.29E-2, 2.28E-2, 2.28E-2, 2.27E-2, 2.26E-2, 2.26E-2, 2.25E-2, 2.24E-2},
83  { 7.50, 2.05E-2, 2.04E-2, 2.04E-2, 2.03E-2, 2.02E-2, 2.01E-2, 2.01E-2, 2.00E-2, 1.99E-2, 1.99E-2, 1.98E-2},
84  { 8.00, 1.83E-2, 1.82E-2, 1.81E-2, 1.81E-2, 1.80E-2, 1.79E-2, 1.78E-2, 1.78E-2, 1.77E-2, 1.76E-2, 1.75E-2},
85  { 8.50, 1.64E-2, 1.63E-2, 1.62E-2, 1.62E-2, 1.61E-2, 1.60E-2, 1.59E-2, 1.59E-2, 1.58E-2, 1.57E-2, 1.56E-2},
86  { 9.00, 1.47E-2, 1.47E-2, 1.46E-2, 1.45E-2, 1.45E-2, 1.44E-2, 1.43E-2, 1.42E-2, 1.42E-2, 1.41E-2, 1.40E-2},
87  { 9.50, 1.33E-2, 1.33E-2, 1.32E-2, 1.31E-2, 1.30E-2, 1.30E-2, 1.29E-2, 1.28E-2, 1.27E-2, 1.27E-2, 1.26E-2},
88  {10.00, 1.21E-2, 1.20E-2, 1.20E-2, 1.19E-2, 1.18E-2, 1.17E-2, 1.17E-2, 1.16E-2, 1.15E-2, 1.14E-2, 1.14E-2},
89  {11.00, 1.01E-2, 1.00E-2, 9.94E-3, 9.87E-3, 9.80E-3, 9.73E-3, 9.66E-3, 9.59E-3, 9.51E-3, 9.44E-3, 9.37E-3},
90  {12.00, 8.51E-3, 8.44E-3, 8.37E-3, 8.31E-3, 8.24E-3, 8.17E-3, 8.10E-3, 8.03E-3, 7.96E-3, 7.89E-3, 7.82E-3},
91  {13.00, 7.26E-3, 7.20E-3, 7.14E-3, 7.07E-3, 7.00E-3, 6.94E-3, 6.87E-3, 6.81E-3, 6.74E-3, 6.67E-3, 6.61E-3},
92  {14.00, 6.27E-3, 6.20E-3, 6.14E-3, 6.08E-3, 6.02E-3, 5.95E-3, 5.89E-3, 5.83E-3, 5.76E-3, 5.70E-3, 5.64E-3},
93  {15.00, 5.46E-3, 5.40E-3, 5.34E-3, 5.28E-3, 5.22E-3, 5.16E-3, 5.10E-3, 5.04E-3, 4.97E-3, 4.91E-3, 4.85E-3},
94  {16.00, 4.79E-3, 4.73E-3, 4.67E-3, 4.62E-3, 4.56E-3, 4.50E-3, 4.44E-3, 4.39E-3, 4.33E-3, 4.27E-3, 4.21E-3},
95  {17.00, 4.23E-3, 4.18E-3, 4.12E-3, 4.07E-3, 4.01E-3, 3.96E-3, 3.90E-3, 3.85E-3, 3.79E-3, 3.73E-3, 3.68E-3},
96  {18.00, 3.77E-3, 3.72E-3, 3.66E-3, 3.61E-3, 3.56E-3, 3.50E-3, 3.45E-3, 3.40E-3, 3.34E-3, 3.29E-3, 3.24E-3},
97  {19.00, 3.37E-3, 3.32E-3, 3.27E-3, 3.22E-3, 3.17E-3, 3.12E-3, 3.07E-3, 3.02E-3, 2.97E-3, 2.91E-3, 2.86E-3},
98  {20.00, 3.04E-3, 2.99E-3, 2.94E-3, 2.89E-3, 2.84E-3, 2.79E-3, 2.74E-3, 2.69E-3, 2.64E-3, 2.60E-3, 2.55E-3},
99  {22.00, 2.49E-3, 2.45E-3, 2.40E-3, 2.36E-3, 2.31E-3, 2.27E-3, 2.22E-3, 2.18E-3, 2.13E-3, 2.09E-3, 2.04E-3},
100  {24.00, 2.08E-3, 2.04E-3, 2.00E-3, 1.96E-3, 1.92E-3, 1.88E-3, 1.83E-3, 1.79E-3, 1.75E-3, 1.71E-3, 1.66E-3},
101  {26.00, 1.77E-3, 1.73E-3, 1.69E-3, 1.65E-3, 1.61E-3, 1.57E-3, 1.53E-3, 1.49E-3, 1.45E-3, 1.41E-3, 1.37E-3},
102  {28.00, 1.51E-3, 1.48E-3, 1.44E-3, 1.41E-3, 1.37E-3, 1.33E-3, 1.30E-3, 1.26E-3, 1.22E-3, 1.18E-3, 1.15E-3},
103  {30.00, 1.31E-3, 1.28E-3, 1.24E-3, 1.21E-3, 1.18E-3, 1.14E-3, 1.11E-3, 1.07E-3, 1.04E-3, 1.00E-3, 9.68E-4},
104  {32.00, 1.15E-3, 1.12E-3, 1.08E-3, 1.05E-3, 1.02E-3, 9.87E-4, 9.54E-4, 9.22E-4, 8.89E-4, 8.56E-4, 8.24E-4},
105  {34.00, 1.01E-3, 9.81E-4, 9.51E-4, 9.21E-4, 8.90E-4, 8.60E-4, 8.29E-4, 7.98E-4, 7.68E-4, 7.37E-4, 7.06E-4},
106  {36.00, 8.98E-4, 8.70E-4, 8.41E-4, 8.12E-4, 7.83E-4, 7.54E-4, 7.25E-4, 6.96E-4, 6.67E-4, 6.38E-4, 6.09E-4}};
107 
108 
109 //kappa = 0.04
110 static double vavilovPdfValues1[42][12]={
111  {-3.00, 7.01E-4, 7.18E-4, 7.34E-4, 7.52E-4, 7.69E-4, 7.87E-4, 8.06E-4, 8.25E-4, 8.44E-4, 8.64E-4, 8.84E-4},
112  {-2.50, 1.00E-2, 1.02E-2, 1.05E-2, 1.07E-2, 1.09E-2, 1.12E-2, 1.14E-2, 1.16E-2, 1.19E-2, 1.21E-2, 1.24E-2},
113  {-2.00, 4.58E-2, 4.67E-2, 4.76E-2, 4.85E-2, 4.94E-2, 5.04E-2, 5.13E-2, 5.23E-2, 5.33E-2, 5.44E-2, 5.54E-2},
114  {-1.50, 1.05E-1, 1.06E-1, 1.08E-1, 1.10E-1, 1.12E-1, 1.14E-1, 1.16E-1, 1.18E-1, 1.20E-1, 1.22E-1, 1.24E-1},
115  {-1.00, 1.58E-1, 1.60E-1, 1.62E-1, 1.65E-1, 1.67E-1, 1.70E-1, 1.73E-1, 1.75E-1, 1.78E-1, 1.81E-1, 1.83E-1},
116  {-0.50, 1.85E-1, 1.87E-1, 1.89E-1, 1.92E-1, 1.95E-1, 1.97E-1, 2.00E-1, 2.02E-1, 2.05E-1, 2.08E-1, 2.10E-1},
117  { 0.00, 1.86E-1, 1.88E-1, 1.90E-1, 1.92E-1, 1.95E-1, 1.97E-1, 1.99E-1, 2.01E-1, 2.03E-1, 2.06E-1, 2.08E-1},
118  { 0.50, 1.72E-1, 1.74E-1, 1.75E-1, 1.77E-1, 1.78E-1, 1.80E-1, 1.82E-1, 1.83E-1, 1.85E-1, 1.87E-1, 1.88E-1},
119  { 1.00, 1.51E-1, 1.52E-1, 1.53E-1, 1.54E-1, 1.55E-1, 1.57E-1, 1.58E-1, 1.59E-1, 1.60E-1, 1.61E-1, 1.62E-1},
120  { 1.50, 1.29E-1, 1.30E-1, 1.31E-1, 1.31E-1, 1.32E-1, 1.33E-1, 1.33E-1, 1.34E-1, 1.34E-1, 1.35E-1, 1.36E-1},
121  { 2.00, 1.09E-1, 1.10E-1, 1.10E-1, 1.10E-1, 1.11E-1, 1.11E-1, 1.11E-1, 1.11E-1, 1.12E-1, 1.12E-1, 1.12E-1},
122  { 2.50, 9.18E-2, 9.19E-2, 9.20E-2, 9.21E-2, 9.22E-2, 9.22E-2, 9.23E-2, 9.23E-2, 9.23E-2, 9.24E-2, 9.24E-2},
123  { 3.00, 7.73E-2, 7.72E-2, 7.71E-2, 7.70E-2, 7.69E-2, 7.68E-2, 7.67E-2, 7.66E-2, 7.64E-2, 7.62E-2, 7.61E-2},
124  { 3.50, 6.53E-2, 6.51E-2, 6.49E-2, 6.47E-2, 6.45E-2, 6.42E-2, 6.40E-2, 6.37E-2, 6.34E-2, 6.32E-2, 6.29E-2},
125  { 4.00, 5.54E-2, 5.52E-2, 5.49E-2, 5.46E-2, 5.43E-2, 5.40E-2, 5.36E-2, 5.33E-2, 5.29E-2, 5.26E-2, 5.22E-2},
126  { 4.50, 4.74E-2, 4.71E-2, 4.67E-2, 4.64E-2, 4.60E-2, 4.56E-2, 4.53E-2, 4.49E-2, 4.45E-2, 4.40E-2, 4.36E-2},
127  { 5.00, 4.08E-2, 4.04E-2, 4.00E-2, 3.96E-2, 3.92E-2, 3.88E-2, 3.84E-2, 3.80E-2, 3.76E-2, 3.71E-2, 3.67E-2},
128  { 5.50, 3.53E-2, 3.49E-2, 3.45E-2, 3.41E-2, 3.37E-2, 3.33E-2, 3.28E-2, 3.24E-2, 3.19E-2, 3.15E-2, 3.10E-2},
129  { 6.00, 3.08E-2, 3.04E-2, 3.00E-2, 2.95E-2, 2.91E-2, 2.87E-2, 2.82E-2, 2.78E-2, 2.73E-2, 2.69E-2, 2.64E-2},
130  { 6.50, 2.70E-2, 2.66E-2, 2.62E-2, 2.57E-2, 2.53E-2, 2.49E-2, 2.44E-2, 2.40E-2, 2.35E-2, 2.31E-2, 2.26E-2},
131  { 7.00, 2.38E-2, 2.34E-2, 2.30E-2, 2.26E-2, 2.21E-2, 2.17E-2, 2.13E-2, 2.08E-2, 2.04E-2, 1.99E-2, 1.94E-2},
132  { 7.50, 2.11E-2, 2.07E-2, 2.03E-2, 1.99E-2, 1.95E-2, 1.90E-2, 1.86E-2, 1.82E-2, 1.77E-2, 1.73E-2, 1.68E-2},
133  { 8.00, 1.88E-2, 1.84E-2, 1.80E-2, 1.76E-2, 1.72E-2, 1.68E-2, 1.64E-2, 1.59E-2, 1.55E-2, 1.51E-2, 1.46E-2},
134  { 8.50, 1.69E-2, 1.65E-2, 1.61E-2, 1.57E-2, 1.53E-2, 1.49E-2, 1.45E-2, 1.40E-2, 1.36E-2, 1.32E-2, 1.27E-2},
135  { 9.00, 1.52E-2, 1.48E-2, 1.44E-2, 1.40E-2, 1.36E-2, 1.32E-2, 1.28E-2, 1.24E-2, 1.20E-2, 1.16E-2, 1.12E-2},
136  { 9.50, 1.37E-2, 1.34E-2, 1.30E-2, 1.26E-2, 1.22E-2, 1.18E-2, 1.14E-2, 1.10E-2, 1.06E-2, 1.02E-2, 9.81E-3},
137  {10.00, 1.25E-2, 1.21E-2, 1.17E-2, 1.14E-2, 1.10E-2, 1.06E-2, 1.02E-2, 9.84E-3, 9.45E-3, 9.05E-3, 8.65E-3},
138  {11.00, 1.04E-2, 1.00E-2, 9.70E-3, 9.35E-3, 8.99E-3, 8.63E-3, 8.27E-3, 7.91E-3, 7.54E-3, 7.17E-3, 6.79E-3},
139  {12.00, 8.77E-3, 8.44E-3, 8.11E-3, 7.78E-3, 7.45E-3, 7.11E-3, 6.77E-3, 6.43E-3, 6.08E-3, 5.73E-3, 5.38E-3},
140  {13.00, 7.49E-3, 7.18E-3, 6.87E-3, 6.56E-3, 6.24E-3, 5.92E-3, 5.60E-3, 5.28E-3, 4.95E-3, 4.63E-3, 4.30E-3},
141  {14.00, 6.46E-3, 6.17E-3, 5.87E-3, 5.58E-3, 5.28E-3, 4.98E-3, 4.68E-3, 4.38E-3, 4.07E-3, 3.76E-3, 3.45E-3},
142  {15.00, 5.62E-3, 5.35E-3, 5.07E-3, 4.79E-3, 4.51E-3, 4.22E-3, 3.94E-3, 3.65E-3, 3.36E-3, 3.07E-3, 2.78E-3},
143  {16.00, 4.93E-3, 4.67E-3, 4.41E-3, 4.14E-3, 3.88E-3, 3.61E-3, 3.34E-3, 3.07E-3, 2.80E-3, 2.52E-3, 2.24E-3},
144  {17.00, 4.36E-3, 4.11E-3, 3.86E-3, 3.61E-3, 3.36E-3, 3.10E-3, 2.85E-3, 2.59E-3, 2.33E-3, 2.07E-3, 1.81E-3},
145  {18.00, 3.88E-3, 3.65E-3, 3.41E-3, 3.17E-3, 2.93E-3, 2.69E-3, 2.44E-3, 2.20E-3, 1.95E-3, 1.71E-3, 1.46E-3},
146  {19.00, 3.48E-3, 3.25E-3, 3.02E-3, 2.79E-3, 2.57E-3, 2.34E-3, 2.10E-3, 1.87E-3, 1.64E-3, 1.41E-3, 1.17E-3},
147  {20.00, 3.13E-3, 2.91E-3, 2.70E-3, 2.48E-3, 2.26E-3, 2.04E-3, 1.82E-3, 1.60E-3, 1.38E-3, 1.16E-3, 9.32E-4},
148  {22.00, 2.57E-3, 2.37E-3, 2.17E-3, 1.98E-3, 1.78E-3, 1.58E-3, 1.37E-3, 1.17E-3, 9.71E-4, 7.69E-4, 5.66E-4},
149  {24.00, 1.97E-3, 1.80E-3, 1.63E-3, 1.47E-3, 1.30E-3, 1.13E-3, 9.69E-4, 8.04E-4, 6.39E-4, 4.75E-4, 3.11E-4},
150  {26.00, 1.14E-3, 1.04E-3, 9.34E-4, 8.32E-4, 7.31E-4, 6.31E-4, 5.34E-4, 4.38E-4, 3.44E-4, 2.52E-4, 1.61E-4},
151  {28.00, 6.50E-4, 5.85E-4, 5.22E-4, 4.61E-4, 4.02E-4, 3.45E-4, 2.89E-4, 2.35E-4, 1.84E-4, 1.34E-4, 8.60E-5},
152  {30.00, 3.95E-4, 3.53E-4, 3.12E-4, 2.73E-4, 2.36E-4, 2.00E-4, 1.66E-4, 1.34E-4, 1.03E-4, 7.46E-5, 4.76E-5}};
153 
154 
155 //kappa = 0.07
156 static double vavilovPdfValues2[41][12]={
157  {-3.50, 0, 0, 0, 0, 0, 0, 1.01E-5, 1.06E-5, 1.10E-5, 1.14E-5, 1.19E-5},
158  {-3.00, 7.23E-4, 7.50E-4, 7.78E-4, 8.07E-4, 8.37E-4, 8.68E-4, 9.00E-4, 9.34E-4, 9.69E-4, 1.01E-3, 1.04E-3},
159  {-2.50, 1.03E-2, 1.07E-2, 1.10E-2, 1.14E-2, 1.18E-2, 1.22E-2, 1.26E-2, 1.30E-2, 1.35E-2, 1.39E-2, 1.44E-2},
160  {-2.00, 4.72E-2, 4.86E-2, 5.00E-2, 5.15E-2, 5.31E-2, 5.47E-2, 5.63E-2, 5.80E-2, 5.98E-2, 6.15E-2, 6.34E-2},
161  {-1.50, 1.08E-1, 1.11E-1, 1.14E-1, 1.17E-1, 1.20E-1, 1.23E-1, 1.26E-1, 1.29E-1, 1.33E-1, 1.36E-1, 1.40E-1},
162  {-1.00, 1.62E-1, 1.66E-1, 1.70E-1, 1.74E-1, 1.78E-1, 1.82E-1, 1.86E-1, 1.90E-1, 1.94E-1, 1.99E-1, 2.03E-1},
163  {-0.50, 1.90E-1, 1.94E-1, 1.98E-1, 2.01E-1, 2.05E-1, 2.09E-1, 2.13E-1, 2.17E-1, 2.21E-1, 2.25E-1, 2.30E-1},
164  { 0.00, 1.92E-1, 1.95E-1, 1.98E-1, 2.01E-1, 2.04E-1, 2.07E-1, 2.10E-1, 2.14E-1, 2.17E-1, 2.20E-1, 2.23E-1},
165  { 0.50, 1.77E-1, 1.79E-1, 1.82E-1, 1.84E-1, 1.86E-1, 1.88E-1, 1.90E-1, 1.92E-1, 1.95E-1, 1.97E-1, 1.99E-1},
166  { 1.00, 1.56E-1, 1.57E-1, 1.58E-1, 1.60E-1, 1.61E-1, 1.62E-1, 1.64E-1, 1.65E-1, 1.66E-1, 1.67E-1, 1.68E-1},
167  { 1.50, 1.33E-1, 1.34E-1, 1.35E-1, 1.35E-1, 1.36E-1, 1.36E-1, 1.37E-1, 1.37E-1, 1.38E-1, 1.38E-1, 1.39E-1},
168  { 2.00, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1, 1.13E-1},
169  { 2.50, 9.46E-2, 9.44E-2, 9.42E-2, 9.40E-2, 9.37E-2, 9.33E-2, 9.30E-2, 9.26E-2, 9.21E-2, 9.16E-2, 9.11E-2},
170  { 3.00, 7.96E-2, 7.92E-2, 7.87E-2, 7.82E-2, 7.77E-2, 7.71E-2, 7.65E-2, 7.58E-2, 7.51E-2, 7.44E-2, 7.36E-2},
171  { 3.50, 6.73E-2, 6.67E-2, 6.60E-2, 6.53E-2, 6.46E-2, 6.39E-2, 6.31E-2, 6.23E-2, 6.14E-2, 6.06E-2, 5.96E-2},
172  { 4.00, 5.71E-2, 5.64E-2, 5.57E-2, 5.49E-2, 5.41E-2, 5.32E-2, 5.23E-2, 5.14E-2, 5.05E-2, 4.95E-2, 4.85E-2},
173  { 4.50, 4.88E-2, 4.80E-2, 4.72E-2, 4.64E-2, 4.55E-2, 4.46E-2, 4.37E-2, 4.27E-2, 4.17E-2, 4.07E-2, 3.97E-2},
174  { 5.00, 4.20E-2, 4.12E-2, 4.03E-2, 3.95E-2, 3.86E-2, 3.76E-2, 3.67E-2, 3.57E-2, 3.47E-2, 3.36E-2, 3.26E-2},
175  { 5.50, 3.64E-2, 3.55E-2, 3.47E-2, 3.38E-2, 3.29E-2, 3.19E-2, 3.10E-2, 3.00E-2, 2.90E-2, 2.80E-2, 2.69E-2},
176  { 6.00, 3.17E-2, 3.09E-2, 3.00E-2, 2.91E-2, 2.82E-2, 2.73E-2, 2.63E-2, 2.53E-2, 2.44E-2, 2.33E-2, 2.23E-2},
177  { 6.50, 2.78E-2, 2.70E-2, 2.61E-2, 2.52E-2, 2.43E-2, 2.34E-2, 2.25E-2, 2.15E-2, 2.06E-2, 1.96E-2, 1.86E-2},
178  { 7.00, 2.45E-2, 2.37E-2, 2.29E-2, 2.20E-2, 2.11E-2, 2.02E-2, 1.93E-2, 1.84E-2, 1.74E-2, 1.65E-2, 1.55E-2},
179  { 7.50, 2.18E-2, 2.09E-2, 2.01E-2, 1.93E-2, 1.84E-2, 1.76E-2, 1.67E-2, 1.58E-2, 1.49E-2, 1.39E-2, 1.30E-2},
180  { 8.00, 1.94E-2, 1.86E-2, 1.78E-2, 1.70E-2, 1.62E-2, 1.53E-2, 1.45E-2, 1.36E-2, 1.27E-2, 1.18E-2, 1.09E-2},
181  { 8.50, 1.74E-2, 1.66E-2, 1.58E-2, 1.50E-2, 1.42E-2, 1.34E-2, 1.26E-2, 1.18E-2, 1.09E-2, 1.00E-2, 9.17E-3},
182  { 9.00, 1.57E-2, 1.49E-2, 1.41E-2, 1.34E-2, 1.26E-2, 1.18E-2, 1.10E-2, 1.02E-2, 9.39E-3, 8.56E-3, 7.72E-3},
183  { 9.50, 1.42E-2, 1.34E-2, 1.27E-2, 1.19E-2, 1.12E-2, 1.04E-2, 9.66E-3, 8.88E-3, 8.09E-3, 7.30E-3, 6.49E-3},
184  {10.00, 1.28E-2, 1.21E-2, 1.14E-2, 1.07E-2, 9.99E-3, 9.25E-3, 8.51E-3, 7.75E-3, 6.99E-3, 6.23E-3, 5.45E-3},
185  {11.00, 1.07E-2, 1.00E-2, 9.37E-3, 8.70E-3, 8.02E-3, 7.34E-3, 6.64E-3, 5.95E-3, 5.24E-3, 4.53E-3, 3.81E-3},
186  {12.00, 9.01E-3, 8.39E-3, 7.77E-3, 7.14E-3, 6.50E-3, 5.87E-3, 5.22E-3, 4.58E-3, 3.93E-3, 3.27E-3, 2.61E-3},
187  {13.00, 7.35E-3, 6.79E-3, 6.23E-3, 5.68E-3, 5.12E-3, 4.55E-3, 3.99E-3, 3.43E-3, 2.86E-3, 2.29E-3, 1.73E-3},
188  {14.00, 5.54E-3, 5.08E-3, 4.63E-3, 4.18E-3, 3.73E-3, 3.28E-3, 2.84E-3, 2.40E-3, 1.97E-3, 1.54E-3, 1.11E-3},
189  {15.00, 3.97E-3, 3.61E-3, 3.27E-3, 2.92E-3, 2.59E-3, 2.26E-3, 1.93E-3, 1.62E-3, 1.31E-3, 1.00E-3, 7.05E-4},
190  {16.00, 2.81E-3, 2.54E-3, 2.28E-3, 2.02E-3, 1.78E-3, 1.54E-3, 1.30E-3, 1.08E-3, 8.59E-4, 6.49E-4, 4.48E-4},
191  {17.00, 2.00E-3, 1.80E-3, 1.60E-3, 1.41E-3, 1.23E-3, 1.05E-3, 8.82E-4, 7.21E-4, 5.68E-4, 4.23E-4, 2.86E-4},
192  {18.00, 1.44E-3, 1.29E-3, 1.14E-3, 9.93E-4, 8.56E-4, 7.26E-4, 6.03E-4, 4.87E-4, 3.79E-4, 2.77E-4, 1.83E-4},
193  {19.00, 1.05E-3, 9.32E-4, 8.17E-4, 7.08E-4, 6.05E-4, 5.08E-4, 4.17E-4, 3.33E-4, 2.55E-4, 1.83E-4, 1.18E-4},
194  {20.00, 7.77E-4, 6.83E-4, 5.94E-4, 5.10E-4, 4.32E-4, 3.59E-4, 2.91E-4, 2.29E-4, 1.72E-4, 1.21E-4, 7.58E-5},
195  {22.00, 4.31E-4, 3.73E-4, 3.20E-4, 2.70E-4, 2.24E-4, 1.82E-4, 1.44E-4, 1.10E-4, 7.96E-5, 5.34E-5, 3.11E-5},
196  {24.00, 2.40E-4, 2.05E-4, 1.73E-4, 1.43E-4, 1.16E-4, 9.22E-5, 7.08E-5, 5.32E-5, 3.63E-5, 2.30E-5, 1.23E-5},
197  {26.00, 1.30E-4, 1.09E-4, 9.02E-5, 7.34E-5, 5.84E-5, 4.51E-5, 3.37E-5, 2.39E-5, 1.58E-5, 0, 0}};
198 
199 
200 //kappa = 0.10
201 static double vavilovPdfValues3[41][12]={
202  {-3.50, 0, 0, 0, 0, 1.00E-5, 1.06E-5, 1.11E-5, 1.18E-5, 1.24E-5, 1.31E-5, 1.38E-5},
203  {-3.00, 7.45E-4, 7.82E-4, 8.21E-4, 8.62E-4, 9.05E-4, 9.50E-4, 9.98E-4, 1.05E-3, 1.10E-3, 1.15E-3, 1.21E-3},
204  {-2.50, 1.07E-2, 1.11E-2, 1.16E-2, 1.21E-2, 1.27E-2, 1.33E-2, 1.38E-2, 1.45E-2, 1.51E-2, 1.58E-2, 1.65E-2},
205  {-2.00, 4.86E-2, 5.05E-2, 5.25E-2, 5.46E-2, 5.67E-2, 5.90E-2, 6.13E-2, 6.37E-2, 6.62E-2, 6.88E-2, 7.15E-2},
206  {-1.50, 1.11E-1, 1.15E-1, 1.19E-1, 1.23E-1, 1.27E-1, 1.32E-1, 1.36E-1, 1.41E-1, 1.45E-1, 1.50E-1, 1.55E-1},
207  {-1.00, 1.67E-1, 1.72E-1, 1.77E-1, 1.82E-1, 1.88E-1, 1.93E-1, 1.99E-1, 2.04E-1, 2.10E-1, 2.16E-1, 2.22E-1},
208  {-0.50, 1.96E-1, 2.01E-1, 2.06E-1, 2.10E-1, 2.15E-1, 2.20E-1, 2.26E-1, 2.31E-1, 2.36E-1, 2.42E-1, 2.47E-1},
209  { 0.00, 1.98E-1, 2.01E-1, 2.05E-1, 2.09E-1, 2.13E-1, 2.17E-1, 2.21E-1, 2.25E-1, 2.28E-1, 2.32E-1, 2.37E-1},
210  { 0.50, 1.83E-1, 1.85E-1, 1.88E-1, 1.90E-1, 1.93E-1, 1.95E-1, 1.98E-1, 2.00E-1, 2.02E-1, 2.05E-1, 2.07E-1},
211  { 1.00, 1.60E-1, 1.62E-1, 1.63E-1, 1.65E-1, 1.66E-1, 1.67E-1, 1.68E-1, 1.69E-1, 1.70E-1, 1.71E-1, 1.72E-1},
212  { 1.50, 1.37E-1, 1.38E-1, 1.38E-1, 1.39E-1, 1.39E-1, 1.39E-1, 1.39E-1, 1.39E-1, 1.39E-1, 1.39E-1, 1.39E-1},
213  { 2.00, 1.16E-1, 1.16E-1, 1.16E-1, 1.15E-1, 1.15E-1, 1.14E-1, 1.14E-1, 1.13E-1, 1.13E-1, 1.12E-1, 1.11E-1},
214  { 2.50, 9.75E-2, 9.69E-2, 9.62E-2, 9.54E-2, 9.45E-2, 9.36E-2, 9.26E-2, 9.16E-2, 9.04E-2, 8.92E-2, 8.79E-2},
215  { 3.00, 8.21E-2, 8.11E-2, 8.01E-2, 7.90E-2, 7.79E-2, 7.66E-2, 7.54E-2, 7.40E-2, 7.26E-2, 7.10E-2, 6.94E-2},
216  { 3.50, 6.93E-2, 6.82E-2, 6.70E-2, 6.57E-2, 6.43E-2, 6.29E-2, 6.15E-2, 5.99E-2, 5.83E-2, 5.67E-2, 5.49E-2},
217  { 4.00, 5.89E-2, 5.76E-2, 5.63E-2, 5.49E-2, 5.34E-2, 5.19E-2, 5.04E-2, 4.88E-2, 4.71E-2, 4.53E-2, 4.35E-2},
218  { 4.50, 5.03E-2, 4.90E-2, 4.76E-2, 4.61E-2, 4.46E-2, 4.31E-2, 4.15E-2, 3.99E-2, 3.82E-2, 3.64E-2, 3.46E-2},
219  { 5.00, 4.33E-2, 4.19E-2, 4.05E-2, 3.90E-2, 3.75E-2, 3.60E-2, 3.44E-2, 3.27E-2, 3.11E-2, 2.93E-2, 2.75E-2},
220  { 5.50, 3.75E-2, 3.61E-2, 3.47E-2, 3.32E-2, 3.17E-2, 3.02E-2, 2.86E-2, 2.70E-2, 2.54E-2, 2.37E-2, 2.20E-2},
221  { 6.00, 3.27E-2, 3.13E-2, 2.99E-2, 2.85E-2, 2.70E-2, 2.55E-2, 2.40E-2, 2.24E-2, 2.08E-2, 1.92E-2, 1.75E-2},
222  { 6.50, 2.86E-2, 2.73E-2, 2.59E-2, 2.45E-2, 2.31E-2, 2.17E-2, 2.02E-2, 1.87E-2, 1.71E-2, 1.55E-2, 1.39E-2},
223  { 7.00, 2.53E-2, 2.40E-2, 2.26E-2, 2.13E-2, 1.99E-2, 1.85E-2, 1.70E-2, 1.56E-2, 1.41E-2, 1.26E-2, 1.11E-2},
224  { 7.50, 2.24E-2, 2.11E-2, 1.98E-2, 1.85E-2, 1.72E-2, 1.58E-2, 1.45E-2, 1.31E-2, 1.16E-2, 1.02E-2, 8.73E-3},
225  { 8.00, 1.98E-2, 1.86E-2, 1.74E-2, 1.61E-2, 1.48E-2, 1.35E-2, 1.22E-2, 1.09E-2, 9.57E-3, 8.21E-3, 6.83E-3},
226  { 8.50, 1.74E-2, 1.62E-2, 1.51E-2, 1.39E-2, 1.27E-2, 1.15E-2, 1.03E-2, 9.04E-3, 7.80E-3, 6.55E-3, 5.29E-3},
227  { 9.00, 1.50E-2, 1.39E-2, 1.28E-2, 1.18E-2, 1.07E-2, 9.57E-3, 8.47E-3, 7.37E-3, 6.27E-3, 5.16E-3, 4.06E-3},
228  { 9.50, 1.27E-2, 1.17E-2, 1.07E-2, 9.77E-3, 8.80E-3, 7.84E-3, 6.88E-3, 5.92E-3, 4.97E-3, 4.03E-3, 3.09E-3},
229  {10.00, 1.05E-2, 9.69E-3, 8.84E-3, 7.99E-3, 7.16E-3, 6.33E-3, 5.51E-3, 4.70E-3, 3.90E-3, 3.11E-3, 2.33E-3},
230  {11.00, 7.12E-3, 6.48E-3, 5.85E-3, 5.23E-3, 4.62E-3, 4.03E-3, 3.45E-3, 2.89E-3, 2.35E-3, 1.83E-3, 1.32E-3},
231  {12.00, 4.77E-3, 4.30E-3, 3.84E-3, 3.39E-3, 2.96E-3, 2.54E-3, 2.15E-3, 1.77E-3, 1.41E-3, 1.07E-3, 7.45E-4},
232  {13.00, 3.21E-3, 2.86E-3, 2.53E-3, 2.21E-3, 1.90E-3, 1.61E-3, 1.34E-3, 1.08E-3, 8.42E-4, 6.21E-4, 4.18E-4},
233  {14.00, 2.18E-3, 1.92E-3, 1.68E-3, 1.45E-3, 1.23E-3, 1.03E-3, 8.38E-4, 6.65E-4, 5.06E-4, 3.62E-4, 2.34E-4},
234  {15.00, 1.49E-3, 1.30E-3, 1.12E-3, 9.53E-4, 7.99E-4, 6.57E-4, 5.27E-4, 4.09E-4, 3.04E-4, 2.11E-4, 1.30E-4},
235  {16.00, 1.01E-3, 8.76E-4, 7.47E-4, 6.28E-4, 5.19E-4, 4.20E-4, 3.31E-4, 2.51E-4, 1.81E-4, 1.21E-4, 7.12E-5},
236  {17.00, 6.88E-4, 5.88E-4, 4.96E-4, 4.12E-4, 3.35E-4, 2.67E-4, 2.06E-4, 1.53E-4, 1.07E-4, 6.91E-5, 3.84E-5},
237  {18.00, 4.60E-4, 3.89E-4, 3.24E-4, 2.66E-4, 2.13E-4, 1.67E-4, 1.26E-4, 9.16E-5, 6.25E-5, 3.87E-5, 2.03E-5},
238  {19.00, 3.01E-4, 2.52E-4, 2.07E-4, 1.68E-4, 1.33E-4, 1.02E-4, 7.62E-5, 5.41E-5, 3.59E-5, 2.15E-5, 1.07E-5},
239  {20.00, 1.94E-4, 1.61E-4, 1.31E-4, 1.05E-4, 8.18E-5, 6.21E-5, 4.55E-5, 3.17E-5, 2.06E-5, 1.20E-5, 0},
240  {22.00, 7.97E-5, 6.47E-5, 5.16E-5, 4.03E-5, 3.07E-5, 2.27E-5, 1.61E-5, 1.08E-5, 0, 0, 0},
241  {24.00, 3.26E-5, 2.59E-5, 2.01E-5, 1.53E-5, 1.12E-5, 0, 0, 0, 0, 0, 0},
242  {26.00, 1.33E-5, 1.04E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
243 
244 
245 // kappa = 0.40
246 static double vavilovPdfValues4[28][12]={
247  {-3.50, 1.07E-5, 1.26E-5, 1.47E-5, 1.73E-5, 2.03E-5, 2.37E-5, 2.78E-5, 3.26E-5, 3.82E-5, 4.48E-5, 5.24E-5},
248  {-3.00, 1.00E-3, 1.16E-3, 1.33E-3, 1.53E-3, 1.75E-3, 2.02E-3, 2.32E-3, 2.66E-3, 3.05E-3, 3.50E-3, 4.02E-3},
249  {-2.50, 1.44E-2, 1.62E-2, 1.83E-2, 2.06E-2, 2.32E-2, 2.61E-2, 2.93E-2, 3.30E-2, 3.71E-2, 4.17E-2, 4.68E-2},
250  {-2.00, 6.56E-2, 7.25E-2, 8.00E-2, 8.83E-2, 9.74E-2, 1.07E-1, 1.18E-1, 1.30E-1, 1.43E-1, 1.57E-1, 1.73E-1},
251  {-1.50, 1.50E-1, 1.62E-1, 1.76E-1, 1.90E-1, 2.05E-1, 2.21E-1, 2.38E-1, 2.57E-1, 2.76E-1, 2.97E-1, 3.19E-1},
252  {-1.00, 2.26E-1, 2.40E-1, 2.54E-1, 2.69E-1, 2.84E-1, 3.00E-1, 3.16E-1, 3.32E-1, 3.49E-1, 3.66E-1, 3.83E-1},
253  {-0.50, 2.65E-1, 2.75E-1, 2.85E-1, 2.96E-1, 3.05E-1, 3.15E-1, 3.24E-1, 3.32E-1, 3.40E-1, 3.47E-1, 3.52E-1},
254  { 0.00, 2.66E-1, 2.71E-1, 2.76E-1, 2.79E-1, 2.82E-1, 2.84E-1, 2.84E-1, 2.84E-1, 2.82E-1, 2.78E-1, 2.73E-1},
255  { 0.50, 2.44E-1, 2.43E-1, 2.42E-1, 2.39E-1, 2.36E-1, 2.31E-1, 2.25E-1, 2.18E-1, 2.09E-1, 1.99E-1, 1.88E-1},
256  { 1.00, 2.07E-1, 2.03E-1, 1.97E-1, 1.91E-1, 1.83E-1, 1.75E-1, 1.65E-1, 1.55E-1, 1.44E-1, 1.31E-1, 1.18E-1},
257  { 1.50, 1.66E-1, 1.59E-1, 1.51E-1, 1.43E-1, 1.34E-1, 1.24E-1, 1.14E-1, 1.03E-1, 9.21E-2, 8.05E-2, 6.86E-2},
258  { 2.00, 1.26E-1, 1.18E-1, 1.10E-1, 1.01E-1, 9.25E-2, 8.35E-2, 7.43E-2, 6.51E-2, 5.58E-2, 4.66E-2, 3.75E-2},
259  { 2.50, 9.11E-2, 8.37E-2, 7.62E-2, 6.87E-2, 6.11E-2, 5.36E-2, 4.63E-2, 3.91E-2, 3.22E-2, 2.56E-2, 1.94E-2},
260  { 3.00, 6.35E-2, 5.71E-2, 5.09E-2, 4.48E-2, 3.88E-2, 3.31E-2, 2.77E-2, 2.26E-2, 1.78E-2, 1.35E-2, 9.60E-3},
261  { 3.50, 4.28E-2, 3.77E-2, 3.29E-2, 2.82E-2, 2.39E-2, 1.98E-2, 1.60E-2, 1.26E-2, 9.52E-3, 6.84E-3, 4.55E-3},
262  { 4.00, 2.79E-2, 2.41E-2, 2.06E-2, 1.72E-2, 1.42E-2, 1.14E-2, 8.95E-3, 6.78E-3, 4.91E-3, 3.35E-3, 2.08E-3},
263  { 4.50, 1.77E-2, 1.50E-2, 1.25E-2, 1.02E-2, 8.21E-3, 6.42E-3, 4.87E-3, 3.55E-3, 2.46E-3, 1.59E-3, 9.22E-4},
264  { 5.00, 1.10E-2, 9.10E-3, 7.42E-3, 5.93E-3, 4.63E-3, 3.51E-3, 2.58E-3, 1.81E-3, 1.20E-3, 7.34E-4, 3.96E-4},
265  { 5.50, 6.65E-3, 5.39E-3, 4.30E-3, 3.35E-3, 2.55E-3, 1.88E-3, 1.33E-3, 9.02E-4, 5.72E-4, 3.30E-4, 1.66E-4},
266  { 6.00, 3.93E-3, 3.13E-3, 2.44E-3, 1.85E-3, 1.37E-3, 9.83E-4, 6.75E-4, 4.39E-4, 2.66E-4, 1.45E-4, 6.73E-5},
267  { 6.50, 2.28E-3, 1.78E-3, 1.35E-3, 1.01E-3, 7.25E-4, 5.04E-4, 3.34E-4, 2.09E-4, 1.21E-4, 6.24E-5, 2.68E-5},
268  { 7.00, 1.30E-3, 9.91E-4, 7.38E-4, 5.35E-4, 3.76E-4, 2.53E-4, 1.63E-4, 9.79E-5, 5.41E-5, 2.64E-5, 1.05E-5},
269  { 7.50, 7.27E-4, 5.43E-4, 3.96E-4, 2.80E-4, 1.91E-4, 1.25E-4, 7.76E-5, 4.49E-5, 2.36E-5, 1.08E-5, 0},
270  { 8.00, 4.00E-4, 2.92E-4, 2.08E-4, 1.44E-4, 9.57E-5, 6.08E-5, 3.64E-5, 2.02E-5, 1.02E-5, 0, 0},
271  { 8.50, 2.17E-4, 1.55E-4, 1.08E-4, 7.29E-5, 4.72E-5, 2.91E-5, 1.69E-5, 0, 0, 0, 0},
272  { 9.00, 1.16E-4, 8.12E-5, 5.53E-5, 3.63E-5, 2.29E-5, 1.37E-5, 0, 0, 0, 0, 0},
273  { 9.50, 6.08E-5, 4.19E-5, 2.79E-5, 1.79E-5, 1.09E-5, 0, 0, 0, 0, 0, 0},
274  {10.00, 3.17E-5, 2.13E-5, 1.39E-5, 0, 0, 0, 0, 0, 0, 0, 0}};
275 
276 
277 //kappa = 0.70
278 static double vavilovPdfValues5[22][12]={
279  {-3.50, 1.44E-5, 1.83E-5, 2.32E-5, 2.95E-5, 3.75E-5, 4.76E-5, 6.04E-5, 7.69E-5, 9.72E-5, 1.23E-4, 1.56E-4},
280  {-3.00, 1.36E-3, 1.67E-3, 2.04E-3, 2.50E-3, 3.07E-3, 3.76E-3, 4.60E-3, 5.62E-3, 6.87E-3, 8.39E-3, 1.02E-2},
281  {-2.50, 1.94E-2, 2.30E-2, 2.72E-2, 3.22E-2, 3.81E-2, 4.49E-2, 5.29E-2, 6.23E-2, 7.33E-2, 8.61E-2, 1.01E-1},
282  {-2.00, 8.86E-2, 1.01E-1, 1.16E-1, 1.32E-1, 1.50E-1, 1.71E-1, 1.94E-1, 2.19E-1, 2.47E-1, 2.79E-1, 3.13E-1},
283  {-1.50, 2.02E-1, 2.23E-1, 2.46E-1, 2.70E-1, 2.96E-1, 3.23E-1, 3.52E-1, 3.82E-1, 4.13E-1, 4.44E-1, 4.76E-1},
284  {-1.00, 3.03E-1, 3.23E-1, 3.42E-1, 3.62E-1, 3.81E-1, 3.99E-1, 4.15E-1, 4.30E-1, 4.42E-1, 4.52E-1, 4.57E-1},
285  {-0.50, 3.44E-1, 3.54E-1, 3.62E-1, 3.68E-1, 3.71E-1, 3.72E-1, 3.70E-1, 3.64E-1, 3.54E-1, 3.40E-1, 3.22E-1},
286  { 0.00, 3.23E-1, 3.20E-1, 3.15E-1, 3.09E-1, 2.97E-1, 2.85E-1, 2.69E-1, 2.51E-1, 2.30E-1, 2.07E-1, 1.81E-1},
287  { 0.50, 2.60E-1, 2.49E-1, 2.36E-1, 2.21E-1, 2.05E-1, 1.87E-1, 1.68E-1, 1.48E-1, 1.27E-1, 1.06E-1, 8.54E-2},
288  { 1.00, 1.86E-1, 1.72E-1, 1.57E-1, 1.41E-1, 1.25E-1, 1.09E-1, 9.28E-2, 7.71E-2, 6.20E-2, 4.79E-2, 3.50E-2},
289  { 1.50, 1.21E-1, 1.08E-1, 9.44E-2, 8.15E-2, 6.91E-2, 5.73E-2, 4.63E-2, 3.62E-2, 2.72E-2, 1.93E-2, 1.28E-2},
290  { 2.00, 7.20E-2, 6.19E-2, 5.22E-2, 4.33E-2, 3.51E-2, 2.79E-2, 2.11E-2, 1.55E-2, 1.09E-2, 7.11E-3, 4.22E-3},
291  { 2.50, 3.99E-2, 3.31E-2, 2.69E-2, 2.13E-2, 1.65E-2, 1.24E-2, 8.97E-3, 6.19E-3, 4.02E-3, 2.41E-3, 1.28E-3},
292  { 3.00, 2.07E-2, 1.66E-2, 1.30E-2, 9.87E-3, 7.30E-3, 5.21E-3, 3.56E-3, 2.31E-3, 1.39E-3, 7.61E-4, 3.60E-4},
293  { 3.50, 1.02E-2, 7.84E-3, 5.89E-3, 4.31E-3, 3.04E-3, 2.06E-3, 1.33E-3, 8.09E-4, 4.52E-4, 2.26E-4, 9.47E-5},
294  { 4.00, 4.74E-3, 3.52E-3, 2.55E-3, 1.78E-3, 1.20E-3, 7.76E-4, 4.73E-4, 2.69E-4, 1.39E-4, 6.32E-5, 2.34E-5},
295  { 4.50, 2.10E-3, 1.51E-3, 1.05E-3, 7.05E-4, 4.54E-4, 2.78E-4, 1.60E-4, 8.52E-5, 4.08E-5, 1.68E-5, 0},
296  { 5.00, 8.98E-4, 6.21E-4, 4.15E-4, 2.67E-4, 1.64E-4, 9.55E-5, 5.19E-5, 2.58E-5, 1.14E-5, 0, 0},
297  { 5.50, 3.68E-4, 2.45E-4, 1.58E-4, 9.71E-5, 5.70E-5, 3.15E-5, 1.61E-5, 0, 0, 0, 0},
298  { 6.00, 1.45E-4, 9.32E-5, 5.76E-5, 3.41E-5, 1.91E-5, 1.00E-5, 0, 0, 0, 0, 0},
299  { 6.50, 5.53E-5, 3.43E-5, 2.04E-5, 1.15E-5, 0, 0, 0, 0, 0, 0, 0},
300  { 7.00, 2.04E-5, 1.22E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
301 
302 
303 // kappa = 1.00
304 static double vavilovPdfValues6[19][12]={
305  {-3.50, 1.94E-5, 2.64E-5, 3.59E-5, 4.87E-5, 6.61E-5, 8.96E-5, 1.21E-4, 1.64E-4, 2.22E-4, 3.00E-4, 4.05E-4},
306  {-3.00, 1.83E-3, 2.37E-3, 3.06E-3, 3.94E-3, 5.08E-3, 6.53E-3, 8.39E-3, 1.08E-2, 1.38E-2, 1.76E-2, 2.25E-2},
307  {-2.50, 2.62E-2, 3.22E-2, 3.95E-2, 4.84E-2, 5.91E-2, 7.20E-2, 8.76E-2, 1.06E-1, 1.28E-1, 1.55E-1, 1.86E-1},
308  {-2.00, 1.19E-1, 1.39E-1, 1.63E-1, 1.89E-1, 2.18E-1, 2.51E-1, 2.88E-1, 3.29E-1, 3.74E-1, 4.23E-1, 4.75E-1},
309  {-1.50, 2.69E-1, 2.99E-1, 3.31E-1, 3.64E-1, 3.97E-1, 4.31E-1, 4.65E-1, 4.97E-1, 5.26E-1, 5.52E-1, 5.73E-1},
310  {-1.00, 3.85E-1, 4.07E-1, 4.26E-1, 4.43E-1, 4.56E-1, 4.66E-1, 4.69E-1, 4.67E-1, 4.58E-1, 4.42E-1, 4.17E-1},
311  {-0.50, 4.01E-1, 4.02E-1, 4.00E-1, 3.92E-1, 3.80E-1, 3.63E-1, 3.42E-1, 3.15E-1, 2.84E-1, 2.49E-1, 2.11E-1},
312  { 0.00, 3.29E-1, 3.13E-1, 2.94E-1, 2.73E-1, 2.49E-1, 2.22E-1, 1.94E-1, 1.65E-1, 1.36E-1, 1.08E-1, 8.10E-2},
313  { 0.50, 2.23E-1, 2.02E-1, 1.80E-1, 1.57E-1, 1.34E-1, 1.12E-1, 9.10E-2, 7.13E-2, 5.35E-2, 3.80E-2, 2.50E-2},
314  { 1.00, 1.29E-1, 1.11E-1, 9.38E-2, 7.73E-2, 6.21E-2, 4.84E-2, 3.64E-2, 2.61E-2, 1.78E-2, 1.12E-2, 6.42E-3},
315  { 1.50, 6.60E-2, 5.39E-2, 4.30E-2, 3.34E-2, 2.51E-2, 1.83E-2, 1.27E-2, 8.37E-3, 5.15E-3, 2.88E-3, 1.42E-3},
316  { 2.00, 3.01E-2, 2.33E-2, 1.76E-2, 1.29E-2, 9.10E-3, 6.15E-3, 3.95E-3, 2.38E-3, 1.32E-3, 6.54E-4, 2.74E-4},
317  { 2.50, 1.24E-2, 9.15E-3, 6.53E-3, 4.50E-3, 2.98E-3, 1.88E-3, 1.11E-3, 6.12E-4, 3.05E-4, 1.33E-4, 4.73E-5},
318  { 3.00, 4.71E-3, 3.29E-3, 2.22E-3, 1.44E-3, 8.94E-4, 5.24E-4, 2.87E-4, 1.44E-4, 6.44E-5, 2.46E-5, 0},
319  { 3.50, 1.65E-3, 1.09E-3, 6.99E-4, 4.28E-4, 2.48E-4, 1.35E-4, 6.81E-5, 3.11E-5, 1.55E-5, 0, 0},
320  { 4.00, 5.38E-4, 3.39E-4, 2.05E-4, 1.18E-4, 6.41E-5, 3.24E-5, 1.51E-5, 0, 0, 0, 0},
321  { 4.50, 1.65E-4, 9.86E-5, 5.64E-5, 3.05E-5, 1.55E-5, 0, 0, 0, 0, 0, 0},
322  { 5.00, 4.75E-5, 2.70E-5, 1.46E-5, 0, 0, 0, 0, 0, 0, 0, 0},
323  { 5.50, 1.30E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
324 
325 
326 // kappa = 4.00
327 static double vavilovPdfValues7[20][12]={
328  {-4.00, 0, 0, 0, 0, 0, 1.32E-5, 3.04E-5, 6.90E-5, 1.55E-4, 3.44E-4, 7.55E-4},
329  {-3.75, 1.38E-5, 2.98E-5, 6.38E-5, 1.35E-4, 2.83E-4, 5.85E-4, 1.20E-3, 2.41E-3, 4.78E-3, 9.35E-3, 1.80E-2},
330  {-3.50, 3.81E-4, 7.44E-4, 1.44E-3, 2.73E-3, 5.12E-3, 9.45E-3, 1.71E-2, 3.05E-2, 5.33E-2, 9.11E-2, 1.52E-1},
331  {-3.25, 4.78E-3, 8.44E-3, 1.47E-2, 2.50E-2, 4.19E-2, 6.88E-2, 1.10E-1, 1.73E-1, 2.64E-1, 3.90E-1, 5.57E-1},
332  {-3.00, 3.19E-2, 5.08E-2, 7.95E-2, 1.22E-1, 1.82E-1, 2.64E-1, 3.73E-1, 5.11E-1, 6.75E-1, 8.56E-1, 1.03E+0},
333  {-2.75, 1.27E-1, 1.83E-1, 2.56E-1, 3.51E-1, 4.66E-1, 6.00E-1, 7.44E-1, 8.86E-1, 1.01E+0, 1.08E+0, 1.09E+0},
334  {-2.50, 3.28E-1, 4.26E-1, 5.38E-1, 6.58E-1, 7.77E-1, 8.81E-1, 9.56E-1, 9.85E-1, 9.56E-1, 8.63E-1, 7.14E-1},
335  {-2.25, 5.89E-1, 6.91E-1, 7.84E-1, 8.57E-1, 8.97E-1, 8.95E-1, 8.46E-1, 7.51E-1, 6.18E-1, 4.65E-1, 3.12E-1},
336  {-2.00, 7.75E-1, 8.22E-1, 8.37E-1, 8.15E-1, 7.56E-1, 6.63E-1, 5.44E-1, 4.14E-1, 2.88E-1, 1.78E-1, 9.59E-2},
337  {-1.75, 7.79E-1, 7.45E-1, 6.81E-1, 5.91E-1, 4.85E-1, 3.73E-1, 2.65E-1, 1.72E-1, 1.01E-1, 5.10E-2, 2.17E-2},
338  {-1.50, 6.16E-1, 5.32E-1, 4.36E-1, 3.38E-1, 2.45E-1, 1.64E-1, 1.01E-1, 5.60E-2, 2.73E-2, 1.13E-2, 3.74E-3},
339  {-1.25, 3.95E-1, 3.07E-1, 2.26E-1, 1.56E-1, 9.96E-2, 5.85E-2, 3.11E-2, 1.46E-2, 5.90E-3, 1.97E-3, 5.05E-4},
340  {-1.00, 2.09E-1, 1.47E-1, 9.68E-2, 5.94E-2, 3.35E-2, 1.72E-2, 7.85E-3, 3.12E-3, 1.05E-3, 2.80E-4, 5.49E-5},
341  {-0.75, 9.33E-2, 5.91E-2, 3.49E-2, 1.91E-2, 9.47E-3, 4.23E-3, 1.66E-3, 5.59E-4, 1.54E-4, 3.29E-5, 0},
342  {-0.50, 3.56E-2, 2.04E-2, 1.08E-2, 5.22E-3, 2.29E-3, 8.90E-4, 3.00E-4, 8.49E-5, 1.93E-5, 0, 0},
343  {-0.25, 1.18E-2, 6.07E-3, 2.88E-3, 1.24E-3, 4.78E-4, 1.62E-4, 4.67E-5, 1.12E-5, 0, 0, 0},
344  { 0.00, 3.41E-3, 1.59E-3, 6.74E-4, 2.58E-4, 8.75E-5, 2.57E-5, 0, 0, 0, 0, 0},
345  { 0.25, 8.74E-4, 3.67E-4, 1.40E-4, 4.74E-5, 1.42E-5, 0, 0, 0, 0, 0, 0},
346  { 0.50, 2.00E-4, 7.57E-5, 2.58E-5, 0, 0, 0, 0, 0, 0, 0, 0},
347  { 0.75, 4.11E-5, 1.41E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
348 
349 
350 // kappa = 7.00
351 static double vavilovPdfValues8[21][12]={
352  {-4.40, 0, 0, 0, 0, 0, 0, 0, 0, 1.11E-5, 3.90E-5, 1.33E-4},
353  {-4.20, 0, 0, 0, 0, 0, 2.17E-5, 6.93E-5, 2.15E-4, 6.46E-4, 1.88E-3, 5.28E-3},
354  {-4.00, 0, 1.07E-5, 3.24E-5, 9.54E-5, 2.73E-4, 7.60E-4, 2.04E-3, 5.31E-3, 1.33E-2, 3.18E-2, 7.28E-2},
355  {-3.80, 1.11E-4, 2.99E-4, 7.80E-4, 1.97E-3, 4.82E-3, 1.14E-2, 2.57E-2, 5.57E-2, 1.15E-1, 2.24E-1, 4.12E-1},
356  {-3.60, 1.77E-3, 4.11E-3, 9.26E-3, 2.01E-2, 4.18E-2, 8.32E-2, 1.58E-1, 2.83E-1, 4.78E-1, 7.51E-1, 1.09E+0},
357  {-3.40, 1.54E-2, 3.10E-2, 6.01E-2, 1.12E-1, 1.97E-1, 3.31E-1, 5.23E-1, 7.74E-1, 1.06E+0, 1.33E+0, 1.50E+0},
358  {-3.20, 7.96E-2, 1.39E-1, 2.33E-1, 3.69E-1, 5.54E-1, 7.81E-1, 1.02E+0, 1.24E+0, 1.37E+0, 1.35E+0, 1.17E+0},
359  {-3.00, 2.63E-1, 3.99E-1, 5.74E-1, 7.78E-1, 9.89E-1, 1.17E+0, 1.27E+0, 1.25E+0, 1.10E+0, 8.46E-1, 5.51E-1},
360  {-2.80, 5.86E-1, 7.71E-1, 9.54E-1, 1.10E+0, 1.18E+0, 1.17E+0, 1.04E+0, 8.35E-1, 5.83E-1, 3.46E-1, 1.68E-1},
361  {-2.60, 9.21E-1, 1.05E+0, 1.12E+0, 1.10E+0, 9.97E-1, 8.19E-1, 6.02E-1, 3.88E-1, 2.14E-1, 9.71E-2, 3.44E-2},
362  {-2.40, 1.06E+0, 1.04E+0, 9.55E-1, 8.02E-1, 6.12E-1, 4.18E-1, 2.52E-1, 1.30E-1, 5.63E-2, 1.94E-2, 4.95E-3},
363  {-2.20, 9.18E-1, 7.85E-1, 6.16E-1, 4.40E-1, 2.83E-1, 1.61E-1, 7.89E-2, 3.27E-2, 1.10E-2, 2.84E-3, 5.18E-4},
364  {-2.00, 6.17E-1, 4.57E-1, 3.08E-1, 1.87E-1, 1.01E-1, 4.75E-2, 1.90E-2, 6.29E-3, 1.64E-3, 3.16E-4, 4.05E-5},
365  {-1.80, 3.28E-1, 2.10E-1, 1.22E-1, 6.30E-2, 2.85E-2, 1.11E-2, 3.62E-3, 9.50E-4, 1.91E-4, 2.78E-5, 0},
366  {-1.60, 1.41E-1, 7.84E-2, 3.89E-2, 1.71E-2, 6.49E-3, 2.09E-3, 5.52E-4, 1.15E-4, 1.77E-5, 0, 0},
367  {-1.40, 4.99E-2, 2.40E-2, 1.02E-2, 3.80E-3, 1.21E-3, 3.22E-4, 6.88E-5, 1.13E-5, 0, 0, 0},
368  {-1.20, 1.47E-2, 6.10E-3, 2.23E-3, 7.05E-4, 1.88E-4, 4.12E-5, 0, 0, 0, 0, 0},
369  {-1.00, 3.64E-3, 1.31E-3, 4.11E-4, 1.10E-4, 2.46E-5, 0, 0, 0, 0, 0, 0},
370  {-0.80, 7.71E-4, 2.40E-4, 6.46E-5, 1.47E-5, 0, 0, 0, 0, 0, 0, 0},
371  {-0.60, 1.41E-4, 3.80E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0},
372  {-0.40, 2.24E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
373 
374 
375 // kappa = 10.00
376 static double vavilovPdfValues9[21][12]={
377  {-4.60, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.85E-5, 1.79E-4},
378  {-4.40, 0, 0, 0, 0, 0, 1.18E-5, 5.14E-5, 2.12E-4, 8.32E-4, 3.08E-3, 1.07E-2},
379  {-4.20, 0, 0, 1.41E-5, 5.57E-5, 2.10E-4, 7.50E-4, 2.54E-3, 8.11E-3, 2.42E-2, 6.69E-2, 1.70E-1},
380  {-4.00, 5.27E-5, 1.84E-4, 6.15E-4, 1.95E-3, 5.83E-3, 1.64E-2, 4.32E-2, 1.06E-1, 2.37E-1, 4.82E-1, 8.79E-1},
381  {-3.80, 1.43E-3, 4.07E-3, 1.10E-2, 2.78E-2, 6.61E-2, 1.46E-1, 2.96E-1, 5.49E-1, 9.16E-1, 1.35E+0, 1.73E+0},
382  {-3.60, 1.79E-2, 4.17E-2, 9.09E-2, 1.85E-1, 3.46E-1, 5.96E-1, 9.30E-1, 1.30E+0, 1.59E+0, 1.68E+0, 1.47E+0},
383  {-3.40, 1.16E-1, 2.20E-1, 3.88E-1, 6.29E-1, 9.31E-1, 1.24E+0, 1.48E+0, 1.54E+0, 1.38E+0, 1.02E+0, 5.99E-1},
384  {-3.20, 4.21E-1, 6.50E-1, 9.23E-1, 1.19E+0, 1.39E+0, 1.44E+0, 1.30E+0, 1.01E+0, 6.46E-1, 3.31E-1, 1.28E-1},
385  {-3.00, 9.12E-1, 1.15E+0, 1.31E+0, 1.35E+0, 1.24E+0, 9.87E-1, 6.74E-1, 3.84E-1, 1.76E-1, 6.16E-2, 1.53E-2},
386  {-2.80, 1.25E+0, 1.28E+0, 1.18E+0, 9.66E-1, 6.92E-1, 4.25E-1, 2.18E-1, 9.11E-2, 2.95E-2, 6.96E-3, 1.09E-3},
387  {-2.60, 1.13E+0, 9.45E-1, 7.01E-1, 4.56E-1, 2.55E-1, 1.20E-1, 4.64E-2, 1.41E-2, 3.19E-3, 5.02E-4, 4.87E-5},
388  {-2.40, 7.06E-1, 4.80E-1, 2.87E-1, 1.48E-1, 6.46E-2, 2.33E-2, 6.71E-3, 1.47E-3, 2.33E-4, 2.41E-5, 0},
389  {-2.20, 3.13E-1, 1.73E-1, 8.33E-2, 3.41E-2, 1.16E-2, 3.19E-3, 6.84E-4, 1.08E-4, 1.18E-5, 0, 0},
390  {-2.00, 1.02E-1, 4.59E-2, 1.77E-2, 5.74E-3, 1.52E-3, 3.19E-4, 5.06E-5, 0, 0, 0, 0},
391  {-1.80, 2.48E-2, 9.10E-3, 2.82E-3, 7.24E-4, 1.49E-4, 2.37E-5, 0, 0, 0, 0, 0},
392  {-1.60, 4.64E-3, 1.38E-3, 3.45E-4, 6.99E-5, 1.11E-5, 0, 0, 0, 0, 0, 0},
393  {-1.40, 6.77E-4, 1.64E-4, 3.29E-5, 0, 0, 0, 0, 0, 0, 0, 0},
394  {-1.20, 7.84E-5, 1.55E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
395 
397 
398 static double vavilovKappaValues[10] = {.01, .04, .07, .1, .4, .7, 1, 4, 7, 10};
399 
400 static int vavilovNLambda[10] = {45, 42, 41, 41, 28, 22, 19, 20, 21, 21};
401 
403  return 10;
404 }
405 
406 double VavilovTest::GetSBKappa (int ikappa) {
407  if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return -1;
408  return vavilovKappaValues[ikappa];
409 }
410 
412  return 11;
413 }
414 
415 double VavilovTest::GetSBBeta2 (int ibeta2) {
416  if (ibeta2 < 0 || ibeta2 >= VavilovTest::GetSBNBeta2()) return -1;
417  return 0.1*ibeta2;
418 }
419 
420 int VavilovTest::GetSBNLambda (int ikappa) {
421  if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return -1;
422  return vavilovNLambda[ikappa];
423 }
424 
425 double VavilovTest::GetSBLambda (int ikappa, int ilambda) {
426  if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return 0;
427  if (ilambda < 0 || ilambda >= VavilovTest::GetSBNLambda(ikappa)) return 0;
428  return vavilovPdfValues[ikappa][ilambda][0];
429 }
430 
431 double VavilovTest::GetSBVavilovPdf (int ikappa, int ibeta2, int ilambda) {
432  if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return 0;
433  if (ibeta2 < 0 || ibeta2 >= VavilovTest::GetSBNBeta2()) return 0;
434  if (ilambda < 0 || ilambda >= VavilovTest::GetSBNLambda(ikappa)) return 0;
435  return vavilovPdfValues[ikappa][ilambda][ibeta2+1];
436 }
437 
438 static double myRound (double x, double y, double& xmantissa, int digits) {
439  int exponent;
440  if (y) exponent = std::floor(std::log10(fabs(y)));
441  else if (x) exponent = std::floor(std::log10(fabs(x)));
442  else exponent = 0;
443  double power = std::pow (10.0, exponent);
444  double mantissa = y/power;
445  double dpower = std::pow (10.0, digits-1);
446  mantissa = roundl (mantissa*dpower)/dpower;
447  if (mantissa >= 10) {
448  mantissa *= 0.1;
449  exponent += 1;
450  power = std::pow (10.0, exponent);
451  }
452  xmantissa = x/power;
453  mantissa = roundl (xmantissa*dpower)/dpower;
454  return mantissa*power;
455 }
456 double myRound (double x, double y, int digits) {
457  double xmantissa;
458  return myRound (x, y, xmantissa, digits);
459 
460 }
461 
462 static std::string format (double x, double y, int digits, int width) {
463  int exponent;
464  if (y) exponent = std::floor(std::log10(fabs(y)));
465  else if (x) exponent = std::floor(std::log10(fabs(x)));
466  else exponent = 0;
467  double power = std::pow (10.0, exponent);
468  double mantissa = y/power;
469  double dpower = std::pow (10.0, digits-1);
470  mantissa = roundl (mantissa*dpower)/dpower;
471  if (mantissa >= 10) {
472  mantissa *= 0.1;
473  exponent += 1;
474  }
475  mantissa = roundl (x/power*dpower)/dpower;
476 
477  std::stringstream out;
478  out << std::setw(width-4) << std::fixed << std::setprecision(digits-1) << mantissa;
479  out << "e" << std::showpos << std::setw(3)<< std::internal << std::setfill('0') << exponent;
480  return out.str();
481 
482 }
483 
484 int VavilovTest::PdfTest (ROOT::Math::Vavilov& v, std::ostream& os) {
485  double maxabsdiff, maxdiffmantissa, agreefraction, agreediffmantissa;
486  GetPdfTestParams (v, maxabsdiff, maxdiffmantissa, agreefraction, agreediffmantissa);
487  return VavilovTest::PdfTest (v, os, maxabsdiff, maxdiffmantissa, agreefraction, agreediffmantissa);
488 }
489 
491  double maxabsdiff, double maxdiffmantissa,
492  double agreefraction, double agreediffmantissa) {
493 
494  std::ios::fmtflags defaultflags = os.flags();
495  int defaultwidth = os.width();
496  int defaultprecision = os.precision();
497 
498  int nfail = 0;
499 
500  os << "\n\nTesting Pdf\n\n";
501 
502  for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
503  double kappa = GetSBKappa (ikappa);
504  os << "\n\nkappa = " << kappa << "\n\n";
505 
506  double maxreldiff = 0;
507  double maxmandiff = 0;
508  bool pass = true;
509  int agree = 0;
510  int disagree = 0;
511 
512  for (int ilambda = 0; ilambda < GetSBNLambda (ikappa); ++ilambda) {
513  double lambda = GetSBLambda (ikappa, ilambda);
514  os << std::setw(5) << std::fixed << std::setprecision(2) << lambda;
515  double x = lambda;
516  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
517  double beta2 = GetSBBeta2 (ibeta);
518  v.SetKappaBeta2 (kappa, beta2);
519  double pdf = v.Pdf(x);
520  double val = GetSBVavilovPdf (ikappa, ibeta, ilambda);
521  double diffmantissa;
522  myRound (pdf-val, val, diffmantissa, 3);
523 
524  if (val > 0 && pdf > 0) {
525  double absdiff = fabs(pdf-val);
526  if (absdiff > maxabsdiff && fabs(diffmantissa) > maxdiffmantissa) {
527  //os << "FAIL";
528  pass = false;
529  }
530  if (fabs(diffmantissa) > 0.006) {
531  if (absdiff > maxabsdiff) maxabsdiff = absdiff;
532  if (absdiff/val > maxreldiff) maxreldiff = absdiff/val;
533  if (fabs(diffmantissa) > maxmandiff) maxmandiff = fabs(diffmantissa);
534  }
535  if (fabs(diffmantissa) > agreediffmantissa) ++disagree; else ++agree;
536  }
537  if (val == 0)
538  os << " ";
539  else if (pdf == 0)
540  os << " --- ";
541  else if (fabs(diffmantissa) > 0.006)
542  os << format (pdf-val, val, 3, 10);
543  else
544  os << " 0 ";
545  }
546  os << "\n";
547  }
548  os.flags (defaultflags);
549  if (agree < disagree*agreefraction) {
550  pass = false;
551  //os << "agreefraction test failed.\n";
552  }
553  if (!pass) ++nfail;
554  os << "Max abs diff: " << maxabsdiff << ", max rel diff: " << maxreldiff
555  << ", max diff mantissa: " << std::fixed << std::setprecision(2) << maxmandiff
556  << ", agree/disagree=" << agree << "/" << disagree
557  << ", pass=" << pass << std::endl;
558  }
559  os << "\n\nNumber of failed tests: " << nfail << std::endl;
560 
561  os.flags (defaultflags);
562  os.width(defaultwidth);
563  os.precision(defaultprecision);
564 
565  return nfail;
566 }
567 
568 static void moments (ROOT::Math::Vavilov& v, double& integral,
569  double& mean, double& variance,
570  double& skewness, double& kurtosis) {
571  int nsteps = 10000;
572  double t0 = v.GetLambdaMin();
573  double t1 = v.GetLambdaMax();
574  double t = t1 - t0;
575  double dt = t/nsteps;
576 
577  double sum = 0;
578  double sumx = 0;
579  for (int i = 0; i < nsteps; ++i) {
580  double x = (i+0.5)*dt + t0;
581  double pdf = v.Pdf(x);
582  sum += pdf;
583  sumx += x*pdf;
584  }
585  integral = sum*dt;
586  mean = sumx/sum;
587  double sumx2 = 0;
588  double sumx3 = 0;
589  double sumx4 = 0;
590  for (int i = 0; i < nsteps; ++i) {
591  double x = (i+0.5)*dt + t0;
592  double pdf = v.Pdf(x);
593  sumx2 += pow(x-mean, 2)*pdf;
594  sumx3 += pow(x-mean, 3)*pdf;
595  sumx4 += pow(x-mean, 4)*pdf;
596  }
597  variance = sumx2/sum;
598  skewness = sumx3/sum*pow (variance, -1.5);
599  kurtosis = sumx4/sum/(variance*variance)-3;
600 }
601 
602 void VavilovTest::PrintPdfTable (ROOT::Math::Vavilov& v, std::ostream& os, int digits) {
603  std::ios::fmtflags defaultflags = os.flags();
604  int defaultwidth = os.width();
605  int defaultprecision = os.precision();
606 
607  for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
608  double kappa = GetSBKappa (ikappa);
609  os << "\n\nkappa = " << kappa << "\n\n";
610 
611  for (int ilambda = 0; ilambda < GetSBNLambda (ikappa); ++ilambda) {
612  double lambda = GetSBLambda (ikappa, ilambda);
613  os << std::setw(5) << std::fixed << std::setprecision(2) << lambda;
614  double x = lambda;
615  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
616  double beta2 = GetSBBeta2 (ibeta);
617  v.SetKappaBeta2 (kappa, beta2);
618  double pdf = v.Pdf(x);
619  if (pdf > 0)
620  os << std::setw(digits+7) << std::scientific << std::setprecision(digits-1) << pdf;
621  else
622  os << std::setw(digits+7) << " ";
623  }
624  os << "\n";
625  }
626  os << "Xmin:";
627  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
628  double beta2 = GetSBBeta2 (ibeta);
629  v.SetKappaBeta2 (kappa, beta2);
630  os << std::setw(digits+7) << std::fixed << std::setprecision(digits-1) << v.GetLambdaMin();
631  }
632  os << "\n";
633  os << "Xmax:";
634  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
635  double beta2 = GetSBBeta2 (ibeta);
636  v.SetKappaBeta2 (kappa, beta2);
637  os << std::setw(digits+7) << std::fixed << std::setprecision(digits-1) << v.GetLambdaMax();
638  }
639  os << "\n";
640  os << "int: ";
641  double integral[11], calcmean[11], calcvariance[11], calcskewness[11], calckurtosis[11];
642  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
643  double beta2 = GetSBBeta2 (ibeta);
644  v.SetKappaBeta2 (kappa, beta2);
645  moments (v, integral[ibeta], calcmean[ibeta], calcvariance[ibeta], calcskewness[ibeta], calckurtosis[ibeta]);
646  os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << integral[ibeta];
647  }
648  os << "\nmean:";
649  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
650  os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calcmean[ibeta];
651  }
652  os << "\n ";
653  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
654  double beta2 = GetSBBeta2 (ibeta);
655  os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Mean (kappa, beta2);
656  }
657  os << "\nvar: ";
658  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
659  os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calcvariance[ibeta];
660  }
661  os << "\n ";
662  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
663  double beta2 = GetSBBeta2 (ibeta);
664  os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Variance (kappa, beta2);
665  }
666  os << "\nskew:";
667  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
668  os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calcskewness[ibeta];
669  }
670  os << "\n ";
671  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
672  double beta2 = GetSBBeta2 (ibeta);
673  os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Skewness (kappa, beta2);
674  }
675  os << "\nkurt:";
676  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
677  os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calckurtosis[ibeta];
678  }
679  os << "\n ";
680  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
681  double beta2 = GetSBBeta2 (ibeta);
682  os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Kurtosis (kappa, beta2);
683  }
684  os << "\n";
685  }
686 
687  os.flags (defaultflags);
688  os.width(defaultwidth);
689  os.precision(defaultprecision);
690 }
691 
692 int VavilovTest::CdfTest (ROOT::Math::Vavilov& v, std::ostream& os) {
693  double maxabsdiff, maxcdfdiff;
694  GetCdfTestParams (v, maxabsdiff, maxcdfdiff);
695  return VavilovTest::CdfTest (v, os, maxabsdiff, maxcdfdiff);
696 }
697 
698 int VavilovTest::CdfTest (ROOT::Math::Vavilov& v, std::ostream& os, double maxabsdiff, double maxcdfdiff) {
699  std::ios::fmtflags defaultflags = os.flags();
700  int defaultwidth = os.width();
701  int defaultprecision = os.precision();
702 
703  int nfail = 0;
704 
705  os << "\n\nTesting Cdf and Cdf_c\n\n";
706 
707  for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
708  double kappa = GetSBKappa (ikappa);
709  os << "\n\nkappa = " << kappa << "\n\n";
710 
711  double absdiffmax = 0;
712  double reldiffmax = 0;
713  double cdfdiffmax = 0;
714  bool pass = true;
715 
716  double cdf_calc[11];
717  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
718  cdf_calc[ibeta] = 0;
719  }
720 
721  for (int ilambda = 0; ilambda < GetSBNLambda (ikappa); ++ilambda) {
722  double lambda1 = GetSBLambda (ikappa, ilambda);
723  os << std::setw(5) << std::fixed << std::setprecision(2) << lambda1;
724  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
725  double beta2 = GetSBBeta2 (ibeta);
726  v.SetKappaBeta2 (kappa, beta2);
727  double lambda0 = (ilambda > 0) ? GetSBLambda (ikappa, ilambda-1) : v.GetLambdaMin();
728  if (lambda1 > lambda0) {
729  int n = 100;
730  double dlambda = (lambda1 - lambda0)/n;
731  for (int i = 0; i < n; ++i) {
732  double lambda = lambda0 + (i+0.5)*dlambda;
733  cdf_calc[ibeta] += v.Pdf(lambda)*dlambda;
734  }
735  }
736 
737  double cdf = v.Cdf(lambda1);
738  double cdf_c = v.Cdf_c(lambda1);
739  double val = cdf_calc[ibeta];
740 
741  if (fabs(cdf-val) > absdiffmax) absdiffmax = fabs(cdf-val);
742  if (fabs(cdf+cdf_c-1) > cdfdiffmax) cdfdiffmax = cdf+cdf_c-1;
743  if (val > 0 && fabs(cdf/val-1) > reldiffmax) reldiffmax = fabs(cdf/val-1);
744 
745  if (val == 0)
746  os << " ";
747  else
748  os << std::scientific << std::setw(10) << std::setprecision(2) << cdf-val;
749  }
750  os << "\n";
751  }
752  os.flags (defaultflags);
753  if (absdiffmax > maxabsdiff) pass = false;
754  if (cdfdiffmax > maxcdfdiff) pass = false;
755  if (!pass) ++nfail;
756  os << "Max abs diff: " << absdiffmax << ", max rel diff: " << reldiffmax
757  << ", max diff cdf+cdf_c-1: " << cdfdiffmax
758  << ", pass=" << pass << std::endl;
759  }
760  os << "\n\nNumber of failed tests: " << nfail << std::endl;
761 
762  os.flags (defaultflags);
763  os.width(defaultwidth);
764  os.precision(defaultprecision);
765 
766  return nfail;
767 }
768 
770  double maxabsdiff;
771  GetQuantileTestParams (v, maxabsdiff);
772  return VavilovTest::QuantileTest (v, os, maxabsdiff);
773 }
774 
775 int VavilovTest::QuantileTest (ROOT::Math::Vavilov& v, std::ostream& os, double maxabsdiff) {
776  double qmin = 0;
777 
778  std::ios::fmtflags defaultflags = os.flags();
779  int defaultwidth = os.width();
780  int defaultprecision = os.precision();
781 
782  int nfail = 0;
783 
784  static const double qvalues[45] = {0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009,
785  0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09,
786  0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
787  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99,
788  0.991, 0.992, 0.993, 0.994, 0.995, 0.996, 0.997, 0.998, 0.999};
789 
790  os << "\n\nTesting Quantile\n\n";
791 
792  for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
793  double kappa = GetSBKappa (ikappa);
794  os << "\n\nkappa = " << kappa << "\n\n";
795 
796  double absdiffmax = 0;
797  bool pass = true;
798 
799  for (int iq = 0; iq < 45; ++iq) {
800  double qval = qvalues[iq];
801  os << std::setw(5) << std::fixed << std::setprecision(3) << qval;
802  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
803  double beta2 = GetSBBeta2 (ibeta);
804  v.SetKappaBeta2 (kappa, beta2);
805  double lambda = v.Quantile (qval);
806  double cdfval = v.Cdf(lambda);
807  if (qval > qmin && (1-qval) > qmin) {
808  if (fabs(cdfval-qval) > absdiffmax) absdiffmax = fabs(cdfval-qval);
809  os << " " << std::fixed << std::setw(7) << std::setprecision(4) << cdfval-qval << " ";
810  }
811  else {
812  os << " (" << std::fixed << std::setw(7) << std::setprecision(4) << cdfval-qval << ")";
813  }
814  }
815  os << "\n";
816  }
817  os.flags (defaultflags);
818  if (absdiffmax > maxabsdiff) pass = false;
819  if (!pass) ++nfail;
820  os << "Max abs diff: " << absdiffmax
821  << ", pass=" << pass << std::endl;
822  }
823 
824  os << "\n\nTesting Quantile_c\n\n";
825 
826  for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
827  double kappa = GetSBKappa (ikappa);
828  os << "\n\nkappa = " << kappa << "\n\n";
829 
830  double absdiffmax = 0;
831  bool pass = true;
832 
833  for (int iq = 0; iq < 45; ++iq) {
834  double qval = qvalues[iq];
835  os << std::setw(5) << std::fixed << std::setprecision(3) << qval;
836  for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
837  double beta2 = GetSBBeta2 (ibeta);
838  v.SetKappaBeta2 (kappa, beta2);
839  double lambda_c = v.Quantile_c (qval);
840  double cdf_c_val = v.Cdf_c(lambda_c);
841  if (qval > qmin && (1-qval) > qmin) {
842  if (fabs(cdf_c_val-qval) > absdiffmax) absdiffmax = fabs(cdf_c_val-qval);
843 
844  os << " " << std::fixed << std::setw(7) << std::setprecision(4) << cdf_c_val-qval << " ";
845  }
846  else {
847  os << " (" << std::fixed << std::setw(7) << std::setprecision(4) << cdf_c_val-qval << ")";
848  }
849  }
850  os << "\n";
851  }
852  os.flags (defaultflags);
853  if (absdiffmax > maxabsdiff) pass = false;
854  if (!pass) ++nfail;
855  os << "Max abs diff: " << absdiffmax
856  << ", pass=" << pass << std::endl;
857  }
858  os << "\n\nNumber of failed tests: " << nfail << std::endl;
859 
860  os.flags (defaultflags);
861  os.width(defaultwidth);
862  os.precision(defaultprecision);
863 
864  return nfail;
865 }
866 
867 void VavilovTest::GetPdfTestParams (const Vavilov& v, double& maxabsdiff, double& maxdiffmantissa, double& agreefraction, double& agreediffmantissa) {
868  if (dynamic_cast <const VavilovFast *>(&v)) {
869  maxabsdiff = 0.08;
870  maxdiffmantissa = 0.1;
871  agreefraction = 1;
872  agreediffmantissa = 0.9;
873  }
874  else {
875  maxabsdiff = 2E-3;
876  maxdiffmantissa = 0.03;
877  agreefraction = 5;
878  agreediffmantissa = 0.015;
879  }
880 }
881 
882 void VavilovTest::GetCdfTestParams (const Vavilov& v, double& maxabsdiff, double& maxcdfdiff) {
883  if (dynamic_cast <const VavilovFast *>(&v)) {
884  maxabsdiff = 0.018;
885  maxcdfdiff = 1E-16;
886  }
887  else {
888  maxabsdiff = 1E-5;
889  maxcdfdiff = 5E-15;
890  }
891 }
892 
893 void VavilovTest::GetQuantileTestParams (const Vavilov& v, double& maxabsdiff) {
894  if (dynamic_cast <const VavilovFast *>(&v)) {
895  maxabsdiff = 0.03;
896  }
897  else {
898  maxabsdiff = 2E-10;
899  }
900 }
901 
902 
903 } // namespace Math
904 } // namespace ROOT
Base class describing a Vavilov distribution.
Definition: Vavilov.h:120
static long int sum(long int i)
Definition: Factory.cxx:2162
static void GetCdfTestParams(const Vavilov &v, double &maxabsdiff, double &maxcdfdiff)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
static double vavilovPdfValues7[20][12]
static int PdfTest(Vavilov &v, std::ostream &os, double maxabsdiff, double maxdiffmantissa, double agreefraction, double agreediffmantissa)
Test the pdf values against the tables of Seltzer and Berger.
static int vavilovNLambda[10]
static double vavilovPdfValues6[19][12]
virtual double Variance() const
Return the theoretical variance .
Definition: Vavilov.cxx:88
static int GetSBNKappa()
Return the number of values in the tables of Seltzer and Berger.
static double(*[10] vavilovPdfValues)[12]
static double vavilovPdfValues8[21][12]
virtual double Cdf_c(double x) const =0
Evaluate the Vavilov complementary cumulative probability density function.
static double GetSBLambda(int ikappa, int ilambda)
Return the value in the tables of Seltzer and Berger.
static double vavilovKappaValues[10]
static double vavilovPdfValues1[42][12]
static std::string format(double x, double y, int digits, int width)
virtual double Pdf(double x) const =0
Evaluate the Vavilov probability density function.
static double vavilovPdfValues4[28][12]
TLatex * t1
Definition: textangle.C:20
virtual double GetLambdaMin() const =0
Return the minimum value of for which is nonzero in the current approximation.
#define roundl(x)
Definition: VavilovTest.cxx:51
static void PrintPdfTable(Vavilov &v, std::ostream &os, int digits=3)
Print a table of the pdf values to stream os.
static double myRound(double x, double y, double &xmantissa, int digits)
Double_t x[n]
Definition: legend1.C:17
double cdf(double *x, double *p)
Definition: unuranDistr.cxx:44
double log10(double)
static double vavilovPdfValues2[41][12]
double pow(double, double)
static void GetPdfTestParams(const Vavilov &v, double &maxabsdiff, double &maxdiffmantissa, double &agreefraction, double &agreediffmantissa)
static double GetSBVavilovPdf(int ikappa, int ibeta2, int ilambda)
Return the value of in the tables of Seltzer and Berger.
virtual double Kurtosis() const
Return the theoretical kurtosis .
Definition: Vavilov.cxx:105
static int GetSBNBeta2()
Return the number of values in the tables of Seltzer and Berger.
virtual double Cdf(double x) const =0
Evaluate the Vavilov cumulative probability density function.
virtual double Quantile(double z) const =0
Evaluate the inverse of the Vavilov cumulative probability density function.
static double vavilovPdfValues3[41][12]
static void GetQuantileTestParams(const Vavilov &v, double &maxabsdiff)
virtual double Skewness() const
Return the theoretical skewness .
Definition: Vavilov.cxx:96
static double vavilovPdfValues0[45][12]
Definition: VavilovTest.cxx:61
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
SVector< double, 2 > v
Definition: Dict.h:5
static int QuantileTest(Vavilov &v, std::ostream &os, double maxabsdiff)
Test the quantile values against the cdf Returns 0 if the test is passed.
double floor(double)
static void moments(ROOT::Math::Vavilov &v, double &integral, double &mean, double &variance, double &skewness, double &kurtosis)
static int CdfTest(Vavilov &v, std::ostream &os, double maxabsdiff, double maxcdfdiff)
Test the cdf values against the integral of the pdf Returns 0 if the test is passed.
constexpr Double_t E()
Definition: TMath.h:74
virtual double Quantile_c(double z) const =0
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
Double_t y[n]
Definition: legend1.C:17
static const double eu
Definition: Vavilov.cxx:44
Namespace for new Math classes and functions.
static double GetSBKappa(int ikappa)
Return the value for ikappa in the tables of Seltzer and Berger.
static double vavilovPdfValues9[21][12]
virtual double GetLambdaMax() const =0
Return the maximum value of for which is nonzero in the current approximation.
static int GetSBNLambda(int ikappa)
Return the number of values in the tables of Seltzer and Berger.
static double GetSBBeta2(int ibeta2)
Return the value in the tables of Seltzer and Berger.
int * iq
Definition: THbookFile.cxx:86
virtual void SetKappaBeta2(double kappa, double beta2)=0
Change and and recalculate coefficients if necessary.
const Int_t n
Definition: legend1.C:16
virtual double Mean() const
Return the theoretical mean , where is Euler&#39;s constant.
Definition: Vavilov.cxx:80
static double vavilovPdfValues5[22][12]