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     : set_name(set_name_), pass(pass_), sel(sel_), pt(pt_), eta(eta_), phi(phi_), type(type_) {}
0189
0190 void initializeInputsAndOutputs() {
0191   // Create the TChain that holds the TTree's of the baby ntuples
0192   ana.events_tchain = RooUtil::FileUtil::createTChain(ana.input_tree_name, ana.input_file_list_tstring);
0193
0194   // This is a standard thing SNT does pretty much every looper we use
0195   ana.looper.init(ana.events_tchain, &lstEff, ana.n_events);
0196
0197   // Set the cutflow object output file
0198   ana.cutflow.setTFile(ana.output_tfile);
0199
0200   // Copy the information of the code
0201   std::vector<TString> vstr = RooUtil::StringUtil::split(ana.input_file_list_tstring, ",");
0202   TFile* firstFile = new TFile(vstr[0]);
0203
0204   // Obtain the info from the input file
0205   TString code_tag_data = firstFile->Get("code_tag_data")->GetTitle();
0206   TString gitdiff = firstFile->Get("gitdiff")->GetTitle();
0207   TString input = firstFile->Get("input")->GetTitle();
0208
0209   // cd to output file
0210   ana.output_tfile->cd();
0211
0212   // Write the githash after parsing whether there is any gitdiff
0213   TString githash = RooUtil::StringUtil::split(code_tag_data, "\n")[0];
0214   githash = githash(0, 6);
0215   std::cout << " gitdiff.Length(): " << gitdiff.Length() << std::endl;
0216   if (gitdiff.Length() > 0)
0217     githash += "D";
0218   std::cout << " githash: " << githash << std::endl;
0219   TNamed githash_tnamed("githash", githash.Data());
0220   githash_tnamed.Write();
0221
0222   // Write the sample information
0223   if (input.Contains("PU200"))
0224     input = "PU200";
0225   std::cout << " input: " << input << std::endl;
0226   TNamed input_tnamed("input", input.Data());
0227   input_tnamed.Write();
0228
0229   ana.do_lower_level = false;  // default is false
0230   TObjArray* brobjArray = ana.events_tchain->GetListOfBranches();
0231   for (unsigned int ibr = 0; ibr < (unsigned int)brobjArray->GetEntries(); ++ibr) {
0232     TString brname = brobjArray->At(ibr)->GetName();
0233     if (brname.EqualTo("t5_pt"))
0234       ana.do_lower_level = true;  // if it has the branch it is set to true
0235   }
0236 }
0237
0238 std::vector<float> getPtBounds(int mode) {
0239   std::vector<float> pt_boundaries;
0240   if (mode == 0)
0241     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};
0242   else if (mode == 1)
0243     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
0244   else if (mode == 2)
0245     pt_boundaries = {0.955, 0.96, 0.965, 0.97, 0.975, 0.98, 0.985, 0.99, 0.995, 1.00,
0246                      1.005, 1.01, 1.015, 1.02, 1.025, 1.03, 1.035, 1.04, 1.045, 1.05};  // pt 0p95 1p05
0247   else if (mode == 3)
0248     pt_boundaries = {
0249         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
0250   else if (mode == 4)
0251     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,
0252                      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,
0253                      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,
0254                      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
0255   else if (mode == 5)
0256     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,
0257                      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,
0258                      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,
0259                      1.36, 1.4,  1.6,  1.8,  2.0,  2.2,  2.4,  2.6,  2.8,  3.0};  // lowpt
0260   else if (mode == 6)
0261     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
0262   else if (mode == 7)
0263     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,
0264                      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,
0265                      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,
0266                      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,
0267                      7,   8,    9,    10,   15,   20,  25,   30,   35,   40,   45,  50};  // lowpt
0268   else if (mode == 8)
0269     pt_boundaries = {0, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0, 10, 15., 25, 50};
0270   else if (mode == 9) {
0271     for (int i = 0; i < 41; ++i) {
0272       pt_boundaries.push_back(pow(10., -1. + 4. * i / 40.));
0273     }
0274   }
0275   return pt_boundaries;
0276 }