File indexing completed on 2023-03-17 11:16:55
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 }
0253
0254 #endif
