Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-06-18 02:12:29

0001 #include <string>
0002 #include <vector>
0003 #include <cassert>
0004 #include <iomanip>
0005 #include <iostream>
0006 #include <sstream>
0007 #include <fstream>
0008 #include <map>
0009 #include <cmath>
0010 #include <TFile.h>
0011 #include <TH1F.h>
0012 #include <TF1.h>
0013 #include <TStyle.h>
0014 #include <TMath.h>
0015 
0016 void GetMPV(char name[100], TH1F* histo, TDirectory* Dir, double& peak, double& error, double& sigma, double& err_sigma);
0017 void GetMEAN(TH1F* histo, double& peak, double& error, double& sigma);
0018 void CalculateResponse(bool UseRatioForResponse, double x, double ex, double y, double ey, double& r, double& e);
0019 void CalculateCorrection(bool UseRatioForResponse, double x, double ex, double y, double ey, double& c, double& e);
0020 void Invert(TF1* f, double Min, double Max, double y, double& x);
0021 bool HistoExists(std::vector<std::string> LIST, std::string hname);
0022 int getBin(double x, std::vector<double> boundaries);
0023 
0024 class CommandLine {
0025 public:
0026   //
0027   // construction/destruction
0028   //
0029   CommandLine();
0030   ~CommandLine();
0031 
0032   //
0033   // member functions
0034   //
0035   bool parse(int argc, char** argv);
0036   bool check();
0037   void print();
0038 
0039   template <class T>
0040   T getValue(const std::string& name);
0041   template <class T>
0042   T getValue(const std::string& name, T default_value);
0043 
0044   template <class T>
0045   std::vector<T> getVector(const std::string& name);
0046   template <class T>
0047   std::vector<T> getVector(const std::string& name, const std::string& default_as_string);
0048 
0049 private:
0050   bool parse_file(const std::string& file_name);
0051 
0052 private:
0053   //
0054   // internal typedefs
0055   //
0056   typedef std::map<std::string, std::pair<std::string, bool> > OptionMap_t;
0057   typedef std::vector<std::string> StrVec_t;
0058 
0059   //
0060   // member data
0061   //
0062   std::string _exe;
0063   OptionMap_t _options;
0064   StrVec_t _ordered_options;
0065   StrVec_t _unknowns;
0066 };
0067 
0068 //
0069 // implemenentation of inline functions
0070 //
0071 
0072 //______________________________________________________________________________
0073 template <class T>
0074 T CommandLine::getValue(const std::string& name) {
0075   T result = T();
0076   OptionMap_t::iterator it = _options.find(name);
0077   if (it != _options.end()) {
0078     it->second.second = true;
0079     _ordered_options.push_back(name);
0080     std::stringstream ss;
0081     ss << it->second.first;
0082     ss >> result;
0083     return result;
0084   }
0085   _unknowns.push_back(name);
0086   return result;
0087 }
0088 
0089 //______________________________________________________________________________
0090 template <class T>
0091 T CommandLine::getValue(const std::string& name, T default_value) {
0092   OptionMap_t::const_iterator it = _options.find(name);
0093   if (it != _options.end())
0094     return getValue<T>(name);
0095   std::string default_as_string;
0096   std::stringstream ss;
0097   ss << default_value;
0098   ss >> default_as_string;
0099   _options[name] = std::make_pair(default_as_string, true);
0100   _ordered_options.push_back(name);
0101   return default_value;
0102 }
0103 
0104 //______________________________________________________________________________
0105 template <>
0106 inline bool CommandLine::getValue<bool>(const std::string& name) {
0107   OptionMap_t::iterator it = _options.find(name);
0108   if (it != _options.end()) {
0109     it->second.second = true;
0110     _ordered_options.push_back(name);
0111     std::string val_as_string = it->second.first;
0112     if (val_as_string == "true")
0113       return true;
0114     if (val_as_string == "false")
0115       return false;
0116     int val_as_int;
0117     std::stringstream ss;
0118     ss << val_as_string;
0119     ss >> val_as_int;
0120     return val_as_int;
0121   }
0122   _unknowns.push_back(name);
0123   return false;
0124 }
0125 
0126 //______________________________________________________________________________
0127 template <>
0128 inline bool CommandLine::getValue(const std::string& name, bool default_value) {
0129   OptionMap_t::const_iterator it = _options.find(name);
0130   if (it != _options.end())
0131     return getValue<bool>(name);
0132   _options[name] = (default_value) ? std::make_pair("true", true) : std::make_pair("false", true);
0133   _ordered_options.push_back(name);
0134   return default_value;
0135 }
0136 
0137 //______________________________________________________________________________
0138 template <class T>
0139 std::vector<T> CommandLine::getVector(const std::string& name) {
0140   std::vector<T> result;
0141   OptionMap_t::iterator it = _options.find(name);
0142   if (it != _options.end()) {
0143     it->second.second = true;
0144     _ordered_options.push_back(name);
0145     std::string tmp = it->second.first;
0146     std::string::size_type pos;
0147     if (!tmp.empty()) {
0148       do {
0149         pos = tmp.find(',');
0150         std::stringstream ss;
0151         ss << tmp.substr(0, pos);
0152         tmp.erase(0, pos + 1);
0153         T element;
0154         ss >> element;
0155         result.push_back(element);
0156       } while (pos != std::string::npos);
0157     }
0158   } else {
0159     _unknowns.push_back(name);
0160   }
0161   return result;
0162 }
0163 
0164 //______________________________________________________________________________
0165 template <class T>
0166 std::vector<T> CommandLine::getVector(const std::string& name, const std::string& default_as_string) {
0167   OptionMap_t::iterator it = _options.find(name);
0168   if (it == _options.end())
0169     _options[name] = std::make_pair(default_as_string, false);
0170   return getVector<T>(name);
0171 }
0172 
0173 ////////////////////////////////////////////////////////////////////////////////
0174 // construction / destruction
0175 ////////////////////////////////////////////////////////////////////////////////
0176 //______________________________________________________________________________
0177 inline CommandLine::CommandLine() {}
0178 //______________________________________________________________________________
0179 inline CommandLine::~CommandLine() {}
0180 ////////////////////////////////////////////////////////////////////////////////
0181 // implementation of member functions
0182 ////////////////////////////////////////////////////////////////////////////////
0183 //______________________________________________________________________________
0184 inline bool CommandLine::parse(int argc, char** argv) {
0185   _exe = argv[0];
0186   _options.clear();
0187   _ordered_options.clear();
0188   _unknowns.clear();
0189 
0190   for (int i = 1; i < argc; i++) {
0191     std::string opt = argv[i];
0192     if (0 != opt.find('-')) {
0193       if (i == 1) {
0194         bool success = parse_file(opt);
0195         if (!success)
0196           return false;
0197         continue;
0198       } else {
0199         std::cout << "CommandLine ERROR: options must start with '-'!" << std::endl;
0200         return false;
0201       }
0202     }
0203     opt.erase(0, 1);
0204     std::string next = argv[i + 1];
0205     if (/*0==next.find("-")||*/ i + 1 >= argc) {
0206       std::cout << "ERROR: option '" << opt << "' requires value!" << std::endl;
0207       return false;
0208     }
0209     _options[opt] = std::make_pair(next, false);
0210     i++;
0211     if (i < argc - 1) {
0212       next = argv[i + 1];
0213       while (next.find('-') != 0) {
0214         _options[opt].first += "," + next;
0215         i++;
0216         next = (i < argc - 1) ? argv[i + 1] : "-";
0217       }
0218     }
0219   }
0220 
0221   return true;
0222 }
0223 //______________________________________________________________________________
0224 inline bool CommandLine::check() {
0225   bool result = true;
0226   OptionMap_t::const_iterator it;
0227   for (it = _options.begin(); it != _options.end(); ++it) {
0228     if (!it->second.second) {
0229       std::cout << "CommandLine WARNING: unused option '" << it->first << "'!" << std::endl;
0230       result = false;
0231     }
0232   }
0233 
0234   if (!_unknowns.empty()) {
0235     result = false;
0236     std::cout << "\nCommandLine WARNING: " << _unknowns.size()
0237               << " the followingparameters *must* be provided:" << std::endl;
0238     for (StrVec_t::const_iterator it = _unknowns.begin(); it != _unknowns.end(); ++it)
0239       std::cout << (*it) << std::endl;
0240     std::cout << std::endl;
0241   }
0242   return result;
0243 }
0244 //______________________________________________________________________________
0245 inline void CommandLine::print() {
0246   std::cout << "------------------------------------------------------------" << std::endl;
0247   std::cout << _exe << " options:" << std::endl;
0248   std::cout << "------------------------------------------------------------" << std::endl;
0249   for (StrVec_t::const_iterator itvec = _ordered_options.begin(); itvec != _ordered_options.end(); ++itvec) {
0250     OptionMap_t::const_iterator it = _options.find(*itvec);
0251     assert(it != _options.end());
0252     if (it->second.first.find(",") < std::string::npos) {
0253       std::string tmp = it->second.first;
0254       std::string::size_type length = tmp.length();
0255       std::string::size_type pos;
0256       do {
0257         pos = tmp.find(',');
0258         if (tmp.length() == length) {
0259           std::cout << std::setiosflags(std::ios::left) << std::setw(22) << it->first
0260                     << std::resetiosflags(std::ios::left) << std::setw(3) << "=" << std::setiosflags(std::ios::right)
0261                     << std::setw(35) << tmp.substr(0, pos) << std::resetiosflags(std::ios::right) << std::endl;
0262         } else {
0263           std::cout << std::setiosflags(std::ios::right) << std::setw(60) << tmp.substr(0, pos)
0264                     << std::resetiosflags(std::ios::right) << std::endl;
0265         }
0266         tmp.erase(0, pos + 1);
0267       } while (pos != std::string::npos);
0268     } else {
0269       std::cout << std::setiosflags(std::ios::left) << std::setw(22) << it->first << std::resetiosflags(std::ios::left)
0270                 << std::setw(3) << "=" << std::setiosflags(std::ios::right) << std::setw(35) << it->second.first
0271                 << std::resetiosflags(std::ios::right) << std::endl;
0272     }
0273   }
0274   std::cout << "------------------------------------------------------------" << std::endl;
0275 }
0276 //______________________________________________________________________________
0277 inline bool CommandLine::parse_file(const std::string& file_name) {
0278   std::ifstream fin(file_name.c_str());
0279   if (!fin.is_open()) {
0280     std::cout << "Can't open configuration file " << file_name << std::endl;
0281     return false;
0282   }
0283 
0284   std::stringstream ss;
0285   bool filter(false);
0286   while (!fin.eof()) {
0287     char next;
0288     fin.get(next);
0289     if (!filter && next == '$')
0290       filter = true;
0291     if (!filter) {
0292       if (next == '=')
0293         ss << " " << next << " ";
0294       else
0295         ss << next;
0296     }
0297     if (filter && next == '\n')
0298       filter = false;
0299   }
0300 
0301   std::string token, last_token, key, value;
0302   ss >> token;
0303   while (!ss.eof()) {
0304     if (token == "=") {
0305       if (!key.empty() && !value.empty())
0306         _options[key] = std::make_pair(value, false);
0307       key = last_token;
0308       last_token = "";
0309       value = "";
0310     } else if (!last_token.empty()) {
0311       if (last_token.find('\"') == 0) {
0312         if (last_token.rfind('\"') == last_token.length() - 1) {
0313           last_token = last_token.substr(1, last_token.length() - 2);
0314           value += (!value.empty()) ? "," + last_token : last_token;
0315           last_token = token;
0316         } else
0317           last_token += " " + token;
0318       } else {
0319         value += (!value.empty()) ? "," + last_token : last_token;
0320         last_token = (token == "=") ? "" : token;
0321       }
0322     } else
0323       last_token = (token == "=") ? "" : token;
0324     ss >> token;
0325   }
0326   if (!last_token.empty()) {
0327     if (last_token.find('\"') == 0 && last_token.rfind('\"') == last_token.length() - 1)
0328       last_token = last_token.substr(1, last_token.length() - 2);
0329     value += (!value.empty()) ? "," + last_token : last_token;
0330   }
0331   if (!key.empty() && !value.empty())
0332     _options[key] = std::make_pair(value, false);
0333 
0334   return true;
0335 }
0336 //////////////////////////////////////////////////////////////////////
0337 inline void GetMPV(
0338     char name[100], TH1F* histo, TDirectory* Dir, double& peak, double& error, double& sigma, double& err_sigma) {
0339   double norm, mean, rms, integral, lowlimit, highlimit, LowResponse, HighResponse, a;
0340   int k;
0341   LowResponse = histo->GetXaxis()->GetXmin();
0342   HighResponse = histo->GetXaxis()->GetXmax();
0343   Dir->cd();
0344   TF1* g;
0345   TStyle* myStyle = new TStyle("mystyle", "mystyle");
0346   myStyle->Reset();
0347   myStyle->SetOptFit(1111);
0348   myStyle->SetOptStat(2200);
0349   myStyle->SetStatColor(0);
0350   myStyle->SetTitleFillColor(0);
0351   myStyle->cd();
0352   integral = histo->Integral();
0353   mean = histo->GetMean();
0354   rms = histo->GetRMS();
0355   a = 1.5;
0356   if (integral > 0) {
0357     lowlimit = TMath::Max(LowResponse, mean - a * rms);
0358     highlimit = TMath::Min(mean + a * rms, HighResponse);
0359     norm = histo->GetMaximumStored();
0360     peak = mean;
0361     sigma = rms;
0362     for (k = 0; k < 3; k++) {
0363       g = new TF1("g", "gaus", lowlimit, highlimit);
0364       g->SetParNames("N", "#mu", "#sigma");
0365       g->SetParameter(0, norm);
0366       g->SetParameter(1, peak);
0367       g->SetParameter(2, sigma);
0368       lowlimit = TMath::Max(LowResponse, peak - a * sigma);
0369       highlimit = TMath::Min(peak + a * sigma, HighResponse);
0370       g->SetRange(lowlimit, highlimit);
0371       histo->Fit(g, "RQ");
0372       norm = g->GetParameter(0);
0373       peak = g->GetParameter(1);
0374       sigma = g->GetParameter(2);
0375     }
0376     if (g->GetNDF() > 5) {
0377       peak = g->GetParameter(1);
0378       sigma = g->GetParameter(2);
0379       error = g->GetParError(1);
0380       err_sigma = g->GetParError(2);
0381     } else {
0382       std::cout << "FIT FAILURE: histogram " << name << "...Using MEAN and RMS." << std::endl;
0383       peak = mean;
0384       sigma = rms;
0385       error = histo->GetMeanError();
0386       err_sigma = histo->GetRMSError();
0387     }
0388   } else {
0389     peak = 0;
0390     sigma = 0;
0391     error = 0;
0392     err_sigma = 0;
0393   }
0394   histo->Write();
0395 }
0396 //////////////////////////////////////////////////////////////////////
0397 inline void GetMEAN(TH1F* histo, double& peak, double& error, double& sigma) {
0398   double N = histo->Integral();
0399   if (N > 2) {
0400     peak = histo->GetMean();
0401     sigma = histo->GetRMS();
0402     error = histo->GetMeanError();
0403   } else {
0404     peak = 0;
0405     sigma = 0;
0406     error = 0;
0407   }
0408 }
0409 ///////////////////////////////////////////////////////////////////////
0410 inline void CalculateResponse(bool UseRatioForResponse, double x, double ex, double y, double ey, double& r, double& e) {
0411   if (x > 0 && fabs(y) > 0) {
0412     if (UseRatioForResponse) {
0413       r = y;
0414       e = ey;
0415     } else {
0416       r = (x + y) / x;
0417       e = fabs(r - 1.) * sqrt(pow(ey / y, 2) + pow(ex / x, 2));
0418     }
0419   } else {
0420     r = 0;
0421     e = 0;
0422   }
0423 }
0424 ///////////////////////////////////////////////////////////////////////
0425 inline void CalculateCorrection(
0426     bool UseRatioForResponse, double x, double ex, double y, double ey, double& c, double& e) {
0427   if (x > 0 && fabs(y) > 0) {
0428     if (UseRatioForResponse) {
0429       c = 1. / y;
0430       e = ey / (y * y);
0431     } else {
0432       c = x / (x + y);
0433       e = (fabs(x * y) / pow(x + y, 2)) * sqrt(pow(ey / y, 2) + pow(ex / x, 2));
0434     }
0435   } else {
0436     c = 0;
0437     e = 0;
0438   }
0439 }
0440 ///////////////////////////////////////////////////////////////////////
0441 inline bool HistoExists(std::vector<std::string> LIST, std::string hname) {
0442   unsigned int i, N;
0443   bool found(false);
0444   N = LIST.size();
0445   if (N == 0)
0446     std::cout << "WARNING: empty file histogram list!!!!" << std::endl;
0447   else
0448     for (i = 0; i < N; i++)
0449       if (hname == LIST[i])
0450         found = true;
0451   if (!found)
0452     std::cout << "Histogram: " << hname << " NOT FOUND!!! Check list of existing objects." << std::endl;
0453   return found;
0454 }
0455 ///////////////////////////////////////////////////////////////////////
0456 inline int getBin(double x, std::vector<double> boundaries) {
0457   int i;
0458   int n = boundaries.size() - 1;
0459   if (n <= 0)
0460     return -1;
0461   if (x < boundaries[0] || x >= boundaries[n])
0462     return -1;
0463   for (i = 0; i < n; i++) {
0464     if (x >= boundaries[i] && x < boundaries[i + 1])
0465       return i;
0466   }
0467   return 0;
0468 }
0469 ///////////////////////////////////////////////////////////////////////