00001 00012 #ifdef _MSC_VER 00013 #include "msdevstudio/MSconfig.h" 00014 #endif 00015 00016 #include "NTupleLikeliHoodFCN.h" 00017 00018 #include "datasrcs/DataPointTuple.h" 00019 #include "datasrcs/DataSource.h" 00020 #include "functions/FunctionBase.h" 00021 00022 #include <cmath> 00023 #include <cassert> 00024 00025 using std::vector; 00026 using std::log; 00027 using std::pow; 00028 using std::exp; 00029 00030 NTupleLikeliHoodFCN:: 00031 NTupleLikeliHoodFCN () 00032 { 00033 } 00034 00035 NTupleLikeliHoodFCN:: 00036 NTupleLikeliHoodFCN ( const NTupleLikeliHoodFCN & fcn ) 00037 : NTupleFCN ( fcn ) 00038 { 00039 } 00040 00041 StatedFCN * 00042 NTupleLikeliHoodFCN:: 00043 clone () const 00044 { 00045 return new NTupleLikeliHoodFCN ( *this ); 00046 } 00047 00048 namespace dp = hippodraw::DataPoint2DTuple; 00049 00050 double 00051 NTupleLikeliHoodFCN:: 00052 objectiveValue () const 00053 { 00054 double result = 0.0; 00055 double sum1 = 0.0; 00056 double sum2 = 0.0; 00057 00058 int ix = m_indices [ dp::X ]; 00059 int iy = m_indices [ dp::Y ]; 00060 int ie = m_indices [ dp::XERR ]; 00061 00062 int num_bins = m_ntuple -> rows (); 00063 for( int i = 0; i < num_bins; i++ ) 00064 { 00065 if ( acceptRow ( i ) ) { 00066 const vector < double > & row = m_ntuple -> getRow ( i ); 00067 00068 // Get the edges of the bin. 00069 double x = row [ ix ]; 00070 double hwidth = 0.5 * row [ ie ]; 00071 double edgeL = x - hwidth; 00072 double edgeR = x + hwidth; 00073 00074 // Integrate function over the entire i-th bin 00075 double theta = m_function-> integrate( edgeL, edgeR ); 00076 00077 // Get number of enteries in the ith bin 00078 double n = 2.0 * hwidth * row [ iy ]; 00079 00080 // Computing the factorial of n and taking log of it. 00081 // Exact computation is done when N is small ( < 100) else 00082 // we use sterliing's approximations for it. 00083 // CAUTION: Exactly calculating factorial itself usually leads to 00084 // a very large and unrepresentable number so we avoid it. 00085 //double logfactn = 0.0; 00086 //if( n < 100 ) 00087 //{ 00088 // 00089 // for( int j = 2; j <= n; j++) 00090 // logfactn += log( (double) j ); 00091 // 00092 // result += ( n * log ( theta ) - theta - logfactn ); 00093 //} 00094 //else 00095 //{ 00096 // logfactn = n * log( (double) n ) - n; 00097 // result +=( n * ( log ( theta ) - log( (double) n ) + 1 ) - theta ); 00098 //} 00099 00100 //result += ( n * log ( theta ) - theta ); 00101 00102 sum1 += n * log ( theta ); 00103 sum2 += theta; 00104 00105 } 00106 } 00107 result = sum1 - sum2; 00108 // The quantity that is (asymptotically) equivalent to 00109 // chi-square is -2 log( L ) i.e we return -2 * result 00110 00111 return ( -2 * result ); 00112 } 00113 00114 bool 00115 NTupleLikeliHoodFCN:: 00116 needsIntegrated () const 00117 { 00118 return true; 00119 }