File indexing completed on 2024-04-06 12:24:20
0001 #ifndef PhisycsTools_Utilities_RootMinuit_h
0002 #define PhisycsTools_Utilities_RootMinuit_h
0003
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);
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);
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 }
0247
0248 #endif
