File indexing completed on 2024-04-06 12:19:23
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
0028
0029 CommandLine();
0030 ~CommandLine();
0031
0032
0033
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
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
0061
0062 std::string _exe;
0063 OptionMap_t _options;
0064 StrVec_t _ordered_options;
0065 StrVec_t _unknowns;
0066 };
0067
0068
0069
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
0175
0176
0177 inline CommandLine::CommandLine() {}
0178
0179 inline CommandLine::~CommandLine() {}
0180
0181
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 ( 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
