Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #ifndef PhysicsTools_Utilities_RootMinuitCommands_h
0002 #define PhysicsTools_Utilities_RootMinuitCommands_h
0003 #include "PhysicsTools/Utilities/interface/RootMinuit.h"
0004 #include "PhysicsTools/Utilities/interface/ParameterMap.h"
0005 #include "FWCore/Utilities/interface/EDMException.h"
0006 #include "PhysicsTools/Utilities/interface/Parameter.h"
0007 #include <map>
0008 #include <string>
0009 #include <cstdlib>
0010 #include <fstream>
0011 #include <sstream>
0012 #include <iostream>
0013 #include <boost/tokenizer.hpp>
0014 #include <boost/algorithm/string/join.hpp>
0015 
0016 const char* kParameter = "par";
0017 const char* kFix = "fix";
0018 const char* kRelease = "release";
0019 const char* kSet = "set";
0020 const char* kMinimize = "minimize";
0021 const char* kMigrad = "migrad";
0022 const char* kPrintAll = "print_all";
0023 
0024 namespace fit {
0025 
0026   struct RootMinuitCommand {
0027     std::string name;
0028     std::vector<std::string> stringArgs;
0029     std::vector<double> doubleArgs;
0030     void print(std::ostream& cout) const {
0031       cout << name;
0032       if (!stringArgs.empty()) {
0033         for (size_t i = 0; i != stringArgs.size(); ++i) {
0034           if (i != 0)
0035             cout << ",";
0036           cout << " \"" << stringArgs[i] << "\"";
0037         }
0038       }
0039       if (!doubleArgs.empty()) {
0040         for (size_t i = 0; i != doubleArgs.size(); ++i) {
0041           if (i != 0)
0042             cout << ",";
0043           cout << " " << doubleArgs[i];
0044         }
0045       }
0046     }
0047   };
0048 
0049   template <class Function>
0050   class RootMinuitCommands {
0051   public:
0052     typedef RootMinuit<Function> minuit;
0053     typedef RootMinuitCommand command;
0054     RootMinuitCommands(bool verbose = true) : verbose_(verbose) {}
0055     RootMinuitCommands(const char* fileName, bool verbose = true) : verbose_(verbose) { init(fileName); }
0056     void init(const char* fileName);
0057     double par(const std::string& name) { return parameter(name).val; }
0058     double err(const std::string& name) { return parameter(name).err; }
0059     double min(const std::string& name) { return parameter(name).min; }
0060     double max(const std::string& name) { return parameter(name).max; }
0061     bool fixed(const std::string& name) { return parameter(name).fixed; }
0062     void add(RootMinuit<Function>& minuit, funct::Parameter& p) const {
0063       const std::string& name = p.name();
0064       const parameter_t& par = parameter(name);
0065       minuit.addParameter(p, par.err, par.min, par.max);
0066       if (par.fixed)
0067         minuit.fixParameter(name);
0068     }
0069     void run(RootMinuit<Function>& minuit) const;
0070 
0071   private:
0072     typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
0073     bool verbose_;
0074     unsigned int lineNumber_;
0075     parameterVector_t pars_;
0076     std::map<std::string, size_t> parIndices_;
0077     std::vector<command> commands_;
0078     double string2double(const std::string& str) const {
0079       const char* begin = str.c_str();
0080       char* end;
0081       double val = strtod(begin, &end);
0082       size_t s = end - begin;
0083       if (s < str.size()) {
0084         throw edm::Exception(edm::errors::Configuration) << "RootMinuitCommands: invalid double value: " << str << "\n";
0085       }
0086       return val;
0087     }
0088     const parameter_t& parameter(const std::string& name) const {
0089       typename std::map<std::string, size_t>::const_iterator p = parIndices_.find(name);
0090       if (p == parIndices_.end())
0091         throw edm::Exception(edm::errors::Configuration) << "RootMinuit: can't find parameter " << name << "\n";
0092       return pars_[p->second].second;
0093     }
0094     std::string errorHeader() const {
0095       std::ostringstream out;
0096       out << "RootMinuitCommands config. error, line " << lineNumber_ << ": ";
0097       return out.str();
0098     }
0099     std::string nextToken(typename tokenizer::iterator& i, const typename tokenizer::iterator& end) const {
0100       ++i;
0101       if (i == end)
0102         throw edm::Exception(edm::errors::Configuration) << errorHeader() << "missing parameter\n";
0103       return *i;
0104     }
0105   };
0106 
0107   template <typename Function>
0108   void RootMinuitCommands<Function>::init(const char* fileName) {
0109     using namespace std;
0110     string cmssw_release_base = std::getenv("CMSSW_RELEASE_BASE");
0111     string cmssw_base = std::getenv("CMSSW_BASE");
0112     vector<string> directories;
0113     directories.reserve(3);
0114     directories.emplace_back(".");
0115     if (!cmssw_release_base.empty()) {
0116       directories.emplace_back(cmssw_release_base + "/src");
0117     }
0118     if (!cmssw_base.empty()) {
0119       directories.emplace_back(cmssw_base + "/src");
0120     }
0121     ifstream file;
0122     for (auto const& d : directories) {
0123       std::ifstream f{d + "/" + fileName};
0124       if (f.good()) {
0125         file = std::move(f);
0126         break;
0127       }
0128     }
0129     if (!file.is_open()) {
0130       throw edm::Exception(edm::errors::Configuration)
0131           << "RootMinuitCommands: can't open file: " << fileName
0132           << " in path: " << boost::algorithm::join(directories, ":") << "\n";
0133     }
0134     if (verbose_)
0135       cout << ">>> configuration file: " << fileName << endl;
0136     string line;
0137     lineNumber_ = 0;
0138     bool commands = false;
0139     while (getline(file, line)) {
0140       ++lineNumber_;
0141       if (line.empty())
0142         continue;
0143       char last = *line.rbegin();
0144       if (!(last >= '0' && last <= 'z'))
0145         line.erase(line.end() - 1);
0146       boost::char_separator<char> sep(" ");
0147       tokenizer tokens(line, sep);
0148       tokenizer::iterator i = tokens.begin(), e = tokens.end();
0149       if (tokens.begin() == tokens.end())
0150         continue;
0151       if (*(i->begin()) != '#') {
0152         if (*i == kParameter) {
0153           if (commands)
0154             throw edm::Exception(edm::errors::Configuration)
0155                 << errorHeader() << "please, declare all parameter before all other minuit commands.\n";
0156           string name = nextToken(i, e);
0157           parameter_t par;
0158           par.val = string2double(nextToken(i, e));
0159           par.err = string2double(nextToken(i, e));
0160           par.min = string2double(nextToken(i, e));
0161           par.max = string2double(nextToken(i, e));
0162           tokenizer::iterator j = i;
0163           ++j;
0164           if (j != e) {
0165             string fixed = nextToken(i, e);
0166             if (fixed == "fixed")
0167               par.fixed = true;
0168             else if (fixed == "free")
0169               par.fixed = false;
0170             else
0171               throw edm::Exception(edm::errors::Configuration)
0172                   << errorHeader() << "fix parameter option unknown: " << *i << "\n"
0173                   << "valid options are: fixed, free.\n";
0174           } else {
0175             par.fixed = false;
0176           }
0177           pars_.push_back(std::make_pair(name, par));
0178           size_t s = parIndices_.size();
0179           parIndices_[name] = s;
0180           if (verbose_)
0181             cout << ">>> " << kParameter << " " << name << " " << par.val << " [" << par.min << ", " << par.max << "],"
0182                  << " err: " << par.err << endl;
0183         } else if (*i == kFix || *i == kRelease) {
0184           commands = true;
0185           command com;
0186           com.name = *i;
0187           string arg = nextToken(i, e);
0188           com.stringArgs.push_back(arg);
0189           commands_.push_back(com);
0190           if (verbose_) {
0191             cout << ">>> ";
0192             com.print(cout);
0193             cout << endl;
0194           }
0195         } else if (*i == kSet) {
0196           commands = true;
0197           command com;
0198           com.name = *i;
0199           string arg = nextToken(i, e);
0200           com.stringArgs.push_back(arg);
0201           com.doubleArgs.push_back(string2double(nextToken(i, e)));
0202           commands_.push_back(com);
0203           if (verbose_) {
0204             cout << ">>> ";
0205             com.print(cout);
0206             cout << endl;
0207           }
0208         } else if (*i == kMinimize || *i == kMigrad || *i == kPrintAll) {
0209           commands = true;
0210           command com;
0211           com.name = *i;
0212           commands_.push_back(com);
0213           if (verbose_) {
0214             cout << ">>> ";
0215             com.print(cout);
0216             cout << endl;
0217           }
0218         } else {
0219           throw edm::Exception(edm::errors::Configuration) << errorHeader() << "unkonwn command:: " << *i << "\n";
0220         }
0221       }
0222     }
0223     if (verbose_)
0224       cout << ">>> end configuration" << endl;
0225   }
0226 
0227   template <typename Function>
0228   void RootMinuitCommands<Function>::run(RootMinuit<Function>& minuit) const {
0229     using namespace std;
0230     typename vector<command>::const_iterator c = commands_.begin(), end = commands_.end();
0231     for (; c != end; ++c) {
0232       if (verbose_) {
0233         cout << ">>> minuit command: ";
0234         c->print(cout);
0235         cout << endl;
0236       }
0237       if (c->name == kMinimize)
0238         minuit.minimize();
0239       else if (c->name == kMigrad)
0240         minuit.migrad();
0241       else if (c->name == kPrintAll)
0242         minuit.printFitResults();
0243       else if (c->name == kFix)
0244         minuit.fixParameter(c->stringArgs[0]);
0245       else if (c->name == kRelease)
0246         minuit.releaseParameter(c->stringArgs[0]);
0247       else if (c->name == kSet)
0248         minuit.setParameter(c->stringArgs[0], c->doubleArgs[0]);
0249     }
0250   }
0251 
0252 }  // namespace fit

0253 
0254 #endif