CMSSW Reference Manual

00001 #ifndef PhisycsTools_Utilities_RootMinuit_h
00002 #define PhisycsTools_Utilities_RootMinuit_h
00003 
00006 #include "PhysicsTools/Utilities/interface/Parameter.h"
00007 #include "PhysicsTools/Utilities/interface/ParameterMap.h"
00008 #include "FWCore/Utilities/interface/EDMException.h"
00009 #include "TMinuit.h"
00010 #include "TMath.h"
00011 #include "Math/SMatrix.h"
00012 #include <boost/shared_ptr.hpp>
00013 #include <iostream>
00014 #include <vector>
00015 #include <string>
00016 #include <memory>
00017 
00018 namespace fit {
00019   
00020   template<class Function>
00021   class RootMinuit {
00022   public:
00023     RootMinuit(Function f, bool verbose = false) : 
00024       initialized_(false), minValue_(0), verbose_(verbose) { 
00025       f_ = f;
00026     }
00027     void addParameter(const std::string & name, boost::shared_ptr<double> val, double err, double min, double max) {
00028       if(initialized_)
00029         throw edm::Exception(edm::errors::Configuration)
00030           << "RootMinuit: can't add parameter " << name 
00031           << " after minuit initialization\n";      
00032       pars_.push_back(val);
00033       parameter_t par;
00034       par.val = *val;
00035       par.err = err;
00036       par.min = min;
00037       par.max = max;
00038       par.fixed = false;
00039       parMap_.push_back(std::make_pair(name, par));
00040       size_t s = parIndices_.size();
00041       parIndices_[name] = s;
00042     }
00043     void addParameter(const funct::Parameter & par, double err, double min, double max) {
00044       return addParameter(par.name(), par, err, min, max);
00045     }
00046     double getParameter(const std::string & name, double & err) {
00047       double val;
00048       init();
00049       minuit_->GetParameter(parameterIndex(name), val, err);
00050       return val;
00051     }
00052     double getParameter(const std::string & name) {
00053       double val, err;
00054       init();
00055       minuit_->GetParameter(parameterIndex(name), val, err);
00056       return val;
00057     }
00058     double getParameterError(const std::string & name, double & val) {
00059       double err;
00060       init();
00061       minuit_->GetParameter(parameterIndex(name), val, err);
00062       return err;
00063     }
00064     double getParameterError(const std::string & name) {
00065       double val, err;
00066       init();
00067       minuit_->GetParameter(parameterIndex(name), val, err);
00068       return err;
00069     }
00070     template<unsigned int N>
00071     void getErrorMatrix(ROOT::Math::SMatrix<double, N, N, ROOT::Math::MatRepSym<double, N> > & err) {
00072       init();
00073       if(N != numberOfParameters())
00074         throw edm::Exception(edm::errors::Configuration)
00075           << "RootMinuit: can't call getErrorMatrix passing an SMatrix of dimension " << N
00076           << " while the number of parameters is " << numberOfParameters() << "\n";
00077       double * e = new double[N*N];
00078       minuit_->mnemat(e, numberOfParameters());
00079       for(size_t i = 0; i < N; ++i) {
00080         for(size_t j = 0; j <= i; ++j) {
00081           err(i, j) = e[i + N*j];
00082         }
00083       }
00084       delete [] e;
00085       setParameters();
00086     }
00087     void fixParameter(const std::string & name) {
00088       size_t i = parameterIndex(name);
00089       parMap_[i].second.fixed = true;
00090       if(initialized_) {
00091         minuit_->FixParameter(i);
00092       }
00093     }
00094     void releaseParameter(const std::string & name) {
00095       size_t i = parameterIndex(name);
00096       parMap_[i].second.fixed = false;
00097       if(initialized_) {
00098         minuit_->Release(i);
00099       }
00100     }
00101     void setParameter(const std::string & name, double val) {
00102       size_t i = parameterIndex(name);
00103       parameter_t & par = parMap_[i].second;
00104       par.val = val;
00105       if(initialized_) {
00106         int ierflg = 0;      
00107         minuit_->mnparm(i, name, par.val, par.err, par.min, par.max, ierflg);
00108         if(ierflg != 0)
00109           throw edm::Exception(edm::errors::Configuration)
00110             << "RootMinuit: error in setting parameter " << i 
00111             << " value = " << par.val << " error = " << par.err
00112             << " range = [" << par.min << ", " << par.max << "]\n";
00113       }
00114     }
00115     void setParameters() {
00116       std::map<std::string, size_t>::const_iterator i = parIndices_.begin(), end = parIndices_.end();
00117       double val, err;
00118       for(; i != end; ++i) {
00119         size_t index = i->second;
00120         minuit_->GetParameter(index, val, err);
00121         *pars_[index] = val;
00122       }
00123     }
00124     int numberOfParameters() { 
00125       init();
00126       return minuit_->GetNumPars();
00127     }
00128     int numberOfFreeParameters() { 
00129       init();
00130       return minuit_->GetNumFreePars();
00131     }
00132     double minimize() {
00133       init();
00134       double arglist[10];
00135       arglist[0] = 5000;
00136       arglist[1] = 0.1;
00137       int ierflag;
00138       minuit_->mnexcm("MINIMIZE", arglist, 2, ierflag);
00139       if ( ierflag != 0 ) std::cerr << "ERROR in minimize!!" << std::endl;
00140       if(verbose_) minuit_->mnmatu(1); //Prints the covariance matrix
00141       double m = minValue();
00142       if(verbose_) minuit_->mnprin(3, m);
00143       setParameters();
00144       return m;
00145     }
00146     double migrad() {
00147       init();
00148       double arglist[10];
00149       arglist[0] = 5000;
00150       arglist[1] = 0.1;
00151       int ierflag;
00152       minuit_->mnexcm("MIGRAD", arglist, 2, ierflag);
00153       if ( ierflag != 0 ) std::cerr << "ERROR in migrad!!" << std::endl;
00154       if(verbose_) minuit_->mnmatu(1); //Prints the covariance matrix
00155       double m = minValue();
00156       if(verbose_) minuit_->mnprin(3, m);
00157       setParameters();
00158       return m;
00159     }
00160     double minValue() {
00161       init();
00162       int ierflag;
00163       double edm, errdef;
00164       int nvpar, nparx;
00165       minuit_->mnstat(minValue_, edm, errdef, nvpar, nparx, ierflag);
00166       return minValue_;
00167     }
00168     void printParameters(std::ostream& cout = std::cout) {
00169       std::map<std::string, size_t>::const_iterator i = parIndices_.begin(), end = parIndices_.end();
00170       for(; i != end; ++i) {
00171         cout << i->first << " = " << *pars_[i->second] 
00172              << " +/- " << getParameterError(i->first) << std::endl;
00173       }
00174     }
00175     void printFitResults(std::ostream& cout = std::cout) {
00176       double amin = minValue();
00177       int ndof = f_.degreesOfFreedom() - numberOfFreeParameters();
00178         cout << "chi-squared/n.d.o.f. = " << amin << "/" << ndof << " = " << amin/ndof 
00179            << "; prob: " << TMath::Prob(amin, ndof)
00180            << std::endl;
00181       printParameters(cout);
00182     }
00183   private:
00184     parameterVector_t parMap_;
00185     std::map<std::string, size_t> parIndices_;
00186     bool initialized_;
00187     double minValue_;
00188     std::auto_ptr<TMinuit> minuit_;
00189     std::vector<boost::shared_ptr<double> > pars_;
00190     static std::vector<boost::shared_ptr<double> > *fPars_;
00191     bool verbose_;
00192     static Function f_;
00193     static void fcn_(int &, double *, double &f, double *par, int) {
00194       size_t size = fPars_->size();
00195       for(size_t i = 0; i < size; ++i) 
00196         *((*fPars_)[i]) = par[i];
00197       f = f_();
00198     }
00199     size_t parameterIndex(const std::string &name) const {
00200       typename std::map<std::string, size_t>::const_iterator p = parIndices_.find(name);
00201       if(p == parIndices_.end())
00202         throw edm::Exception(edm::errors::Configuration)
00203           << "RootMinuit: can't find parameter " << name << "\n";
00204       return p->second;
00205     }
00206     void init() {
00207       if(initialized_) return;
00208       minuit_.reset(new TMinuit(parMap_.size()));
00209       double arglist[10];
00210       int ierflg = 0;      
00211       if (! verbose_) {
00212         arglist[0] = -1;
00213         minuit_->mnexcm("SET PRINT", arglist, 1, ierflg); 
00214         if (ierflg != 0) 
00215           throw edm::Exception(edm::errors::Configuration)
00216             << "RootMinuit: error in calling SET PRINT\n";
00217       }     
00218       arglist[0] = 1;
00219       minuit_->mnexcm("SET ERR", arglist, 1, ierflg);
00220       if (ierflg != 0) 
00221         throw edm::Exception(edm::errors::Configuration)
00222           << "RootMinuit: error in calling SET ERR\n";
00223 
00224       size_t i = 0;
00225       typename parameterVector_t::const_iterator p = parMap_.begin(), end = parMap_.end();
00226       for(; p != end; ++p, ++i) {
00227         const std::string & name = p->first;
00228         const parameter_t & par = p->second;
00229         minuit_->mnparm(i, name, par.val, par.err, par.min, par.max, ierflg);
00230         if(ierflg != 0)
00231           throw edm::Exception(edm::errors::Configuration)
00232             << "RootMinuit: error in setting parameter " << i 
00233             << " value = " << par.val << " error = " << par.err
00234             << " range = [" << par.min << ", " << par.max << "]\n";
00235       }
00236       initialized_ = true;
00237       for(i = 0, p = parMap_.begin(); p != end; ++p, ++i)
00238         if(p->second.fixed)
00239           minuit_->FixParameter(i);
00240       fPars_= & pars_; 
00241       minuit_->SetFCN(fcn_);
00242     }
00243   };
00244   
00245   template<class Function>
00246   Function RootMinuit<Function>::f_;
00247 
00248   template<class Function>
00249   std::vector<boost::shared_ptr<double> > * RootMinuit<Function>::fPars_ = 0;
00250 }
00251 
00252 #endif
RootMinuit.h
