efficiency/src/helper.cc

0001 #include "helper.h"
0002
0003 // Global variable
0004 AnalysisConfig ana;
0005
0006 void parseArguments(int argc, char** argv) {
0007   //********************************************************************************
0008   //
0009   // 1. Parsing options
0010   //
0011   //********************************************************************************
0012
0013   // cxxopts is just a tool to parse argc, and argv easily
0014
0015   // Grand option setting
0016   cxxopts::Options options("\n  $ doAnalysis",
0017                            "\n         **********************\n         *                    *\n         *       "
0018                            "Looper       *\n         *                    *\n         **********************\n");
0019
0020   // Read the options
0021   options.add_options()("i,input",
0022                         "Comma separated input file list OR if just a directory is provided it will glob all in the "
0023                         "directory BUT must end with '/' for the path",
0024                         cxxopts::value<std::string>())("o,output", "Output file name", cxxopts::value<std::string>())(
0025       "n,nevents", "N events to loop over", cxxopts::value<int>()->default_value("-1"))(
0026       "j,nsplit_jobs", "Enable splitting jobs by N blocks (--job_index must be set)", cxxopts::value<int>())(
0027       "I,job_index",
0028       "job_index of split jobs (--nsplit_jobs must be set. index starts from 0. i.e. 0, 1, 2, 3, etc...)",
0029       cxxopts::value<int>())(
0030       "g,pdgid", "additional pdgid filtering to use (must be comma separated)", cxxopts::value<std::string>())(
0031       "p,pt_cut", "Transverse momentum cut", cxxopts::value<float>()->default_value("0.9"))(
0032       "e,eta_cut", "Pseudorapidity cut", cxxopts::value<float>()->default_value("4.5"))(
0033       "d,debug", "Run debug job. i.e. overrides output option to 'debug.root' and 'recreate's the file.")("h,help",
0034                                                                                                           "Print help");
0035
0036   auto result = options.parse(argc, argv);
0037
0038   // NOTE: When an option was provided (e.g. -i or --input), then the result.count("<option name>") is more than 0
0039   // Therefore, the option can be parsed easily by asking the condition if (result.count("<option name>");
0040   // That's how the several options are parsed below
0041
0042   //_______________________________________________________________________________
0043   // --help
0044   if (result.count("help")) {
0045     std::cout << options.help() << std::endl;
0046     exit(1);
0047   }
0048
0049   //_______________________________________________________________________________
0050   // --input
0051   if (result.count("input")) {
0052     ana.input_file_list_tstring = result["input"].as<std::string>();
0053   } else {
0054     std::cout << options.help() << std::endl;
0055     std::cout << "ERROR: Input list is not provided! Check your arguments" << std::endl;
0056     exit(1);
0057   }
0058
0059   ana.input_tree_name = "tree";
0060
0061   //_______________________________________________________________________________
0062   // --pdgid
0063   if (result.count("pdgid")) {
0064     std::vector<TString> pdgid_strs = RooUtil::StringUtil::split(result["pdgid"].as<std::string>(), ",");
0065     for (auto& pdgid_str : pdgid_strs) {
0066       ana.pdgids.push_back(pdgid_str.Atoi());
0067     }
0068   }
0069
0070   //_______________________________________________________________________________
0071   // --pt_cut
0072   ana.pt_cut = result["pt_cut"].as<float>();
0073
0074   //_______________________________________________________________________________
0075   // --eta_cut
0076   ana.eta_cut = result["eta_cut"].as<float>();
0077
0078   //_______________________________________________________________________________
0079   // --debug
0080   if (result.count("debug")) {
0081     ana.output_tfile = new TFile("debug.root", "recreate");
0082   } else {
0083     //_______________________________________________________________________________
0084     // --output
0085     if (result.count("output")) {
0086       ana.output_tfile = new TFile(result["output"].as<std::string>().c_str(), "create");
0087       if (not ana.output_tfile->IsOpen()) {
0088         std::cout << options.help() << std::endl;
0089         std::cout << "ERROR: output already exists! provide new output name or delete old file. OUTPUTFILE="
0090                   << result["output"].as<std::string>() << std::endl;
0091         exit(1);
0092       }
0093     } else {
0094       std::cout << options.help() << std::endl;
0095       std::cout << "ERROR: Output file name is not provided! Check your arguments" << std::endl;
0096       exit(1);
0097     }
0098   }
0099
0100   //_______________________________________________________________________________
0101   // --nevents
0102   ana.n_events = result["nevents"].as<int>();
0103
0104   //_______________________________________________________________________________
0105   // --nsplit_jobs
0106   if (result.count("nsplit_jobs")) {
0107     ana.nsplit_jobs = result["nsplit_jobs"].as<int>();
0108     if (ana.nsplit_jobs <= 0) {
0109       std::cout << options.help() << std::endl;
0110       std::cout << "ERROR: option string --nsplit_jobs" << ana.nsplit_jobs << " has zero or negative value!"
0111                 << std::endl;
0112       std::cout << "I am not sure what this means..." << std::endl;
0113       exit(1);
0114     }
0115   } else {
0116     ana.nsplit_jobs = -1;
0117   }
0118
0119   //_______________________________________________________________________________
0120   // --nsplit_jobs
0121   if (result.count("job_index")) {
0122     ana.job_index = result["job_index"].as<int>();
0123     if (ana.job_index < 0) {
0124       std::cout << options.help() << std::endl;
0125       std::cout << "ERROR: option string --job_index" << ana.job_index << " has negative value!" << std::endl;
0126       std::cout << "I am not sure what this means..." << std::endl;
0127       exit(1);
0128     }
0129   } else {
0130     ana.job_index = -1;
0131   }
0132
0133   // Sanity check for split jobs (if one is set the other must be set too)
0134   if (result.count("job_index") or result.count("nsplit_jobs")) {
0135     // If one is not provided then throw error
0136     if (not(result.count("job_index") and result.count("nsplit_jobs"))) {
0137       std::cout << options.help() << std::endl;
0138       std::cout << "ERROR: option string --job_index and --nsplit_jobs must be set at the same time!" << std::endl;
0139       exit(1);
0140     }
0141     // If it is set then check for sanity
0142     else {
0143       if (ana.job_index >= ana.nsplit_jobs) {
0144         std::cout << options.help() << std::endl;
0145         std::cout << "ERROR: --job_index >= --nsplit_jobs ! This does not make sense..." << std::endl;
0146         exit(1);
0147       }
0148     }
0149   }
0150
0151   //
0152   // Printing out the option settings overview
0153   //
0154   std::cout << "=========================================================" << std::endl;
0155   std::cout << " Setting of the analysis job based on provided arguments " << std::endl;
0156   std::cout << "---------------------------------------------------------" << std::endl;
0157   std::cout << " ana.input_file_list_tstring: " << ana.input_file_list_tstring << std::endl;
0158   std::cout << " ana.output_tfile: " << ana.output_tfile->GetName() << std::endl;
0159   std::cout << " ana.n_events: " << ana.n_events << std::endl;
0160   std::cout << " ana.pt_cut: " << ana.pt_cut << std::endl;
0161   std::cout << " ana.eta_cut: " << ana.eta_cut << std::endl;
0162   std::cout << " ana.nsplit_jobs: " << ana.nsplit_jobs << std::endl;
0163   std::cout << " ana.job_index: " << ana.job_index << std::endl;
0164   std::cout << "=========================================================" << std::endl;
0165 }
0166
0167 AnalysisConfig::AnalysisConfig() : tx("variable", "variable") {}
0168
0169 SimTrackSetDefinition::SimTrackSetDefinition(TString set_name_,
0170                                              int pdgid_,
0171                                              int q_,
0172                                              std::function<bool(unsigned int)> pass_,
0173                                              std::function<bool(unsigned int)> sel_) {
0174   set_name = set_name_;
0175   pdgid = pdgid_;
0176   q = q_;
0177   pass = pass_;
0178   sel = sel_;
0179 }
0180
0181 RecoTrackSetDefinition::RecoTrackSetDefinition(TString set_name_,
0182                                                std::function<bool(unsigned int)> pass_,
0183                                                std::function<bool(unsigned int)> sel_,
0184                                                std::function<const std::vector<float>()> pt_,
0185                                                std::function<const std::vector<float>()> eta_,
0186                                                std::function<const std::vector<float>()> phi_,
0187                                                std::function<const std::vector<int>()> type_)
0188     : pt(pt_), eta(eta_), phi(phi_), type(type_) {
0189   set_name = set_name_;
0190   pass = pass_;
0191   sel = sel_;
0192 }
0193
0194 void initializeInputsAndOutputs() {
0195   // Create the TChain that holds the TTree's of the baby ntuples
0196   ana.events_tchain = RooUtil::FileUtil::createTChain(ana.input_tree_name, ana.input_file_list_tstring);
0197
0198   // This is a standard thing SNT does pretty much every looper we use
0199   ana.looper.init(ana.events_tchain, &lstEff, ana.n_events);
0200
0201   // Set the cutflow object output file
0202   ana.cutflow.setTFile(ana.output_tfile);
0203
0204   // Copy the information of the code
0205   std::vector<TString> vstr = RooUtil::StringUtil::split(ana.input_file_list_tstring, ",");
0206   TFile* firstFile = new TFile(vstr[0]);
0207
0208   // Obtain the info from the input file
0209   TString code_tag_data = firstFile->Get("code_tag_data")->GetTitle();
0210   TString gitdiff = firstFile->Get("gitdiff")->GetTitle();
0211   TString input = firstFile->Get("input")->GetTitle();
0212
0213   // cd to output file
0214   ana.output_tfile->cd();
0215
0216   // Write the githash after parsing whether there is any gitdiff
0217   TString githash = RooUtil::StringUtil::split(code_tag_data, "\n")[0];
0218   githash = githash(0, 6);
0219   std::cout << " gitdiff.Length(): " << gitdiff.Length() << std::endl;
0220   if (gitdiff.Length() > 0)
0221     githash += "D";
0222   std::cout << " githash: " << githash << std::endl;
0223   TNamed githash_tnamed("githash", githash.Data());
0224   githash_tnamed.Write();
0225
0226   // Write the sample information
0227   if (input.Contains("PU200"))
0228     input = "PU200";
0229   std::cout << " input: " << input << std::endl;
0230   TNamed input_tnamed("input", input.Data());
0231   input_tnamed.Write();
0232
0233   ana.do_lower_level = false;  // default is false
0234   TObjArray* brobjArray = ana.events_tchain->GetListOfBranches();
0235   for (unsigned int ibr = 0; ibr < (unsigned int)brobjArray->GetEntries(); ++ibr) {
0236     TString brname = brobjArray->At(ibr)->GetName();
0237     if (brname.EqualTo("t5_pt"))
0238       ana.do_lower_level = true;  // if it has the branch it is set to true
0239   }
0240 }
0241
0242 std::vector<float> getPtBounds(int mode) {
0243   std::vector<float> pt_boundaries;
0244   if (mode == 0)
0245     pt_boundaries = {0, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.5, 2.0, 3.0, 5.0, 10, 15., 25, 50};
0246   else if (mode == 1)
0247     pt_boundaries = {0.988, 0.99, 0.992, 0.994, 0.996, 0.998, 1.0, 1.002, 1.004, 1.006, 1.008, 1.01, 1.012};  // lowpt
0248   else if (mode == 2)
0249     pt_boundaries = {0.955, 0.96, 0.965, 0.97, 0.975, 0.98, 0.985, 0.99, 0.995, 1.00,
0250                      1.005, 1.01, 1.015, 1.02, 1.025, 1.03, 1.035, 1.04, 1.045, 1.05};  // pt 0p95 1p05
0251   else if (mode == 3)
0252     pt_boundaries = {
0253         0.5, 0.6, 0.7, 0.8, 0.9, 0.92, 0.94, 0.96, 0.98, 1.0, 1.02, 1.04, 1.06, 1.08, 1.1, 1.2, 1.5};  // lowpt
0254   else if (mode == 4)
0255     pt_boundaries = {0.5,  0.52, 0.54, 0.56, 0.58, 0.6,  0.62, 0.64, 0.66, 0.68, 0.7,  0.72, 0.74,
0256                      0.76, 0.78, 0.8,  0.82, 0.84, 0.86, 0.88, 0.9,  0.92, 0.94, 0.96, 0.98, 1.0,
0257                      1.02, 1.04, 1.06, 1.08, 1.1,  1.12, 1.14, 1.16, 1.18, 1.2,  1.22, 1.24, 1.26,
0258                      1.28, 1.3,  1.32, 1.34, 1.36, 1.38, 1.4,  1.5,  1.6,  1.7,  1.8,  1.9,  2.0};  // lowpt
0259   else if (mode == 5)
0260     pt_boundaries = {0.5,  0.52, 0.54, 0.56, 0.58, 0.6,  0.62, 0.64, 0.66, 0.68, 0.7,  0.72, 0.74,
0261                      0.76, 0.78, 0.8,  0.82, 0.84, 0.86, 0.88, 0.9,  0.92, 0.94, 0.96, 0.98, 1.0,
0262                      1.02, 1.04, 1.06, 1.08, 1.1,  1.12, 1.14, 1.16, 1.18, 1.2,  1.24, 1.28, 1.32,
0263                      1.36, 1.4,  1.6,  1.8,  2.0,  2.2,  2.4,  2.6,  2.8,  3.0};  // lowpt
0264   else if (mode == 6)
0265     pt_boundaries = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 2.0, 3.0, 4.0, 5.0};  // lowpt
0266   else if (mode == 7)
0267     pt_boundaries = {0.5, 0.52, 0.54, 0.56, 0.58, 0.6, 0.62, 0.64, 0.66, 0.68, 0.7, 0.72, 0.74, 0.76, 0.78,
0268                      0.8, 0.82, 0.84, 0.86, 0.88, 0.9, 0.92, 0.94, 0.96, 0.98, 1.0, 1.02, 1.04, 1.06, 1.08,
0269                      1.1, 1.12, 1.14, 1.16, 1.18, 1.2, 1.22, 1.24, 1.26, 1.28, 1.3, 1.32, 1.34, 1.36, 1.38,
0270                      1.4, 1.5,  1.6,  1.7,  1.8,  1.9, 2.0,  2.2,  2.4,  2.6,  2.8, 3,    4,    5,    6,
0271                      7,   8,    9,    10,   15,   20,  25,   30,   35,   40,   45,  50};  // lowpt
0272   else if (mode == 8)
0273     pt_boundaries = {0, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0, 10, 15., 25, 50};
0274   else if (mode == 9) {
0275     for (int i = 0; i < 41; ++i) {
0276       pt_boundaries.push_back(pow(10., -1. + 4. * i / 40.));
0277     }
0278   }
0279   return pt_boundaries;
0280 }