00001 //-------------------------------------------------------------------------- 00002 // File and Version Information: 00003 // $Id: ChisqConsistency.cxx,v 1.1.1.1 2005/04/21 01:17:17 zhangy Exp $ 00004 // 00005 // Description: 00006 // 00007 // Environment: 00008 // Software developed for the BaBar Detector at the SLAC B-Factory. 00009 // 00010 // Author List: 00011 // Bob Jacobsen, Ed Iskander 00012 // 00013 // Copyright Information: 00014 // Copyright (C) 1996 Lawrence Berkeley Laboratory 00015 // 00016 //------------------------------------------------------------------------ 00017 00018 //----------------------- 00019 // This Class's Header -- 00020 //----------------------- 00021 //#include "BaBar/BaBar.hh" 00022 #include "ProbTools/ChisqConsistency.h" 00023 00024 #include <iostream> 00025 #include <math.h> 00026 #include <float.h> 00027 00028 #include "ProbTools/NumRecipes.h" 00029 00030 //#include "ErrLogger/ErrLog.hh" 00031 00032 // prototype the cernlib function 00033 extern "C" { 00034 float chisin_(const float&, const int&); 00035 double dgausn_(double& arg); 00036 } 00037 00038 ChisqConsistency::ChisqConsistency() : 00039 _chisq(-1.0), _nDof(0) 00040 {} 00041 00042 ChisqConsistency::ChisqConsistency(double chisq, double nDof) : 00043 _chisq(chisq), _nDof(nDof) 00044 { 00045 double z2 = 0.5*_chisq; 00046 double n2 = 0.5*_nDof; 00047 00048 if (n2<=0 || z2<0) { 00049 std::cout << "ErrMsg(warning)" << " Got unphysical values: chisq = " << chisq 00050 << " #dof = " << nDof << std::endl; 00051 _value=0; 00052 _likelihood=0; 00053 setStatus(Consistency::unPhysical); 00054 return; 00055 } 00056 setStatus(Consistency::OK); 00057 00058 // given that n2>0 && z2>=0, gammq will NOT abort 00059 _value = NumRecipes::gammq(n2,z2); 00060 00061 if (_chisq==0) { 00062 _likelihood=1; 00063 } else { 00064 double loglike=(n2-1)*log(z2)-z2-NumRecipes::gammln(n2); 00065 if ( loglike < DBL_MIN_EXP ) { 00066 _likelihood = 0; 00067 setStatus(Consistency::underFlow); 00068 } else { 00069 _likelihood = 0.5*exp(loglike); 00070 } 00071 } 00072 } 00073 00074 00075 00076 ChisqConsistency::ChisqConsistency(unsigned nDof, double prob) : 00077 _nDof(nDof) 00078 { 00079 if(prob >= 0.0|| prob <= 1.0 || nDof < 0) 00080 _value = prob; 00081 else { 00082 std::cout << "ErrMsg(warning)" << " Got unphysical values: prob = " << prob 00083 << " #dof = " << nDof << std::endl; 00084 _value=0; 00085 _likelihood=0; 00086 setStatus(Consistency::unPhysical); 00087 return; 00088 } 00089 setStatus(Consistency::OK); 00090 if(prob != 1.0){ 00091 // use the cernlib function to get chisq. Note the funny convention on prob!! 00092 float value = 1.0-float(_value); 00093 int ndof = nDof; 00094 if(value < 1.0) 00095 _chisq = chisin_(value,ndof); 00096 else 00097 _chisq = log(double(FLT_MAX)); 00098 // use the same algorithm as above to get loglikelihood 00099 double z2 = 0.5*_chisq; 00100 double n2 = 0.5*_nDof; 00101 if (_chisq==0) { 00102 _likelihood=1; 00103 } else { 00104 double loglike=(n2-1)*log(z2)-z2-NumRecipes::gammln(n2); 00105 if ( loglike < DBL_MIN_EXP ) { 00106 _likelihood = 0; 00107 setStatus(Consistency::underFlow); 00108 } else { 00109 _likelihood = 0.5*exp(loglike); 00110 } 00111 } 00112 } 00113 } 00114 00115 00116 ChisqConsistency::ChisqConsistency(const ChisqConsistency& other) : 00117 Consistency(other), _chisq(other._chisq), _nDof(other._nDof) 00118 {} 00119 00120 ChisqConsistency& 00121 ChisqConsistency::operator =(const ChisqConsistency& other) { 00122 if(this != &other){ 00123 Consistency::operator =(other); 00124 _chisq = other._chisq; 00125 _nDof = other._nDof; 00126 } 00127 return *this; 00128 }