00001 00012 #ifdef _MSC_VER 00013 #include "msdevstudio/MSconfig.h" 00014 #endif 00015 00016 #include "NTupleFCN.h" 00017 00018 #include "datasrcs/DataPointTuple.h" 00019 #include "datasrcs/DataSource.h" 00020 #include "datasrcs/TupleCut.h" 00021 00022 #include "functions/FunctionBase.h" 00023 00024 #include <algorithm> 00025 #include <functional> 00026 00027 using std::bind2nd; 00028 using std::count; 00029 using std::find_if; 00030 using std::not_equal_to; 00031 using std::vector; 00032 00033 using namespace hippodraw; 00034 00035 NTupleFCN:: 00036 NTupleFCN ( ) 00037 : m_fit_cut ( 0 ), 00038 m_ntuple ( 0 ), 00039 m_has_errors ( false ), 00040 m_fit_range ( false ) 00041 { 00042 } 00043 00044 NTupleFCN:: 00045 NTupleFCN ( const NTupleFCN & fcn ) 00046 : StatedFCN ( fcn ), 00047 m_fit_cut ( 0 ), 00048 m_ntuple ( fcn.m_ntuple ), 00049 m_has_errors ( fcn.m_has_errors ), 00050 m_fit_range ( fcn.m_fit_range ) 00051 { 00052 } 00053 00054 void 00055 NTupleFCN:: 00056 copyFrom ( const StatedFCN * base ) 00057 { 00058 StatedFCN::copyFrom ( base ); 00059 00060 const NTupleFCN * fcn = dynamic_cast < const NTupleFCN * > ( base ); 00061 if ( fcn != 0 ) { 00062 m_fit_cut = fcn -> m_fit_cut; 00063 m_fit_range = fcn -> m_fit_range; 00064 m_indices = fcn -> m_indices; 00065 m_ntuple = fcn -> m_ntuple; 00066 m_has_errors = fcn -> m_has_errors; 00067 } 00068 } 00069 00070 namespace dp2 = hippodraw::DataPoint2DTuple; 00071 namespace dp3 = hippodraw::DataPoint3DTuple; 00072 00073 void 00074 NTupleFCN:: 00075 setDataSource ( const DataSource * ntuple ) 00076 { 00077 unsigned int size = ntuple -> columns (); 00078 vector < int > indices ( size ); 00079 00080 for ( unsigned int i = 0; i < size; i++ ) { 00081 indices [i] = i; 00082 } 00083 00084 setDataSource ( ntuple, -1, indices ); 00085 } 00086 00087 void 00088 NTupleFCN:: 00089 setDataSource ( const DataSource * ntuple, 00090 int dimension, 00091 const std::vector < int > & indices ) 00092 { 00093 m_ntuple = ntuple; 00094 m_indices = indices; 00095 } 00096 00097 int 00098 NTupleFCN:: 00099 getErrorColumn () const 00100 { 00101 unsigned int dim = ( m_indices.size() -2 ) / 2; 00102 int ie = m_indices [ 2 * dim + 1 ]; 00103 00104 return ie; 00105 } 00106 00107 bool 00108 NTupleFCN:: 00109 hasErrors ( ) const 00110 { 00111 bool yes = false; 00112 00113 unsigned int ie = getErrorColumn (); 00114 unsigned int cols = m_ntuple -> columns (); 00115 if ( ie < cols ) { 00116 const vector < double > & errors = m_ntuple -> getColumn ( ie ); 00117 00118 if ( errors.empty() ) return false; 00119 00120 vector < double >::const_iterator first 00121 = find_if ( errors.begin(), errors.end (), 00122 bind2nd ( not_equal_to < double > (), 0.0 ) ); 00123 00124 yes = first != errors.end (); 00125 } 00126 00127 return yes; 00128 } 00129 00130 bool 00131 NTupleFCN:: 00132 setUseErrors ( bool yes ) 00133 { 00134 bool didit = false; 00135 if ( yes ) { 00136 if ( hasErrors () ) { 00137 m_has_errors = true; 00138 didit = true; 00139 } 00140 else m_has_errors = false; 00141 } 00142 else { 00143 m_has_errors = false; 00144 didit = true; 00145 } 00146 return didit; 00147 } 00148 00149 bool 00150 NTupleFCN:: 00151 getUseErrors () const 00152 { 00153 return m_has_errors; 00154 } 00155 00156 int 00157 NTupleFCN:: 00158 degreesOfFreedom () const 00159 { 00160 int ie = getErrorColumn (); 00161 00162 const vector < double > & errors = m_ntuple -> getColumn ( ie ); 00163 int number_points = errors.size(); 00164 if ( m_has_errors ) { 00165 int zeros = count ( errors.begin(), errors.end(), 0.0 ); 00166 number_points -= zeros; 00167 } 00168 00169 vector< double > free_parms; 00170 return number_points - getNumberFreeParms (); 00171 } 00172 00173 void 00174 NTupleFCN:: 00175 reset ( std::vector < std::vector < double > > & alpha, 00176 std::vector < double > & beta, 00177 unsigned int size ) 00178 { 00179 beta.clear (); 00180 beta.resize ( size, 0.0 ); 00181 00182 alpha.resize ( size ); 00183 00184 for ( unsigned int i = 0; i < alpha.size (); i++ ) { 00185 alpha[i].clear (); 00186 alpha[i].resize ( size, 0.0 ); 00187 } 00188 } 00189 00192 void 00193 NTupleFCN:: 00194 calcAlphaBeta ( std::vector < std::vector < double > > & alpha, 00195 std::vector < double > & beta ) 00196 { 00197 int ix = m_indices [ dp2::X ]; 00198 int iy = m_indices [ dp2::Y ]; 00199 int ie = m_indices [ dp2::YERR ]; 00200 unsigned int num_parms = getNumberFreeParms (); 00201 reset ( alpha, beta, num_parms ); 00202 00203 unsigned int rows = m_ntuple -> rows (); 00204 for ( unsigned int i = 0; i < rows; i++ ) { 00205 if ( acceptRow ( i ) ) { 00206 const vector < double > & row = m_ntuple -> getRow ( i ); 00207 00208 double err = ie < 0 ? 0. : row [ ie ]; 00209 if ( err == 0.0 && m_has_errors ) continue; 00210 if ( m_has_errors == false ) err = 1.0; 00211 00212 double x = row [ ix ]; 00213 double y = row [ iy ]; 00214 00215 double y_diff = y - m_function -> operator () ( x ); 00216 vector < double > derives; 00217 fillFreeDerivatives ( derives, x ); 00218 00219 for ( unsigned int j = 0; j < num_parms; j++ ) { 00220 double t = derives[j] / ( err * err ); 00221 00222 for ( unsigned int k = 0; k <= j; k++ ) { 00223 alpha[j][k] = alpha[j][k] + t * derives[k]; 00224 } 00225 00226 beta[j] += t * y_diff; 00227 } 00228 } 00229 } 00230 } 00231 00234 void 00235 NTupleFCN:: 00236 setFitCut ( TupleCut * cut ) 00237 { 00238 m_fit_cut = cut; 00239 00240 if ( cut != 0 ) { 00241 int ix = m_indices [ dp2::X ]; 00242 cut -> setColumn ( ix ); 00243 } 00244 } 00245 00246 void 00247 NTupleFCN:: 00248 setFitRange ( bool yes ) 00249 { 00250 m_fit_range = yes; 00251 } 00252 00253 bool 00254 NTupleFCN:: 00255 acceptRow ( unsigned int row ) const 00256 { 00257 bool yes = true; 00258 if ( m_fit_cut != 0 && 00259 m_fit_cut -> isEnabled () && 00260 m_fit_range ) { 00261 yes = m_fit_cut -> acceptRow ( m_ntuple, row ); 00262 } 00263 00264 return yes; 00265 }