Go to the documentation of this file.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 }