Back to home page

Project CMSSW displayed by LXR

 
 

    


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   // 0. Preliminary operations
0013   //
0014   //********************************************************************************
0015 
0016   // Checking the TRACKLOOPERDIR is set
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   // Write the command line used to run it
0025   // N.B. This needs to be before the argument parsing as it will change some values
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   // 1. Parsing options
0035   //
0036   //********************************************************************************
0037 
0038   // cxxopts is just a tool to parse argc, and argv easily
0039 
0040   // Grand option setting
0041   cxxopts::Options options("\n  $ lst",
0042                            "\n         **********************\n         *                    *\n         *       "
0043                            "Looper       *\n         *                    *\n         **********************\n");
0044 
0045   // Read the options
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   // NOTE: When an option was provided (e.g. -i or --input), then the result.count("<option name>") is more than 0
0076   // Therefore, the option can be parsed easily by asking the condition if (result.count("<option name>");
0077   // That's how the several options are parsed below
0078 
0079   //_______________________________________________________________________________
0080   // --help
0081   if (result.count("help")) {
0082     std::cout << options.help() << std::endl;
0083     exit(1);
0084   }
0085 
0086   //_______________________________________________________________________________
0087   // --input
0088   ana.input_raw_string = result["input"].as<std::string>();
0089 
0090   // A default value one
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());  // RelVal files under CMSSW_12_5_0_pre3 directory, not CMSSW_12_2_0_pre2 as is the case for the rest of the samples
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   // --tree
0124   ana.input_tree_name = result["tree"].as<std::string>();
0125 
0126   //_______________________________________________________________________________
0127   // --debug
0128   if (result.count("debug")) {
0129     ana.output_tfile = new TFile("debug.root", "recreate");
0130   } else {
0131     //_______________________________________________________________________________
0132     // --output
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   // --ptCut
0151   ana.ptCut = result["ptCut"].as<float>();
0152 
0153   //_______________________________________________________________________________
0154   // --nmatch
0155   ana.nmatch_threshold = result["nmatch"].as<int>();
0156 
0157   //_______________________________________________________________________________
0158   // --nevents
0159   ana.n_events = result["nevents"].as<int>();
0160   ana.specific_event_index = result["event_index"].as<int>();
0161 
0162   //_______________________________________________________________________________
0163   // --pdg_id
0164   ana.pdg_id = result["pdg_id"].as<int>();
0165 
0166   //_______________________________________________________________________________
0167   // --nsplit_jobs
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   // --job_index
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   // Sanity check for split jobs (if one is set the other must be set too)
0196   if (result.count("job_index") or result.count("nsplit_jobs")) {
0197     // If one is not provided then throw error
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     // If it is set then check for sanity
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   // --verbose
0215   ana.verbose = result["verbose"].as<int>();
0216 
0217   //_______________________________________________________________________________
0218   // --streams
0219   ana.streams = result["streams"].as<int>();
0220 
0221   //_______________________________________________________________________________
0222   // --mode
0223   ana.mode = result["mode"].as<int>();
0224 
0225   //_______________________________________________________________________________
0226   // --write_ntuple
0227   ana.do_write_ntuple = result["write_ntuple"].as<int>();
0228 
0229   //_______________________________________________________________________________
0230   // --optimization
0231 
0232   //_______________________________________________________________________________
0233   // --lower_level
0234   if (result.count("lower_level")) {
0235     ana.do_lower_level = true;
0236   } else {
0237     ana.do_lower_level = false;
0238   }
0239 
0240   //_______________________________________________________________________________
0241   // --gnn_ntuple
0242   if (result.count("gnn_ntuple")) {
0243     ana.gnn_ntuple = true;
0244     // If one is not provided then throw error
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   // --tc_pls_triplets
0256   ana.tc_pls_triplets = result["tc_pls_triplets"].as<bool>();
0257 
0258   //_______________________________________________________________________________
0259   // --no_pls_dupclean
0260   ana.no_pls_dupclean = result["no_pls_dupclean"].as<bool>();
0261 
0262   // Printing out the option settings overview
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   // Create the TChain that holds the TTree's of the baby ntuples
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   // Set the cutflow object output file
0288   ana.cutflow.setTFile(ana.output_tfile);
0289 
0290   // Create ttree to output to the ana.output_tfile
0291   ana.output_ttree = new TTree("tree", "tree");
0292 
0293   // Create TTreeX instance that will take care of the interface part of TTree
0294   ana.tx = new RooUtil::TTreeX(ana.output_ttree);
0295 
0296   // Write metadata related to this run
0297   writeMetaData();
0298 
0299   // Run the code
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   // Load various maps used in the lst reconstruction
0314   TStopwatch full_timer;
0315   full_timer.Start();
0316   // Determine which maps to use based on given pt cut for standalone.
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   // Looping input file
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)  // private(event)
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  // nowait// private(event)
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           // TODO BROKEN //
0465           // printAllObjects(events.at(omp_get_thread_num()));
0466         }
0467       }
0468 
0469       if (ana.verbose == 5) {
0470 #pragma omp critical
0471         {
0472           // TODO: debugPrintOutlierMultiplicities
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       // Clear this event
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;  // for loop has implicit barrier I think. So this stops onces all cpu threads have finished but before the next critical section.
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     // Writing ttree output to file
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 }