File indexing completed on 2024-12-05 02:48:10
0001 #include "lst.h"
0002
0003 #include <typeinfo>
0004
0005 using LSTEvent = ALPAKA_ACCELERATOR_NAMESPACE::lst::LSTEvent;
0006 using namespace ::lst;
0007
0008
0009 int main(int argc, char **argv) {
0010
0011
0012
0013
0014
0015
0016
0017 ana.track_looper_dir_path = gSystem->Getenv("TRACKLOOPERDIR");
0018 if (ana.track_looper_dir_path.IsNull()) {
0019 RooUtil::error(
0020 "TRACKLOOPERDIR is not set! Did you run $ source setup.sh from TrackLooper/ main repository directory?");
0021 }
0022 RooUtil::print(TString::Format("TRACKLOOPERDIR=%s", ana.track_looper_dir_path.Data()));
0023
0024
0025
0026 std::vector<std::string> allArgs(argv, argv + argc);
0027 ana.full_cmd_line = "";
0028 for (auto &str : allArgs) {
0029 ana.full_cmd_line += TString::Format(" %s", str.c_str());
0030 }
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041 cxxopts::Options options("\n $ lst",
0042 "\n **********************\n * *\n * "
0043 "Looper *\n * *\n **********************\n");
0044
0045
0046 options.add_options()("m,mode", "Run mode (NOT DEFINED)", cxxopts::value<int>()->default_value("5"))(
0047 "i,input",
0048 "Comma separated input file list OR if just a directory is provided it will glob all in the directory BUT must "
0049 "end with '/' for the path",
0050 cxxopts::value<std::string>()->default_value("muonGun"))(
0051 "t,tree",
0052 "Name of the tree in the root file to open and loop over",
0053 cxxopts::value<std::string>()->default_value("trackingNtuple/tree"))(
0054 "o,output", "Output file name", cxxopts::value<std::string>())(
0055 "N,nmatch", "N match for MTV-like matching", cxxopts::value<int>()->default_value("9"))(
0056 "p,ptCut", "Min pT cut In GeV", cxxopts::value<float>()->default_value("0.8"))(
0057 "n,nevents", "N events to loop over", cxxopts::value<int>()->default_value("-1"))(
0058 "x,event_index", "specific event index to process", cxxopts::value<int>()->default_value("-1"))(
0059 "g,pdg_id", "The simhit pdgId match option", cxxopts::value<int>()->default_value("0"))(
0060 "v,verbose",
0061 "Verbose mode (0: no print, 1: only final timing, 2: object multiplitcity",
0062 cxxopts::value<int>()->default_value("0"))(
0063 "w,write_ntuple", "Write Ntuple", cxxopts::value<int>()->default_value("1"))(
0064 "s,streams", "Set number of streams", cxxopts::value<int>()->default_value("1"))(
0065 "d,debug", "Run debug job. i.e. overrides output option to 'debug.root' and 'recreate's the file.")(
0066 "l,lower_level", "write lower level objects ntuple results")("G,gnn_ntuple", "write gnn input variable ntuple")(
0067 "j,nsplit_jobs", "Enable splitting jobs by N blocks (--job_index must be set)", cxxopts::value<int>())(
0068 "I,job_index",
0069 "job_index of split jobs (--nsplit_jobs must be set. index starts from 0. i.e. 0, 1, 2, 3, etc...)",
0070 cxxopts::value<int>())("3,tc_pls_triplets", "Allow triplet pLSs in TC collection")(
0071 "2,no_pls_dupclean", "Disable pLS duplicate cleaning (both steps)")("h,help", "Print help");
0072
0073 auto result = options.parse(argc, argv);
0074
0075
0076
0077
0078
0079
0080
0081 if (result.count("help")) {
0082 std::cout << options.help() << std::endl;
0083 exit(1);
0084 }
0085
0086
0087
0088 ana.input_raw_string = result["input"].as<std::string>();
0089
0090
0091 TString TrackingNtupleDir = gSystem->Getenv("TRACKINGNTUPLEDIR");
0092 if (ana.input_raw_string.EqualTo("muonGun"))
0093 ana.input_file_list_tstring = TString::Format("%s/trackingNtuple_10mu_pt_0p5_2.root", TrackingNtupleDir.Data());
0094 else if (ana.input_raw_string.EqualTo("muonGun_highPt"))
0095 ana.input_file_list_tstring = TString::Format("%s/trackingNtuple_10mu_pt_0p5_50.root", TrackingNtupleDir.Data());
0096 else if (ana.input_raw_string.EqualTo("pionGun"))
0097 ana.input_file_list_tstring =
0098 TString::Format("%s/trackingNtuple_6pion_1k_pt_0p5_50.root", TrackingNtupleDir.Data());
0099 else if (ana.input_raw_string.EqualTo("PU200"))
0100 ana.input_file_list_tstring = TString::Format("%s/trackingNtuple_ttbar_PU200.root", TrackingNtupleDir.Data());
0101 else if (ana.input_raw_string.EqualTo("PU200RelVal"))
0102 ana.input_file_list_tstring = TString::Format(
0103 "%s/RelValTTbar_14TeV_CMSSW_12_5_0_pre3/",
0104 (TrackingNtupleDir.Replace(31, 1, "5").Replace(38, 1, "3"))
0105 .Data());
0106 else if (ana.input_raw_string.EqualTo("cube5"))
0107 ana.input_file_list_tstring =
0108 TString::Format("%s/trackingNtuple_10mu_10k_pt_0p5_2_5cm_cube.root", TrackingNtupleDir.Data());
0109 else if (ana.input_raw_string.EqualTo("cube5_highPt"))
0110 ana.input_file_list_tstring =
0111 TString::Format("%s/trackingNtuple_10mu_10k_pt_0p5_50_5cm_cube.root", TrackingNtupleDir.Data());
0112 else if (ana.input_raw_string.EqualTo("cube50"))
0113 ana.input_file_list_tstring =
0114 TString::Format("%s/trackingNtuple_10mu_10k_pt_0p5_2_50cm_cube.root", TrackingNtupleDir.Data());
0115 else if (ana.input_raw_string.EqualTo("cube50_highPt"))
0116 ana.input_file_list_tstring =
0117 TString::Format("%s/trackingNtuple_10mu_10k_pt_0p5_50_50cm_cube.root", TrackingNtupleDir.Data());
0118 else {
0119 ana.input_file_list_tstring = ana.input_raw_string;
0120 }
0121
0122
0123
0124 ana.input_tree_name = result["tree"].as<std::string>();
0125
0126
0127
0128 if (result.count("debug")) {
0129 ana.output_tfile = new TFile("debug.root", "recreate");
0130 } else {
0131
0132
0133 if (result.count("output")) {
0134 ana.output_tfile = new TFile(result["output"].as<std::string>().c_str(), "create");
0135 if (not ana.output_tfile->IsOpen()) {
0136 std::cout << options.help() << std::endl;
0137 std::cout << "ERROR: output already exists! provide new output name or delete old file. OUTPUTFILE="
0138 << result["output"].as<std::string>() << std::endl;
0139 exit(1);
0140 }
0141 } else {
0142 std::cout
0143 << "Warning: Output file name is not provided! Check your arguments. Output file will be set to 'debug.root'"
0144 << std::endl;
0145 ana.output_tfile = new TFile("debug.root", "recreate");
0146 }
0147 }
0148
0149
0150
0151 ana.ptCut = result["ptCut"].as<float>();
0152
0153
0154
0155 ana.nmatch_threshold = result["nmatch"].as<int>();
0156
0157
0158
0159 ana.n_events = result["nevents"].as<int>();
0160 ana.specific_event_index = result["event_index"].as<int>();
0161
0162
0163
0164 ana.pdg_id = result["pdg_id"].as<int>();
0165
0166
0167
0168 if (result.count("nsplit_jobs")) {
0169 ana.nsplit_jobs = result["nsplit_jobs"].as<int>();
0170 if (ana.nsplit_jobs <= 0) {
0171 std::cout << options.help() << std::endl;
0172 std::cout << "ERROR: option string --nsplit_jobs" << ana.nsplit_jobs << " has zero or negative value!"
0173 << std::endl;
0174 std::cout << "I am not sure what this means..." << std::endl;
0175 exit(1);
0176 }
0177 } else {
0178 ana.nsplit_jobs = -1;
0179 }
0180
0181
0182
0183 if (result.count("job_index")) {
0184 ana.job_index = result["job_index"].as<int>();
0185 if (ana.job_index < 0) {
0186 std::cout << options.help() << std::endl;
0187 std::cout << "ERROR: option string --job_index" << ana.job_index << " has negative value!" << std::endl;
0188 std::cout << "I am not sure what this means..." << std::endl;
0189 exit(1);
0190 }
0191 } else {
0192 ana.job_index = -1;
0193 }
0194
0195
0196 if (result.count("job_index") or result.count("nsplit_jobs")) {
0197
0198 if (not(result.count("job_index") and result.count("nsplit_jobs"))) {
0199 std::cout << options.help() << std::endl;
0200 std::cout << "ERROR: option string --job_index and --nsplit_jobs must be set at the same time!" << std::endl;
0201 exit(1);
0202 }
0203
0204 else {
0205 if (ana.job_index >= ana.nsplit_jobs) {
0206 std::cout << options.help() << std::endl;
0207 std::cout << "ERROR: --job_index >= --nsplit_jobs ! This does not make sense..." << std::endl;
0208 exit(1);
0209 }
0210 }
0211 }
0212
0213
0214
0215 ana.verbose = result["verbose"].as<int>();
0216
0217
0218
0219 ana.streams = result["streams"].as<int>();
0220
0221
0222
0223 ana.mode = result["mode"].as<int>();
0224
0225
0226
0227 ana.do_write_ntuple = result["write_ntuple"].as<int>();
0228
0229
0230
0231
0232
0233
0234 if (result.count("lower_level")) {
0235 ana.do_lower_level = true;
0236 } else {
0237 ana.do_lower_level = false;
0238 }
0239
0240
0241
0242 if (result.count("gnn_ntuple")) {
0243 ana.gnn_ntuple = true;
0244
0245 if (not ana.do_write_ntuple) {
0246 std::cout << options.help() << std::endl;
0247 std::cout << "ERROR: option string --write_ntuple 1 and --gnn_ntuple must be set at the same time!" << std::endl;
0248 exit(1);
0249 }
0250 } else {
0251 ana.gnn_ntuple = false;
0252 }
0253
0254
0255
0256 ana.tc_pls_triplets = result["tc_pls_triplets"].as<bool>();
0257
0258
0259
0260 ana.no_pls_dupclean = result["no_pls_dupclean"].as<bool>();
0261
0262
0263 std::cout << "=========================================================" << std::endl;
0264 std::cout << " Running for Acc = " << alpaka::getAccName<ALPAKA_ACCELERATOR_NAMESPACE::Acc3D>() << std::endl;
0265 std::cout << " Setting of the analysis job based on provided arguments " << std::endl;
0266 std::cout << "---------------------------------------------------------" << std::endl;
0267 std::cout << " ana.input_file_list_tstring: " << ana.input_file_list_tstring << std::endl;
0268 std::cout << " ana.output_tfile: " << ana.output_tfile->GetName() << std::endl;
0269 std::cout << " ana.n_events: " << ana.n_events << std::endl;
0270 std::cout << " ana.nsplit_jobs: " << ana.nsplit_jobs << std::endl;
0271 std::cout << " ana.job_index: " << ana.job_index << std::endl;
0272 std::cout << " ana.specific_event_index: " << ana.specific_event_index << std::endl;
0273 std::cout << " ana.do_write_ntuple: " << ana.do_write_ntuple << std::endl;
0274 std::cout << " ana.mode: " << ana.mode << std::endl;
0275 std::cout << " ana.streams: " << ana.streams << std::endl;
0276 std::cout << " ana.verbose: " << ana.verbose << std::endl;
0277 std::cout << " ana.nmatch_threshold: " << ana.nmatch_threshold << std::endl;
0278 std::cout << " ana.tc_pls_triplets: " << ana.tc_pls_triplets << std::endl;
0279 std::cout << " ana.no_pls_dupclean: " << ana.no_pls_dupclean << std::endl;
0280 std::cout << "=========================================================" << std::endl;
0281
0282
0283 ana.events_tchain = RooUtil::FileUtil::createTChain(ana.input_tree_name, ana.input_file_list_tstring);
0284 ana.looper.init(ana.events_tchain, &trk, ana.n_events);
0285 ana.looper.setSilent();
0286
0287
0288 ana.cutflow.setTFile(ana.output_tfile);
0289
0290
0291 ana.output_ttree = new TTree("tree", "tree");
0292
0293
0294 ana.tx = new RooUtil::TTreeX(ana.output_ttree);
0295
0296
0297 writeMetaData();
0298
0299
0300 run_lst();
0301
0302 return 0;
0303 }
0304
0305
0306 void run_lst() {
0307 ALPAKA_ACCELERATOR_NAMESPACE::Device devAcc = alpaka::getDevByIdx(ALPAKA_ACCELERATOR_NAMESPACE::Platform{}, 0u);
0308 std::vector<ALPAKA_ACCELERATOR_NAMESPACE::Queue> queues;
0309 for (int s = 0; s < ana.streams; s++) {
0310 queues.push_back(ALPAKA_ACCELERATOR_NAMESPACE::Queue(devAcc));
0311 }
0312
0313
0314 TStopwatch full_timer;
0315 full_timer.Start();
0316
0317 std::string ptCutString = (ana.ptCut >= 0.8) ? "0.8" : "0.6";
0318 auto hostESData = lst::loadAndFillESHost(ptCutString);
0319 auto deviceESData =
0320 cms::alpakatools::CopyToDevice<lst::LSTESData<alpaka_common::DevHost>>::copyAsync(queues[0], *hostESData.get());
0321 float timeForMapLoading = full_timer.RealTime() * 1000;
0322
0323 if (ana.do_write_ntuple) {
0324 createOutputBranches();
0325 if (ana.gnn_ntuple) {
0326 createGnnNtupleBranches();
0327 }
0328 }
0329
0330 std::vector<std::vector<float>> out_trkX;
0331 std::vector<std::vector<float>> out_trkY;
0332 std::vector<std::vector<float>> out_trkZ;
0333 std::vector<std::vector<unsigned int>> out_hitId;
0334 std::vector<std::vector<unsigned int>> out_hitIdxs;
0335 std::vector<std::vector<unsigned int>> out_hitIndices_vec0;
0336 std::vector<std::vector<unsigned int>> out_hitIndices_vec1;
0337 std::vector<std::vector<unsigned int>> out_hitIndices_vec2;
0338 std::vector<std::vector<unsigned int>> out_hitIndices_vec3;
0339 std::vector<std::vector<float>> out_deltaPhi_vec;
0340 std::vector<std::vector<float>> out_ptIn_vec;
0341 std::vector<std::vector<float>> out_ptErr_vec;
0342 std::vector<std::vector<float>> out_px_vec;
0343 std::vector<std::vector<float>> out_py_vec;
0344 std::vector<std::vector<float>> out_pz_vec;
0345 std::vector<std::vector<float>> out_eta_vec;
0346 std::vector<std::vector<float>> out_etaErr_vec;
0347 std::vector<std::vector<float>> out_phi_vec;
0348 std::vector<std::vector<int>> out_charge_vec;
0349 std::vector<std::vector<unsigned int>> out_seedIdx_vec;
0350 std::vector<std::vector<int>> out_superbin_vec;
0351 std::vector<std::vector<PixelType>> out_pixelType_vec;
0352 std::vector<std::vector<char>> out_isQuad_vec;
0353 std::vector<int> evt_num;
0354 std::vector<TString> file_name;
0355
0356
0357 full_timer.Reset();
0358 full_timer.Start();
0359 while (ana.looper.nextEvent()) {
0360 if (ana.verbose >= 1)
0361 std::cout << "PreLoading event number = " << ana.looper.getCurrentEventIndex() << std::endl;
0362
0363 if (not goodEvent())
0364 continue;
0365
0366 addInputsToLineSegmentTrackingPreLoad(out_trkX,
0367 out_trkY,
0368 out_trkZ,
0369 out_hitId,
0370 out_hitIdxs,
0371 out_hitIndices_vec0,
0372 out_hitIndices_vec1,
0373 out_hitIndices_vec2,
0374 out_hitIndices_vec3,
0375 out_deltaPhi_vec,
0376 out_ptIn_vec,
0377 out_ptErr_vec,
0378 out_px_vec,
0379 out_py_vec,
0380 out_pz_vec,
0381 out_eta_vec,
0382 out_etaErr_vec,
0383 out_phi_vec,
0384 out_charge_vec,
0385 out_seedIdx_vec,
0386 out_superbin_vec,
0387 out_pixelType_vec,
0388 out_isQuad_vec);
0389 evt_num.push_back(ana.looper.getCurrentEventIndex());
0390 file_name.push_back(ana.looper.getCurrentFileName());
0391 }
0392 float timeForInputLoading = full_timer.RealTime() * 1000;
0393
0394 full_timer.Reset();
0395 full_timer.Start();
0396 std::vector<LSTEvent *> events;
0397 for (int s = 0; s < ana.streams; s++) {
0398 LSTEvent *event = new LSTEvent(ana.verbose >= 2, ana.ptCut, queues[s], &deviceESData);
0399 events.push_back(event);
0400 }
0401 float timeForEventCreation = full_timer.RealTime() * 1000;
0402
0403 std::vector<std::vector<float>> timevec;
0404 full_timer.Reset();
0405 full_timer.Start();
0406 float full_elapsed = 0;
0407 #pragma omp parallel num_threads(ana.streams)
0408 {
0409 std::vector<std::vector<float>> timing_information;
0410 float timing_input_loading;
0411 float timing_MD;
0412 float timing_LS;
0413 float timing_T3;
0414 float timing_T5;
0415 float timing_pLS;
0416 float timing_pT5;
0417 float timing_pT3;
0418 float timing_TC;
0419
0420 #pragma omp for
0421 for (int evt = 0; evt < static_cast<int>(out_trkX.size()); evt++) {
0422 if (ana.verbose >= 1)
0423 std::cout << "Running Event number = " << evt << " " << omp_get_thread_num() << std::endl;
0424
0425 events.at(omp_get_thread_num())->initSync();
0426 timing_input_loading = addInputsToEventPreLoad(events.at(omp_get_thread_num()),
0427 false,
0428 out_trkX.at(evt),
0429 out_trkY.at(evt),
0430 out_trkZ.at(evt),
0431 out_hitId.at(evt),
0432 out_hitIdxs.at(evt),
0433 out_hitIndices_vec0.at(evt),
0434 out_hitIndices_vec1.at(evt),
0435 out_hitIndices_vec2.at(evt),
0436 out_hitIndices_vec3.at(evt),
0437 out_deltaPhi_vec.at(evt),
0438 out_ptIn_vec.at(evt),
0439 out_ptErr_vec.at(evt),
0440 out_px_vec.at(evt),
0441 out_py_vec.at(evt),
0442 out_pz_vec.at(evt),
0443 out_eta_vec.at(evt),
0444 out_etaErr_vec.at(evt),
0445 out_phi_vec.at(evt),
0446 out_charge_vec.at(evt),
0447 out_seedIdx_vec.at(evt),
0448 out_superbin_vec.at(evt),
0449 out_pixelType_vec.at(evt),
0450 out_isQuad_vec.at(evt));
0451
0452 timing_MD = runMiniDoublet(events.at(omp_get_thread_num()), evt);
0453 timing_LS = runSegment(events.at(omp_get_thread_num()));
0454 timing_T3 = runT3(events.at(omp_get_thread_num()));
0455 timing_T5 = runQuintuplet(events.at(omp_get_thread_num()));
0456 timing_pLS = runPixelLineSegment(events.at(omp_get_thread_num()), ana.no_pls_dupclean);
0457 timing_pT5 = runPixelQuintuplet(events.at(omp_get_thread_num()));
0458 timing_pT3 = runpT3(events.at(omp_get_thread_num()));
0459 timing_TC = runTrackCandidate(events.at(omp_get_thread_num()), ana.no_pls_dupclean, ana.tc_pls_triplets);
0460
0461 if (ana.verbose == 4) {
0462 #pragma omp critical
0463 {
0464
0465
0466 }
0467 }
0468
0469 if (ana.verbose == 5) {
0470 #pragma omp critical
0471 {
0472
0473 }
0474 }
0475
0476 if (ana.do_write_ntuple) {
0477 #pragma omp critical
0478 {
0479 unsigned int trkev = evt_num.at(evt);
0480 TString fname = file_name.at(evt);
0481 TFile *f = TFile::Open(fname.Data(), "open");
0482 TTree *t = (TTree *)f->Get(ana.input_tree_name.Data());
0483 trk.Init(t);
0484 trk.GetEntry(trkev);
0485 fillOutputBranches(events.at(omp_get_thread_num()));
0486 f->Close();
0487 }
0488 }
0489
0490
0491 TStopwatch my_timer;
0492 my_timer.Start();
0493 events.at(omp_get_thread_num())->resetEventSync();
0494 float timing_resetEvent = my_timer.RealTime();
0495
0496 timing_information.push_back({timing_input_loading,
0497 timing_MD,
0498 timing_LS,
0499 timing_T3,
0500 timing_T5,
0501 timing_pLS,
0502 timing_pT5,
0503 timing_pT3,
0504 timing_TC,
0505 timing_resetEvent});
0506 }
0507
0508 full_elapsed =
0509 full_timer.RealTime() *
0510 1000.f;
0511 #pragma omp critical
0512 timevec.insert(timevec.end(), timing_information.begin(), timing_information.end());
0513 }
0514
0515 float avg_elapsed = full_elapsed / out_trkX.size();
0516
0517 std::cout << "Time for map loading = " << timeForMapLoading << " ms\n";
0518 std::cout << "Time for input loading = " << timeForInputLoading << " ms\n";
0519 std::cout << "Time for event creation = " << timeForEventCreation << " ms\n";
0520 printTimingInformation(timevec, full_elapsed, avg_elapsed);
0521
0522 if (ana.do_write_ntuple) {
0523
0524 ana.output_tfile->cd();
0525 ana.cutflow.saveOutput();
0526 ana.output_ttree->Write();
0527 }
0528
0529 for (int s = 0; s < ana.streams; s++) {
0530 delete events.at(s);
0531 }
0532
0533 delete ana.output_tfile;
0534 }