00001 #ifndef __BCMATH__H 00002 #define __BCMATH__H 00003 00004 /*! 00005 * \namespace BCMath 00006 * \brief Some useful mathematic functions. 00007 * \author Daniel Kollar 00008 * \author Kevin Kröninger 00009 * \author Jing Liu 00010 * \version 1.0 00011 * \date 08.2008 00012 * \detail A namespace which encapsulates some mathematical functions 00013 * necessary for BAT. 00014 */ 00015 00016 /* 00017 * Copyright (C) 2008-2010, Daniel Kollar and Kevin Kroeninger. 00018 * All rights reserved. 00019 * 00020 * For the licensing terms see doc/COPYING. 00021 */ 00022 00023 // --------------------------------------------------------- 00024 00025 #define BCMATH_NFACT_ALIMIT 20 00026 00027 // --------------------------------------------------------- 00028 //#include <cstring> 00029 #include <vector> 00030 00031 class TH1D; 00032 00033 namespace BCMath 00034 { 00035 00036 /** 00037 * Calculate the natural logarithm of a gaussian function with mean 00038 * and sigma. If norm=true (default is false) the result is 00039 * multiplied by the normalization constant, i.e. divided by 00040 * sqrt(2*Pi)*sigma. 00041 */ 00042 double LogGaus(double x, double mean = 0, double sigma = 1, bool norm = false); 00043 00044 /** 00045 * Calculate the natural logarithm of a poisson distribution. 00046 */ 00047 double LogPoisson(double x, double par); 00048 00049 /** 00050 * Calculates Binomial probability using approximations for 00051 * factorial calculations if calculation for number greater than 20 00052 * required using the BCMath::ApproxLogFact function. 00053 */ 00054 double ApproxBinomial(int n, int k, double p); 00055 00056 /** 00057 * Calculates natural logarithm of the Binomial probability using 00058 * approximations for factorial calculations if calculation for 00059 * number greater than 20 required using the BCMath::ApproxLogFact 00060 * function. 00061 */ 00062 double LogApproxBinomial(int n, int k, double p); 00063 00064 /** 00065 * Calculates natural logarithm of the Binomial factor "n over k" 00066 * using approximations for factorial calculations if calculation 00067 * for number greater than 20 required using the 00068 * BCMath::ApproxLogFact function. Even for large numbers the 00069 * calculation is performed precisely, if n-k < 5 00070 */ 00071 double LogBinomFactor(int n, int k); 00072 00073 /** 00074 * Calculates natural logarithm of the n-factorial (n!) using 00075 * Srinivasa Ramanujan approximation 00076 * log(n!) = n*log(n) - n + log(n*(1.+4.*n*(1.+2.*n)))/6. + log(PI)/2. 00077 * if n > 20. If n <= 20 it uses BCMath::LogFact to calculate it 00078 * exactly. 00079 */ 00080 double ApproxLogFact(double x); 00081 00082 /** 00083 * Calculates natural logarithm of the Binomial factor "n over k". 00084 */ 00085 double LogNoverK(int n, int k); 00086 00087 /** 00088 * Calculates natural logarithm of the n-factorial (n!) 00089 */ 00090 double LogFact(int n); 00091 00092 /** 00093 * Returns the "greater or equal" of two numbers 00094 */ 00095 inline int Max(int x, int y) 00096 { return x >= y ? x : y; } 00097 00098 inline double Max(double x, double y) 00099 { return x >= y ? x : y; } 00100 00101 /** 00102 * Returns the "less or equal" of two numbers 00103 */ 00104 inline int Min(int x, int y) 00105 { return x <= y ? x : y; } 00106 00107 inline double Min(double x, double y) 00108 { return x <= y ? x : y; } 00109 00110 /** 00111 * Returns the nearest integer of a double number. 00112 */ 00113 int Nint(double x); 00114 00115 /** 00116 * Returns the rms of an array. 00117 */ 00118 double rms(int n, const double * a); 00119 00120 /** 00121 * Calculates the logarithm of the non-relativistic Breit-Wigner 00122 * distribution. 00123 */ 00124 double LogBreitWignerNonRel(double x, double mean, double Gamma, bool norm = false); 00125 double LogBreitWignerRel(double x, double mean, double Gamma); 00126 00127 /** 00128 * Calculates the logarithm of chi square function: 00129 * chi2(double x; size_t n) 00130 */ 00131 double LogChi2(double x, int n); 00132 00133 /** 00134 * Calculates the logarithm of normalized voigtian function: 00135 * voigtian(double x, double sigma, double gamma) 00136 * 00137 * voigtian is a convolution of the following two functions: 00138 * gaussian(x) = 1/(sqrt(2*pi)*sigma) * exp(x*x/(2*sigma*sigma) 00139 * and 00140 * lorentz(x) = (1/pi)*(gamma/2) / (x*x + (gamma/2)*(gamma/2)) 00141 * 00142 * it is singly peaked at x=0. 00143 * The width of the peak is decided by sigma and gamma, 00144 * so they should be positive. 00145 */ 00146 double LogVoigtian(double x, double sigma, double gamma); 00147 00148 /** 00149 * Get N random numbers distributed according to chi square function 00150 * with K degrees of freedom 00151 */ 00152 void RandomChi2(std::vector<double> &randoms, int K); 00153 00154 00155 /** 00156 * Calculate the empirical cumulative distribution function for 00157 * one dimensional data vector. For consistency, the ECDF 00158 * of value smaller than the minimum observed (underflow bin) is zero, and 00159 * for larger than maximum (overflow bin) it is one. 00160 * 00161 * @param data the observations 00162 * @return histogram with normalized ECDF 00163 */ 00164 TH1D* ECDF(const std::vector<double>& data); 00165 00166 /** 00167 * Find the longest runs of zeros and ones in 00168 * the bit stream 00169 * 00170 * @param bitStream input sequence of boolean values 00171 * @return runs first entry the longest zeros run, second entry the longest ones run 00172 */ 00173 00174 std::vector<int> longestRuns(const std::vector<bool>& bitStream); 00175 00176 /** 00177 * Find the longest success/failure run in set of norm. distributed variables. 00178 * Success = observation >= expectation. 00179 * Runs are weighted by the total chi^2 of all elements in the run 00180 * 00181 * @param yMeasured the observations 00182 * @param yExpected the expected values 00183 * @param sigma the theoretical uncertainties on the expectations 00184 * @return runs first entry the max. weight failure run, 00185 * second entry the max. success run */ 00186 00187 std::vector<double> longestRunsChi2(const std::vector<double>& yMeasured, 00188 const std::vector<double>& yExpected, const std::vector<double>& sigma); 00189 00190 /** 00191 * Find the sampling probability that, given n independent Bernoulli 00192 * trials with success rate = failure rate = 1/2, the longest run of 00193 * consecutive successes is greater than the longest observed run. 00194 * Key idea from 00195 * Burr, E.J. & Cane, G. Longest Run of Consecutive Observations Having a Specified Attribute. Biometrika 48, 461-465 (1961). 00196 * 00197 * 00198 * @param longestObserved actual longest run 00199 * @param nTrials number of independent trials 00200 * @return frequency 00201 */ 00202 00203 double longestRunFrequency(unsigned int longestObserved, unsigned int nTrials); 00204 00205 00206 00207 double SplitGaussian(double* x, double* par); 00208 00209 } 00210 00211 // --------------------------------------------------------- 00212 00213 #endif 00214
1.7.1