// @(#)root/mathmore:$Id$
// Authors: B. List 29.4.2010
/**********************************************************************
* *
* Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU General Public License *
* as published by the Free Software Foundation; either version 2 *
* of the License, or (at your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this library (see file COPYING); if not, write *
* to the Free Software Foundation, Inc., 59 Temple Place, Suite *
* 330, Boston, MA 02111-1307 USA, or contact the author. *
* *
**********************************************************************/
// Implementation file for class VavilovTest
//
// Created by: blist at Thu Apr 29 11:19:00 2010
//
// Last update: Thu Apr 29 11:19:00 2010
//
#include "VavilovTest.h"
#include "Math/Vavilov.h"
#include "Math/VavilovAccurate.h"
#include "Math/VavilovFast.h"
//#include "Math/SpecFuncMathCore.h"
//#include "Math/SpecFuncMathMore.h"
#include <cassert>
#include <iostream>
#include <cmath>
#include <iomanip>
#include <cmath>
#include <cstdlib>
#include <string>
#include <sstream>
#ifndef roundl
#define roundl(x) (long double)((long long)((x == 0) ? 0.0 : ( (x) + ( ((x) > 0) ? 0.5 : -0.5) )))
#endif
namespace ROOT {
namespace Math {
static const double eu = 0.577215664901532860606; // Euler's constant
//kappa = 0.01
static double vavilovPdfValues0[45][12]={
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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}};
//kappa = 0.04
static double vavilovPdfValues1[42][12]={
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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}};
//kappa = 0.07
static double vavilovPdfValues2[41][12]={
{-3.50, 0, 0, 0, 0, 0, 0, 1.01E-5, 1.06E-5, 1.10E-5, 1.14E-5, 1.19E-5},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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}};
//kappa = 0.10
static double vavilovPdfValues3[41][12]={
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{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},
{24.00, 3.26E-5, 2.59E-5, 2.01E-5, 1.53E-5, 1.12E-5, 0, 0, 0, 0, 0, 0},
{26.00, 1.33E-5, 1.04E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
// kappa = 0.40
static double vavilovPdfValues4[28][12]={
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 9.50, 6.08E-5, 4.19E-5, 2.79E-5, 1.79E-5, 1.09E-5, 0, 0, 0, 0, 0, 0},
{10.00, 3.17E-5, 2.13E-5, 1.39E-5, 0, 0, 0, 0, 0, 0, 0, 0}};
//kappa = 0.70
static double vavilovPdfValues5[22][12]={
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 6.50, 5.53E-5, 3.43E-5, 2.04E-5, 1.15E-5, 0, 0, 0, 0, 0, 0, 0},
{ 7.00, 2.04E-5, 1.22E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
// kappa = 1.00
static double vavilovPdfValues6[19][12]={
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 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},
{ 4.50, 1.65E-4, 9.86E-5, 5.64E-5, 3.05E-5, 1.55E-5, 0, 0, 0, 0, 0, 0},
{ 5.00, 4.75E-5, 2.70E-5, 1.46E-5, 0, 0, 0, 0, 0, 0, 0, 0},
{ 5.50, 1.30E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
// kappa = 4.00
static double vavilovPdfValues7[20][12]={
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{ 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},
{ 0.25, 8.74E-4, 3.67E-4, 1.40E-4, 4.74E-5, 1.42E-5, 0, 0, 0, 0, 0, 0},
{ 0.50, 2.00E-4, 7.57E-5, 2.58E-5, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0.75, 4.11E-5, 1.41E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
// kappa = 7.00
static double vavilovPdfValues8[21][12]={
{-4.40, 0, 0, 0, 0, 0, 0, 0, 0, 1.11E-5, 3.90E-5, 1.33E-4},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-1.00, 3.64E-3, 1.31E-3, 4.11E-4, 1.10E-4, 2.46E-5, 0, 0, 0, 0, 0, 0},
{-0.80, 7.71E-4, 2.40E-4, 6.46E-5, 1.47E-5, 0, 0, 0, 0, 0, 0, 0},
{-0.60, 1.41E-4, 3.80E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-0.40, 2.24E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
// kappa = 10.00
static double vavilovPdfValues9[21][12]={
{-4.60, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.85E-5, 1.79E-4},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-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},
{-1.60, 4.64E-3, 1.38E-3, 3.45E-4, 6.99E-5, 1.11E-5, 0, 0, 0, 0, 0, 0},
{-1.40, 6.77E-4, 1.64E-4, 3.29E-5, 0, 0, 0, 0, 0, 0, 0, 0},
{-1.20, 7.84E-5, 1.55E-5, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
static double (*vavilovPdfValues[10])[12]={vavilovPdfValues0, vavilovPdfValues1, vavilovPdfValues2, vavilovPdfValues3, vavilovPdfValues4, vavilovPdfValues5, vavilovPdfValues6, vavilovPdfValues7, vavilovPdfValues8, vavilovPdfValues9};
static double vavilovKappaValues[10] = {.01, .04, .07, .1, .4, .7, 1, 4, 7, 10};
static int vavilovNLambda[10] = {45, 42, 41, 41, 28, 22, 19, 20, 21, 21};
int VavilovTest::GetSBNKappa () {
return 10;
}
double VavilovTest::GetSBKappa (int ikappa) {
if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return -1;
return vavilovKappaValues[ikappa];
}
int VavilovTest::GetSBNBeta2 () {
return 11;
}
double VavilovTest::GetSBBeta2 (int ibeta2) {
if (ibeta2 < 0 || ibeta2 >= VavilovTest::GetSBNBeta2()) return -1;
return 0.1*ibeta2;
}
int VavilovTest::GetSBNLambda (int ikappa) {
if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return -1;
return vavilovNLambda[ikappa];
}
double VavilovTest::GetSBLambda (int ikappa, int ilambda) {
if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return 0;
if (ilambda < 0 || ilambda >= VavilovTest::GetSBNLambda(ikappa)) return 0;
return vavilovPdfValues[ikappa][ilambda][0];
}
double VavilovTest::GetSBVavilovPdf (int ikappa, int ibeta2, int ilambda) {
if (ikappa < 0 || ikappa >= VavilovTest::GetSBNKappa()) return 0;
if (ibeta2 < 0 || ibeta2 >= VavilovTest::GetSBNBeta2()) return 0;
if (ilambda < 0 || ilambda >= VavilovTest::GetSBNLambda(ikappa)) return 0;
return vavilovPdfValues[ikappa][ilambda][ibeta2+1];
}
static double myRound (double x, double y, double& xmantissa, int digits) {
int exponent;
if (y) exponent = std::floor(std::log10(fabs(y)));
else if (x) exponent = std::floor(std::log10(fabs(x)));
else exponent = 0;
double power = std::pow (10.0, exponent);
double mantissa = y/power;
double dpower = std::pow (10.0, digits-1);
mantissa = roundl (mantissa*dpower)/dpower;
if (mantissa >= 10) {
mantissa *= 0.1;
exponent += 1;
power = std::pow (10.0, exponent);
}
xmantissa = x/power;
mantissa = roundl (xmantissa*dpower)/dpower;
return mantissa*power;
}
double myRound (double x, double y, int digits) {
double xmantissa;
return myRound (x, y, xmantissa, digits);
}
static std::string format (double x, double y, int digits, int width) {
int exponent;
if (y) exponent = std::floor(std::log10(fabs(y)));
else if (x) exponent = std::floor(std::log10(fabs(x)));
else exponent = 0;
double power = std::pow (10.0, exponent);
double mantissa = y/power;
double dpower = std::pow (10.0, digits-1);
mantissa = roundl (mantissa*dpower)/dpower;
if (mantissa >= 10) {
mantissa *= 0.1;
exponent += 1;
}
mantissa = roundl (x/power*dpower)/dpower;
std::stringstream out;
out << std::setw(width-4) << std::fixed << std::setprecision(digits-1) << mantissa;
out << "e" << std::showpos << std::setw(3)<< std::internal << std::setfill('0') << exponent;
return out.str();
}
int VavilovTest::PdfTest (ROOT::Math::Vavilov& v, std::ostream& os) {
double maxabsdiff, maxdiffmantissa, agreefraction, agreediffmantissa;
GetPdfTestParams (v, maxabsdiff, maxdiffmantissa, agreefraction, agreediffmantissa);
return VavilovTest::PdfTest (v, os, maxabsdiff, maxdiffmantissa, agreefraction, agreediffmantissa);
}
int VavilovTest::PdfTest (ROOT::Math::Vavilov& v, std::ostream& os,
double maxabsdiff, double maxdiffmantissa,
double agreefraction, double agreediffmantissa) {
std::ios::fmtflags defaultflags = os.flags();
int defaultwidth = os.width();
int defaultprecision = os.precision();
int nfail = 0;
os << "\n\nTesting Pdf\n\n";
for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
double kappa = GetSBKappa (ikappa);
os << "\n\nkappa = " << kappa << "\n\n";
double maxreldiff = 0;
double maxmandiff = 0;
bool pass = true;
int agree = 0;
int disagree = 0;
for (int ilambda = 0; ilambda < GetSBNLambda (ikappa); ++ilambda) {
double lambda = GetSBLambda (ikappa, ilambda);
os << std::setw(5) << std::fixed << std::setprecision(2) << lambda;
double x = lambda;
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
v.SetKappaBeta2 (kappa, beta2);
double pdf = v.Pdf(x);
double val = GetSBVavilovPdf (ikappa, ibeta, ilambda);
double diffmantissa;
myRound (pdf-val, val, diffmantissa, 3);
if (val > 0 && pdf > 0) {
double absdiff = fabs(pdf-val);
if (absdiff > maxabsdiff && fabs(diffmantissa) > maxdiffmantissa) {
//os << "FAIL";
pass = false;
}
if (fabs(diffmantissa) > 0.006) {
if (absdiff > maxabsdiff) maxabsdiff = absdiff;
if (absdiff/val > maxreldiff) maxreldiff = absdiff/val;
if (fabs(diffmantissa) > maxmandiff) maxmandiff = fabs(diffmantissa);
}
if (fabs(diffmantissa) > agreediffmantissa) ++disagree; else ++agree;
}
if (val == 0)
os << " ";
else if (pdf == 0)
os << " --- ";
else if (fabs(diffmantissa) > 0.006)
os << format (pdf-val, val, 3, 10);
else
os << " 0 ";
}
os << "\n";
}
os.flags (defaultflags);
if (agree < disagree*agreefraction) {
pass = false;
//os << "agreefraction test failed.\n";
}
if (!pass) ++nfail;
os << "Max abs diff: " << maxabsdiff << ", max rel diff: " << maxreldiff
<< ", max diff mantissa: " << std::fixed << std::setprecision(2) << maxmandiff
<< ", agree/disagree=" << agree << "/" << disagree
<< ", pass=" << pass << std::endl;
}
os << "\n\nNumber of failed tests: " << nfail << std::endl;
os.flags (defaultflags);
os.width(defaultwidth);
os.precision(defaultprecision);
return nfail;
}
static void moments (ROOT::Math::Vavilov& v, double& integral,
double& mean, double& variance,
double& skewness, double& kurtosis) {
int nsteps = 10000;
double t0 = v.GetLambdaMin();
double t1 = v.GetLambdaMax();
double t = t1 - t0;
double dt = t/nsteps;
double sum = 0;
double sumx = 0;
for (int i = 0; i < nsteps; ++i) {
double x = (i+0.5)*dt + t0;
double pdf = v.Pdf(x);
sum += pdf;
sumx += x*pdf;
}
integral = sum*dt;
mean = sumx/sum;
double sumx2 = 0;
double sumx3 = 0;
double sumx4 = 0;
for (int i = 0; i < nsteps; ++i) {
double x = (i+0.5)*dt + t0;
double pdf = v.Pdf(x);
sumx2 += pow(x-mean, 2)*pdf;
sumx3 += pow(x-mean, 3)*pdf;
sumx4 += pow(x-mean, 4)*pdf;
}
variance = sumx2/sum;
skewness = sumx3/sum*pow (variance, -1.5);
kurtosis = sumx4/sum/(variance*variance)-3;
}
void VavilovTest::PrintPdfTable (ROOT::Math::Vavilov& v, std::ostream& os, int digits) {
std::ios::fmtflags defaultflags = os.flags();
int defaultwidth = os.width();
int defaultprecision = os.precision();
for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
double kappa = GetSBKappa (ikappa);
os << "\n\nkappa = " << kappa << "\n\n";
for (int ilambda = 0; ilambda < GetSBNLambda (ikappa); ++ilambda) {
double lambda = GetSBLambda (ikappa, ilambda);
os << std::setw(5) << std::fixed << std::setprecision(2) << lambda;
double x = lambda;
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
v.SetKappaBeta2 (kappa, beta2);
double pdf = v.Pdf(x);
if (pdf > 0)
os << std::setw(digits+7) << std::scientific << std::setprecision(digits-1) << pdf;
else
os << std::setw(digits+7) << " ";
}
os << "\n";
}
os << "Xmin:";
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
v.SetKappaBeta2 (kappa, beta2);
os << std::setw(digits+7) << std::fixed << std::setprecision(digits-1) << v.GetLambdaMin();
}
os << "\n";
os << "Xmax:";
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
v.SetKappaBeta2 (kappa, beta2);
os << std::setw(digits+7) << std::fixed << std::setprecision(digits-1) << v.GetLambdaMax();
}
os << "\n";
os << "int: ";
double integral[11], calcmean[11], calcvariance[11], calcskewness[11], calckurtosis[11];
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
v.SetKappaBeta2 (kappa, beta2);
moments (v, integral[ibeta], calcmean[ibeta], calcvariance[ibeta], calcskewness[ibeta], calckurtosis[ibeta]);
os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << integral[ibeta];
}
os << "\nmean:";
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calcmean[ibeta];
}
os << "\n ";
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Mean (kappa, beta2);
}
os << "\nvar: ";
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calcvariance[ibeta];
}
os << "\n ";
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Variance (kappa, beta2);
}
os << "\nskew:";
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calcskewness[ibeta];
}
os << "\n ";
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Skewness (kappa, beta2);
}
os << "\nkurt:";
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << calckurtosis[ibeta];
}
os << "\n ";
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
os << std::setw(digits+7) << std::fixed << std::setprecision(digits+1) << v.Kurtosis (kappa, beta2);
}
os << "\n";
}
os.flags (defaultflags);
os.width(defaultwidth);
os.precision(defaultprecision);
}
int VavilovTest::CdfTest (ROOT::Math::Vavilov& v, std::ostream& os) {
double maxabsdiff, maxcdfdiff;
GetCdfTestParams (v, maxabsdiff, maxcdfdiff);
return VavilovTest::CdfTest (v, os, maxabsdiff, maxcdfdiff);
}
int VavilovTest::CdfTest (ROOT::Math::Vavilov& v, std::ostream& os, double maxabsdiff, double maxcdfdiff) {
std::ios::fmtflags defaultflags = os.flags();
int defaultwidth = os.width();
int defaultprecision = os.precision();
int nfail = 0;
os << "\n\nTesting Cdf and Cdf_c\n\n";
for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
double kappa = GetSBKappa (ikappa);
os << "\n\nkappa = " << kappa << "\n\n";
double absdiffmax = 0;
double reldiffmax = 0;
double cdfdiffmax = 0;
bool pass = true;
double cdf_calc[11];
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
cdf_calc[ibeta] = 0;
}
for (int ilambda = 0; ilambda < GetSBNLambda (ikappa); ++ilambda) {
double lambda1 = GetSBLambda (ikappa, ilambda);
os << std::setw(5) << std::fixed << std::setprecision(2) << lambda1;
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
v.SetKappaBeta2 (kappa, beta2);
double lambda0 = (ilambda > 0) ? GetSBLambda (ikappa, ilambda-1) : v.GetLambdaMin();
if (lambda1 > lambda0) {
int n = 100;
double dlambda = (lambda1 - lambda0)/n;
for (int i = 0; i < n; ++i) {
double lambda = lambda0 + (i+0.5)*dlambda;
cdf_calc[ibeta] += v.Pdf(lambda)*dlambda;
}
}
double cdf = v.Cdf(lambda1);
double cdf_c = v.Cdf_c(lambda1);
double val = cdf_calc[ibeta];
if (fabs(cdf-val) > absdiffmax) absdiffmax = fabs(cdf-val);
if (fabs(cdf+cdf_c-1) > cdfdiffmax) cdfdiffmax = cdf+cdf_c-1;
if (val > 0 && fabs(cdf/val-1) > reldiffmax) reldiffmax = fabs(cdf/val-1);
if (val == 0)
os << " ";
else
os << std::scientific << std::setw(10) << std::setprecision(2) << cdf-val;
}
os << "\n";
}
os.flags (defaultflags);
if (absdiffmax > maxabsdiff) pass = false;
if (cdfdiffmax > maxcdfdiff) pass = false;
if (!pass) ++nfail;
os << "Max abs diff: " << absdiffmax << ", max rel diff: " << reldiffmax
<< ", max diff cdf+cdf_c-1: " << cdfdiffmax
<< ", pass=" << pass << std::endl;
}
os << "\n\nNumber of failed tests: " << nfail << std::endl;
os.flags (defaultflags);
os.width(defaultwidth);
os.precision(defaultprecision);
return nfail;
}
int VavilovTest::QuantileTest (ROOT::Math::Vavilov& v, std::ostream& os) {
double maxabsdiff;
GetQuantileTestParams (v, maxabsdiff);
return VavilovTest::QuantileTest (v, os, maxabsdiff);
}
int VavilovTest::QuantileTest (ROOT::Math::Vavilov& v, std::ostream& os, double maxabsdiff) {
double qmin = 0;
std::ios::fmtflags defaultflags = os.flags();
int defaultwidth = os.width();
int defaultprecision = os.precision();
int nfail = 0;
static const double qvalues[45] = {0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009,
0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09,
0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99,
0.991, 0.992, 0.993, 0.994, 0.995, 0.996, 0.997, 0.998, 0.999};
os << "\n\nTesting Quantile\n\n";
for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
double kappa = GetSBKappa (ikappa);
os << "\n\nkappa = " << kappa << "\n\n";
double absdiffmax = 0;
bool pass = true;
for (int iq = 0; iq < 45; ++iq) {
double qval = qvalues[iq];
os << std::setw(5) << std::fixed << std::setprecision(3) << qval;
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
v.SetKappaBeta2 (kappa, beta2);
double lambda = v.Quantile (qval);
double cdfval = v.Cdf(lambda);
if (qval > qmin && (1-qval) > qmin) {
if (fabs(cdfval-qval) > absdiffmax) absdiffmax = fabs(cdfval-qval);
os << " " << std::fixed << std::setw(7) << std::setprecision(4) << cdfval-qval << " ";
}
else {
os << " (" << std::fixed << std::setw(7) << std::setprecision(4) << cdfval-qval << ")";
}
}
os << "\n";
}
os.flags (defaultflags);
if (absdiffmax > maxabsdiff) pass = false;
if (!pass) ++nfail;
os << "Max abs diff: " << absdiffmax
<< ", pass=" << pass << std::endl;
}
os << "\n\nTesting Quantile_c\n\n";
for (int ikappa = 0; ikappa < GetSBNKappa(); ++ikappa) {
double kappa = GetSBKappa (ikappa);
os << "\n\nkappa = " << kappa << "\n\n";
double absdiffmax = 0;
bool pass = true;
for (int iq = 0; iq < 45; ++iq) {
double qval = qvalues[iq];
os << std::setw(5) << std::fixed << std::setprecision(3) << qval;
for (int ibeta = 0; ibeta < GetSBNBeta2(); ++ibeta) {
double beta2 = GetSBBeta2 (ibeta);
v.SetKappaBeta2 (kappa, beta2);
double lambda_c = v.Quantile_c (qval);
double cdf_c_val = v.Cdf_c(lambda_c);
if (qval > qmin && (1-qval) > qmin) {
if (fabs(cdf_c_val-qval) > absdiffmax) absdiffmax = fabs(cdf_c_val-qval);
os << " " << std::fixed << std::setw(7) << std::setprecision(4) << cdf_c_val-qval << " ";
}
else {
os << " (" << std::fixed << std::setw(7) << std::setprecision(4) << cdf_c_val-qval << ")";
}
}
os << "\n";
}
os.flags (defaultflags);
if (absdiffmax > maxabsdiff) pass = false;
if (!pass) ++nfail;
os << "Max abs diff: " << absdiffmax
<< ", pass=" << pass << std::endl;
}
os << "\n\nNumber of failed tests: " << nfail << std::endl;
os.flags (defaultflags);
os.width(defaultwidth);
os.precision(defaultprecision);
return nfail;
}
void VavilovTest::GetPdfTestParams (const Vavilov& v, double& maxabsdiff, double& maxdiffmantissa, double& agreefraction, double& agreediffmantissa) {
if (dynamic_cast <const VavilovFast *>(&v)) {
maxabsdiff = 0.08;
maxdiffmantissa = 0.1;
agreefraction = 1;
agreediffmantissa = 0.9;
}
else {
maxabsdiff = 2E-3;
maxdiffmantissa = 0.03;
agreefraction = 5;
agreediffmantissa = 0.015;
}
}
void VavilovTest::GetCdfTestParams (const Vavilov& v, double& maxabsdiff, double& maxcdfdiff) {
if (dynamic_cast <const VavilovFast *>(&v)) {
maxabsdiff = 0.018;
maxcdfdiff = 1E-16;
}
else {
maxabsdiff = 1E-5;
maxcdfdiff = 5E-15;
}
}
void VavilovTest::GetQuantileTestParams (const Vavilov& v, double& maxabsdiff) {
if (dynamic_cast <const VavilovFast *>(&v)) {
maxabsdiff = 0.03;
}
else {
maxabsdiff = 2E-10;
}
}
} // namespace Math
} // namespace ROOT