00001
00012 #include "MinuitMigrad.h"
00013
00014 #include "StatedFCN.h"
00015
00016 #include "pattern/string_convert.h"
00017
00018 #ifdef HAVE_MINUIT
00019 #include <Minuit/FunctionMinimum.h>
00020 #include <Minuit/MinuitParameter.h>
00021 #include <Minuit/MnMigrad.h>
00022 #include <Minuit/MnUserParameters.h>
00023 #endif
00024
00025 #ifdef HAVE_MINUIT2
00026 #include <Minuit2/FunctionMinimum.h>
00027 #include <Minuit2/MinuitParameter.h>
00028 #include <Minuit2/MnMigrad.h>
00029 #include <Minuit2/MnUserParameters.h>
00030 using namespace ROOT::Minuit2;
00031 #endif
00032
00033 #include <stdexcept>
00034
00035 using std::string;
00036 using std::vector;
00037
00038 using namespace hippodraw;
00039
00040 MinuitMigrad::
00041 MinuitMigrad ( const char * name )
00042 : Fitter ( name ),
00043 m_minimizer ( 0 )
00044 {
00045 }
00046
00047 MinuitMigrad::
00048 MinuitMigrad ( const MinuitMigrad & mm )
00049 : Fitter ( mm ),
00050 m_minimizer ( 0 )
00051 {
00052
00053
00054
00055 }
00056
00057 Fitter *
00058 MinuitMigrad::
00059 clone ( ) const
00060 {
00061 return new MinuitMigrad ( *this );
00062 }
00063 void
00064 MinuitMigrad::
00065 initialize ()
00066 {
00067 m_minimizer = 0;
00068 bool yes = m_fcn -> hasFunction ();
00069 if ( yes ) {
00070 const vector < double > & parms = m_fcn -> getParameters ();
00071 const vector < int > & fixes = m_fcn -> getFixedFlags ();
00072 assert ( fixes.size() == parms.size() );
00073 unsigned int size = parms.size ();
00074 MnUserParameters mn_parms;
00075
00076
00077
00078 for ( unsigned int i = 0; i < size; i++ ) {
00079 const string name = String::convert ( i );
00080 double value = parms [ i ];
00081 #ifdef HAVE_MINUIT
00082 mn_parms.add ( name.c_str(), value, 0.1 );
00083 #else
00084 mn_parms.Add ( name.c_str(), value, 0.1 );
00085 #endif
00086
00087 if ( fixes [ i ] != 0 ) {
00088 #ifdef HAVE_MINUIT
00089 mn_parms.fix ( i );
00090 #else
00091 mn_parms.Fix ( i );
00092 #endif
00093 }
00094 }
00095
00096 m_minimizer = new MnMigrad ( *m_fcn, mn_parms );
00097
00098 if ( m_limits.empty () == false ) {
00099 for ( unsigned int i = 0; i < m_limits.size (); i++ ) {
00100 if ( m_limits[i].set == true ) {
00101 #ifdef HAVE_MINUIT
00102 m_minimizer -> setLimits ( i, m_limits[i].lower, m_limits[i].upper );
00103 #else
00104 m_minimizer -> SetLimits ( i, m_limits[i].lower, m_limits[i].upper );
00105 #endif
00106 }
00107 }
00108 }
00109
00110 }
00111 }
00112
00113 void
00114 MinuitMigrad::
00115 checkIndex ( unsigned int i )
00116 {
00117 if ( m_minimizer == 0 ) initialize ();
00118
00119 if ( m_minimizer == 0 ) {
00120 string what ( m_name );
00121 what += ": model function no yet set";
00122 throw std::runtime_error ( what );
00123 }
00124 const vector < double > & parms = m_fcn -> getParameters ();
00125 unsigned int size = parms.size();
00126 if ( i < size == false ) {
00127 string what ( m_name );
00128 what += ": index to parameter out of range";
00129 throw std::runtime_error ( what );
00130 }
00131 }
00132
00133 void
00134 MinuitMigrad::
00135 initLimits ()
00136 {
00137 const vector < double > & parms = m_fcn -> getParameters ();
00138 for ( unsigned int i = 0; i < parms.size(); i ++ ) {
00139 limit l;
00140 l.set = false;
00141 m_limits.push_back ( l );
00142 }
00143 }
00144
00145 void
00146 MinuitMigrad::
00147 setLimits ( unsigned int i, double lower, double upper )
00148 {
00149 if ( m_limits.empty() == true ) {
00150 initLimits ();
00151 }
00152
00153 checkIndex ( i );
00154
00155 m_limits[i].lower = lower;
00156 m_limits[i].upper = upper;
00157 m_limits[i].set = true;
00158 }
00159
00160 void
00161 MinuitMigrad::
00162 setStepSize ( unsigned int i, double size )
00163 {
00164 checkIndex ( i );
00165 #ifdef HAVE_MINUIT2
00166 m_minimizer -> SetError ( i, size );
00167 #else
00168 m_minimizer -> setError ( i, size );
00169 #endif
00170 }
00171
00178 bool
00179 MinuitMigrad::
00180 calcBestFit ()
00181 {
00182 initialize ();
00183 FunctionMinimum fun_min = m_minimizer -> operator () ();
00184
00185
00186 #ifdef HAVE_MINUIT2
00187 bool yes = fun_min.IsValid ();
00188 #else
00189 bool yes = fun_min.isValid ();
00190 #endif
00191
00192 if ( yes ) {
00193
00194 #ifdef HAVE_MINUIT2
00195 std::vector < double > cur_parms = m_minimizer -> Params();
00196 #else
00197 std::vector < double > cur_parms = m_minimizer -> params();
00198 #endif
00199 m_fcn -> setParameters ( cur_parms );
00200 }
00201 return yes;
00202 }
00203
00204 int
00205 MinuitMigrad::
00206 calcCovariance ( std::vector < std::vector < double > >& covar )
00207 {
00208 if ( m_minimizer == 0 ) initialize ();
00209
00210 #ifdef HAVE_MINUIT2
00211 const MnUserCovariance & covar_m = m_minimizer ->Covariance ();
00212 unsigned int size = covar_m.Nrow ();
00213 #else
00214 const MnUserCovariance & covar_m = m_minimizer ->covariance ();
00215 unsigned int size = covar_m.nrow ();
00216 #endif
00217 covar.resize ( size );
00218 for ( unsigned int i = 0; i < size; i++ ) {
00219 covar[i].resize ( size, 0.0 );
00220 }
00221
00222 for ( unsigned int i = 0; i < size; i++ ) {
00223 for ( unsigned int j = 0; j < size; j++ ) {
00224 covar[i][j] = covar_m ( i, j );
00225 }
00226 }
00227
00228 return EXIT_SUCCESS;
00229 }