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