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 int Max(unsigned int x, unsigned int y) 00099 { return x >= y ? x : y; } 00100 00101 inline double Max(double x, double y) 00102 { return x >= y ? x : y; } 00103 00104 inline double Max(float x, float y) 00105 { return x >= y ? x : y; } 00106 00107 /** 00108 * Returns the "less or equal" of two numbers 00109 */ 00110 inline int Min(int x, int y) 00111 { return x <= y ? x : y; } 00112 00113 inline int Min(unsigned int x, unsigned int y) 00114 { return x <= y ? x : y; } 00115 00116 inline double Min(double x, double y) 00117 { return x <= y ? x : y; } 00118 00119 inline double Min(float x, float y) 00120 { return x <= y ? x : y; } 00121 00122 /** 00123 * Returns the nearest integer of a double number. 00124 */ 00125 int Nint(double x); 00126 00127 /** 00128 * Returns the rms of an array. 00129 */ 00130 double rms(int n, const double * a); 00131 00132 /** 00133 * Calculates the logarithm of the non-relativistic Breit-Wigner 00134 * distribution. 00135 */ 00136 double LogBreitWignerNonRel(double x, double mean, double Gamma, bool norm = false); 00137 double LogBreitWignerRel(double x, double mean, double Gamma); 00138 00139 /** 00140 * Calculates the logarithm of chi square function: 00141 * chi2(double x; size_t n) 00142 */ 00143 double LogChi2(double x, int n); 00144 00145 /** 00146 * Calculates the logarithm of normalized voigtian function: 00147 * voigtian(double x, double sigma, double gamma) 00148 * 00149 * voigtian is a convolution of the following two functions: 00150 * gaussian(x) = 1/(sqrt(2*pi)*sigma) * exp(x*x/(2*sigma*sigma) 00151 * and 00152 * lorentz(x) = (1/pi)*(gamma/2) / (x*x + (gamma/2)*(gamma/2)) 00153 * 00154 * it is singly peaked at x=0. 00155 * The width of the peak is decided by sigma and gamma, 00156 * so they should be positive. 00157 */ 00158 double LogVoigtian(double x, double sigma, double gamma); 00159 00160 /** 00161 * Get N random numbers distributed according to chi square function 00162 * with K degrees of freedom 00163 */ 00164 void RandomChi2(std::vector<double> &randoms, int K); 00165 00166 00167 /** 00168 * Calculate the empirical cumulative distribution function for 00169 * one dimensional data vector. For consistency, the ECDF 00170 * of value smaller than the minimum observed (underflow bin) is zero, and 00171 * for larger than maximum (overflow bin) it is one. 00172 * 00173 * @param data the observations 00174 * @return histogram with normalized ECDF 00175 */ 00176 TH1D* ECDF(const std::vector<double>& data); 00177 00178 /** 00179 * Find the longest runs of zeros and ones in 00180 * the bit stream 00181 * 00182 * @param bitStream input sequence of boolean values 00183 * @return runs first entry the longest zeros run, second entry the longest ones run 00184 */ 00185 00186 std::vector<int> longestRuns(const std::vector<bool>& bitStream); 00187 00188 /** 00189 * Find the longest success/failure run in set of norm. distributed variables. 00190 * Success = observation >= expectation. 00191 * Runs are weighted by the total chi^2 of all elements in the run 00192 * 00193 * @param yMeasured the observations 00194 * @param yExpected the expected values 00195 * @param sigma the theoretical uncertainties on the expectations 00196 * @return runs first entry the max. weight failure run, 00197 * second entry the max. success run */ 00198 00199 std::vector<double> longestRunsChi2(const std::vector<double>& yMeasured, 00200 const std::vector<double>& yExpected, const std::vector<double>& sigma); 00201 00202 /** 00203 * Find the sampling probability that, given n independent Bernoulli 00204 * trials with success rate = failure rate = 1/2, the longest run of 00205 * consecutive successes is greater than the longest observed run. 00206 * Key idea from 00207 * Burr, E.J. & Cane, G. Longest Run of Consecutive Observations Having a Specified Attribute. Biometrika 48, 461-465 (1961). 00208 * 00209 * 00210 * @param longestObserved actual longest run 00211 * @param nTrials number of independent trials 00212 * @return frequency 00213 */ 00214 00215 double longestRunFrequency(unsigned int longestObserved, unsigned int nTrials); 00216 00217 00218 double SplitGaussian(double* x, double* par); 00219 00220 00221 /** Cache factorials for first \arg \c n integers. 00222 * The cache is filled upon first call of LogFact(). */ 00223 void CacheFactorial(unsigned int n); 00224 00225 } 00226 00227 // --------------------------------------------------------- 00228 00229 #endif 00230
1.7.1