Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:33:29

0001 #ifndef PhisycsTools_Utilities_RootMinuit_h
0002 #define PhisycsTools_Utilities_RootMinuit_h
0003 /** \class fit::RootMinuit
0004  *
0005  */
0006 #include "PhysicsTools/Utilities/interface/Parameter.h"
0007 #include "PhysicsTools/Utilities/interface/ParameterMap.h"
0008 #include "PhysicsTools/Utilities/interface/RootMinuitResultPrinter.h"
0009 #include "PhysicsTools/Utilities/interface/RootMinuitFuncEvaluator.h"
0010 #include "FWCore/Utilities/interface/EDMException.h"
0011 #include "TMinuit.h"
0012 #include "Math/SMatrix.h"
0013 
0014 #include <vector>
0015 #include <string>
0016 #include <memory>
0017 
0018 namespace fit {
0019 
0020   template <class Function>
0021   class RootMinuit {
0022   public:
0023     RootMinuit(const Function &f, bool verbose = false) : initialized_(false), minValue_(0), verbose_(verbose) {
0024       f_ = f;
0025     }
0026     void addParameter(const std::string &name, std::shared_ptr<double> val, double err, double min, double max) {
0027       if (initialized_)
0028         throw edm::Exception(edm::errors::Configuration)
0029             << "RootMinuit: can't add parameter " << name << " after minuit initialization\n";
0030       pars_.push_back(val);
0031       parameter_t par;
0032       par.val = *val;
0033       par.err = err;
0034       par.min = min;
0035       par.max = max;
0036       par.fixed = false;
0037       parMap_.push_back(std::make_pair(name, par));
0038       size_t s = parIndices_.size();
0039       parIndices_[name] = s;
0040     }
0041     void addParameter(const funct::Parameter &par, double err, double min, double max) {
0042       return addParameter(par.name(), par, err, min, max);
0043     }
0044     double getParameter(const std::string &name, double &err) {
0045       double val;
0046       init();
0047       minuit_->GetParameter(parameterIndex(name), val, err);
0048       return val;
0049     }
0050     double getParameter(const std::string &name) {
0051       double val, err;
0052       init();
0053       minuit_->GetParameter(parameterIndex(name), val, err);
0054       return val;
0055     }
0056     double getParameterError(const std::string &name, double &val) {
0057       double err;
0058       init();
0059       minuit_->GetParameter(parameterIndex(name), val, err);
0060       return err;
0061     }
0062     double getParameterError(const std::string &name) {
0063       double val, err;
0064       init();
0065       minuit_->GetParameter(parameterIndex(name), val, err);
0066       return err;
0067     }
0068     template <unsigned int N>
0069     void getErrorMatrix(ROOT::Math::SMatrix<double, N, N, ROOT::Math::MatRepSym<double, N> > &err) {
0070       init();
0071       if (N != numberOfParameters())
0072         throw edm::Exception(edm::errors::Configuration)
0073             << "RootMinuit: can't call getErrorMatrix passing an SMatrix of dimension " << N
0074             << " while the number of parameters is " << numberOfParameters() << "\n";
0075       double *e = new double[N * N];
0076       minuit_->mnemat(e, numberOfParameters());
0077       for (size_t i = 0; i < N; ++i) {
0078         for (size_t j = 0; j <= i; ++j) {
0079           err(i, j) = e[i + N * j];
0080         }
0081       }
0082       delete[] e;
0083       setParameters();
0084     }
0085     void fixParameter(const std::string &name) {
0086       size_t i = parameterIndex(name);
0087       parMap_[i].second.fixed = true;
0088       if (initialized_) {
0089         minuit_->FixParameter(i);
0090       }
0091     }
0092     void releaseParameter(const std::string &name) {
0093       size_t i = parameterIndex(name);
0094       parMap_[i].second.fixed = false;
0095       if (initialized_) {
0096         minuit_->Release(i);
0097       }
0098     }
0099     void setParameter(const std::string &name, double val) {
0100       size_t i = parameterIndex(name);
0101       parameter_t &par = parMap_[i].second;
0102       par.val = val;
0103       if (initialized_) {
0104         int ierflg = 0;
0105         minuit_->mnparm(i, name, par.val, par.err, par.min, par.max, ierflg);
0106         if (ierflg != 0)
0107           throw edm::Exception(edm::errors::Configuration)
0108               << "RootMinuit: error in setting parameter " << i << " value = " << par.val << " error = " << par.err
0109               << " range = [" << par.min << ", " << par.max << "]\n";
0110       }
0111     }
0112     void setParameters() {
0113       std::map<std::string, size_t>::const_iterator i = parIndices_.begin(), end = parIndices_.end();
0114       double val, err;
0115       for (; i != end; ++i) {
0116         size_t index = i->second;
0117         minuit_->GetParameter(index, val, err);
0118         *pars_[index] = val;
0119       }
0120     }
0121     int numberOfParameters() {
0122       init();
0123       return minuit_->GetNumPars();
0124     }
0125     int numberOfFreeParameters() {
0126       init();
0127       return minuit_->GetNumFreePars();
0128     }
0129     double minimize() {
0130       init();
0131       double arglist[10];
0132       arglist[0] = 5000;
0133       arglist[1] = 0.1;
0134       int ierflag;
0135       minuit_->mnexcm("MINIMIZE", arglist, 2, ierflag);
0136       if (ierflag != 0)
0137         std::cerr << "ERROR in minimize!!" << std::endl;
0138       if (verbose_)
0139         minuit_->mnmatu(1);  //Prints the covariance matrix
0140       double m = minValue();
0141       if (verbose_)
0142         minuit_->mnprin(3, m);
0143       setParameters();
0144       return m;
0145     }
0146     double migrad() {
0147       init();
0148       double arglist[10];
0149       arglist[0] = 5000;
0150       arglist[1] = 0.1;
0151       int ierflag;
0152       minuit_->mnexcm("MIGRAD", arglist, 2, ierflag);
0153       if (ierflag != 0)
0154         std::cerr << "ERROR in migrad!!" << std::endl;
0155       if (verbose_)
0156         minuit_->mnmatu(1);  //Prints the covariance matrix
0157       double m = minValue();
0158       if (verbose_)
0159         minuit_->mnprin(3, m);
0160       setParameters();
0161       return m;
0162     }
0163     double minValue() {
0164       init();
0165       int ierflag;
0166       double edm, errdef;
0167       int nvpar, nparx;
0168       minuit_->mnstat(minValue_, edm, errdef, nvpar, nparx, ierflag);
0169       return minValue_;
0170     }
0171     void printParameters(std::ostream &cout = std::cout) {
0172       std::map<std::string, size_t>::const_iterator i = parIndices_.begin(), end = parIndices_.end();
0173       for (; i != end; ++i) {
0174         cout << i->first << " = " << *pars_[i->second] << " +/- " << getParameterError(i->first) << std::endl;
0175       }
0176     }
0177     void printFitResults(std::ostream &cout = std::cout) {
0178       RootMinuitResultPrinter<Function>::print(minValue(), numberOfFreeParameters(), f_);
0179       printParameters(cout);
0180     }
0181 
0182   private:
0183     parameterVector_t parMap_;
0184     std::map<std::string, size_t> parIndices_;
0185     bool initialized_;
0186     double minValue_;
0187     std::unique_ptr<TMinuit> minuit_;
0188     std::vector<std::shared_ptr<double> > pars_;
0189     static std::vector<std::shared_ptr<double> > *fPars_;
0190     bool verbose_;
0191     static Function f_;
0192     static void fcn_(int &, double *, double &f, double *par, int) {
0193       size_t size = fPars_->size();
0194       for (size_t i = 0; i < size; ++i)
0195         *((*fPars_)[i]) = par[i];
0196       f = RootMinuitFuncEvaluator<Function>::evaluate(f_);
0197     }
0198     size_t parameterIndex(const std::string &name) const {
0199       typename std::map<std::string, size_t>::const_iterator p = parIndices_.find(name);
0200       if (p == parIndices_.end())
0201         throw edm::Exception(edm::errors::Configuration) << "RootMinuit: can't find parameter " << name << "\n";
0202       return p->second;
0203     }
0204     void init() {
0205       if (initialized_)
0206         return;
0207       minuit_.reset(new TMinuit(parMap_.size()));
0208       double arglist[10];
0209       int ierflg = 0;
0210       if (!verbose_) {
0211         arglist[0] = -1;
0212         minuit_->mnexcm("SET PRINT", arglist, 1, ierflg);
0213         if (ierflg != 0)
0214           throw edm::Exception(edm::errors::Configuration) << "RootMinuit: error in calling SET PRINT\n";
0215       }
0216       arglist[0] = 1;
0217       minuit_->mnexcm("SET ERR", arglist, 1, ierflg);
0218       if (ierflg != 0)
0219         throw edm::Exception(edm::errors::Configuration) << "RootMinuit: error in calling SET ERR\n";
0220 
0221       size_t i = 0;
0222       typename parameterVector_t::const_iterator p = parMap_.begin(), end = parMap_.end();
0223       for (; p != end; ++p, ++i) {
0224         const std::string &name = p->first;
0225         const parameter_t &par = p->second;
0226         minuit_->mnparm(i, name, par.val, par.err, par.min, par.max, ierflg);
0227         if (ierflg != 0)
0228           throw edm::Exception(edm::errors::Configuration)
0229               << "RootMinuit: error in setting parameter " << i << " value = " << par.val << " error = " << par.err
0230               << " range = [" << par.min << ", " << par.max << "]\n";
0231       }
0232       initialized_ = true;
0233       for (i = 0, p = parMap_.begin(); p != end; ++p, ++i)
0234         if (p->second.fixed)
0235           minuit_->FixParameter(i);
0236       fPars_ = &pars_;
0237       minuit_->SetFCN(fcn_);
0238     }
0239   };
0240 
0241   template <class Function>
0242   Function RootMinuit<Function>::f_;
0243 
0244   template <class Function>
0245   std::vector<std::shared_ptr<double> > *RootMinuit<Function>::fPars_ = nullptr;
0246 }  // namespace fit
0247 
0248 #endif