diff --git a/examples/aod/createO2tables.C b/examples/aod/createO2tables.C index 5a98014..3546be5 100644 --- a/examples/aod/createO2tables.C +++ b/examples/aod/createO2tables.C @@ -14,6 +14,7 @@ R__LOAD_LIBRARY(libDelphesO2) #include "TClonesArray.h" #include "TRandom3.h" #include "TDatabasePDG.h" +#include "TH1F.h" // Delphes includes #include "ExRootAnalysis/ExRootTreeReader.h" @@ -68,8 +69,9 @@ const double forward_rich_sigma = 7.e-3; // Resolution of the Forward RIC const char* inputFileAccMuonPID = "muonAccEffPID.root"; // Simulation parameters -const bool do_vertexing = true; -const bool enable_nuclei = false; +constexpr bool do_vertexing = true; // Vertexing with the O2 +constexpr bool enable_nuclei = true; // Nuclei LUTs +constexpr bool debug_qa = false; // Debug QA histograms int createO2tables(const char* inputFile = "delphes.root", const char* outputFile = "AODRun5.root", @@ -90,11 +92,20 @@ int createO2tables(const char* inputFile = "delphes.root", TDatabasePDG::Instance()->AddParticle("helium3", "helium3", 2.80839160743, kTRUE, 0.0, 6, "Nucleus", 1000020030); TDatabasePDG::Instance()->AddAntiParticle("anti-helium3", -1000020030); - if (do_vertexing) { + if constexpr (do_vertexing) { // Load files for the vertexing o2::base::GeometryManager::loadGeometry("./", false); o2::base::Propagator::initFieldFromGRP("o2sim_grp.root"); } + // Debug histograms + std::map debugHisto; + std::map debugEffNum; + std::map debugEffDen; + std::map debugEffDenPart; + if constexpr (debug_qa) { // Create histograms for debug QA + debugHisto["Multiplicity"] = new TH1F("Multiplicity", "Multiplicity", 1000, 0, 5000); + } + // Create chain of root trees TChain chain("Delphes"); chain.Add(inputFile); @@ -116,7 +127,7 @@ int createO2tables(const char* inputFile = "delphes.root", mapPdgLut.insert(std::make_pair(211, "lutCovm.pi.dat")); mapPdgLut.insert(std::make_pair(321, "lutCovm.ka.dat")); mapPdgLut.insert(std::make_pair(2212, "lutCovm.pr.dat")); - if (enable_nuclei) { + if constexpr (enable_nuclei) { mapPdgLut.insert(std::make_pair(1000010020, "lutCovm.de.dat")); mapPdgLut.insert(std::make_pair(1000010030, "lutCovm.tr.dat")); mapPdgLut.insert(std::make_pair(1000020030, "lutCovm.he3.dat")); @@ -216,8 +227,8 @@ int createO2tables(const char* inputFile = "delphes.root", // Load selected branches with data from specified event treeReader->ReadEntry(ientry); - const float multEtaRange = 2.f; // Range in eta to count the charged particles - float dNdEta = 0.f; // Charged particle multiplicity to use in the efficiency evaluation + constexpr float multEtaRange = 2.f; // Range in eta to count the charged particles + float dNdEta = 0.f; // Charged particle multiplicity to use in the efficiency evaluation for (Int_t iparticle = 0; iparticle < particles->GetEntries(); ++iparticle) { // Loop over particles auto particle = (GenParticle*)particles->At(iparticle); @@ -259,8 +270,17 @@ int createO2tables(const char* inputFile = "delphes.root", dNdEta += 1.f; } FillTree(kMcParticle); + if constexpr (debug_qa) { + if (!debugEffDenPart[particle->PID]) { + debugEffDenPart[particle->PID] = new TH1F(Form("denPart%i", particle->PID), Form("denPart%i;#it{p}_{T} (GeV/#it{c})", particle->PID), 1000, 0, 10); + } + debugEffDenPart[particle->PID]->Fill(particle->PT); + } } dNdEta = 0.5f * dNdEta / multEtaRange; + if constexpr (debug_qa) { + debugHisto["Multiplicity"]->Fill(dNdEta); + } fOffsetLabel += particles->GetEntries(); // For vertexing @@ -300,9 +320,21 @@ int createO2tables(const char* inputFile = "delphes.root", O2Track o2track; // tracks in internal O2 format o2::delphes::TrackUtils::convertTrackToO2Track(*track, o2track, true); + if constexpr (debug_qa) { + if (!debugEffDen[track->PID]) { + debugEffDen[track->PID] = new TH1F(Form("den%i", track->PID), Form("den%i;#it{p}_{T} (GeV/#it{c})", track->PID), 1000, 0, 10); + } + debugEffDen[track->PID]->Fill(track->PT); + } if (!smearer.smearTrack(o2track, track->PID, dNdEta)) { // Skipping inefficient/not correctly smeared tracks continue; } + if constexpr (debug_qa) { + if (!debugEffNum[track->PID]) { + debugEffNum[track->PID] = new TH1F(Form("num%i", track->PID), Form("num%i;#it{p}_{T} (GeV/#it{c})", track->PID), 1000, 0, 10); + } + debugEffNum[track->PID]->Fill(track->PT); + } o2::delphes::TrackUtils::convertO2TrackToTrack(o2track, *track, true); // fill the label tree @@ -421,7 +453,7 @@ int createO2tables(const char* inputFile = "delphes.root", FillTree(kMID); } } - if (do_vertexing) { + if constexpr (do_vertexing) { const float t = (ir.bc2ns() + gRandom->Gaus(0., 100.)) * 1e-3; tracks_for_vertexing.push_back(TrackAlice3{o2track, t, 100.f * 1e-3, TMath::Abs(alabel)}); } @@ -482,7 +514,7 @@ int createO2tables(const char* inputFile = "delphes.root", // fill collision information collision.fIndexBCs = ientry + eventOffset; bc.fGlobalBC = ientry + eventOffset; - if (do_vertexing) { // Performing vertexing + if constexpr (do_vertexing) { // Performing vertexing std::vector lblTracks; std::vector vertices; std::vector vertexTrackIDs; @@ -582,6 +614,19 @@ int createO2tables(const char* inputFile = "delphes.root", if (Trees[i]) Trees[i]->Write(); } + fout->cd(); + for (auto e : debugHisto) { + e.second->Write(); + } + for (auto e : debugEffNum) { + e.second->Write(); + } + for (auto e : debugEffDen) { + e.second->Write(); + } + for (auto e : debugEffDenPart) { + e.second->Write(); + } fout->ls(); fout->Close(); diff --git a/examples/scripts/createO2tables.py b/examples/scripts/createO2tables.py index 8ba57b7..9bde27f 100755 --- a/examples/scripts/createO2tables.py +++ b/examples/scripts/createO2tables.py @@ -42,7 +42,8 @@ def main(configuration_file, turn_off_vertexing, append_production, use_nuclei, - avoid_file_copy): + avoid_file_copy, + debug_aod): arguments = locals() # List of arguments to put into the log parser = configparser.RawConfigParser() parser.read(configuration_file) @@ -121,6 +122,7 @@ def do_copy(in_file, out_file=None, in_path=None): lut_path = opt("lut_path") lut_tag = opt("lut_tag") + lut_tag = f"rmin{int(float(barrel_radius))}.{lut_tag}" lut_particles = ["el", "mu", "pi", "ka", "pr"] if use_nuclei: lut_particles += ["de", "tr", "he3"] @@ -236,10 +238,13 @@ def set_config(config_file, config, value): set_config("createO2tables.C", "const double Bz = ", f"{bField}""e\-1\;/") if turn_off_vertexing: set_config("createO2tables.C", - "const bool do_vertexing = ", "false\;/") + "constexpr bool do_vertexing = ", "false\;/") if use_nuclei: set_config("createO2tables.C", - "const bool enable_nuclei = ", "true\;/") + "constexpr bool enable_nuclei = ", "true\;/") + if debug_aod: + set_config("createO2tables.C", + "constexpr bool debug_qa = ", "true\;/") else: # Check that the geometry file for the vertexing is there if not os.path.isfile("o2sim_grp.root") or not os.path.isfile("o2sim_geometry.root"): run_cmd("mkdir tmpo2sim && cd tmpo2sim && o2-sim -m PIPE ITS MFT -g boxgen -n 1 -j 1 --configKeyValues 'BoxGun.number=1' && cp o2sim_grp.root .. && cp o2sim_geometry.root .. && cd .. && rm -r tmpo2sim") @@ -493,6 +498,9 @@ def write_config(entry, prefix=""): parser.add_argument("--no_nuclei", "--no-nuclei", action="store_true", help="Option use nuclei LUTs") + parser.add_argument("--debug", "-d", + action="store_true", + help="Option to use the debug flag for the AOD making") parser.add_argument("--avoid-config-copy", "--avoid_config_copy", action="store_true", help="Option to avoid copying the configuration files and to use the ones directly in the current path e.g. for grid use") @@ -518,4 +526,5 @@ def write_config(entry, prefix=""): turn_off_vertexing=args.no_vertexing, append_production=args.append, use_nuclei=not args.no_nuclei, - avoid_file_copy=args.avoid_config_copy) + avoid_file_copy=args.avoid_config_copy, + debug_aod=args.debug) diff --git a/examples/scripts/default_configfile.ini b/examples/scripts/default_configfile.ini index de78eb8..256c669 100644 --- a/examples/scripts/default_configfile.ini +++ b/examples/scripts/default_configfile.ini @@ -37,7 +37,7 @@ propagate_card = propagate.2kG.tcl lut_path = $DELPHESO2_ROOT/examples/scripts/ # tag of the LUTs to use in simulation, if the LUT is not created on the fly the minimum track radius is taken from the preexisting LUTs -lut_tag = scenario3 +lut_tag = geometry_v1 # path of the DELPHES aod utilities aod_path = $DELPHESO2_ROOT/examples/aod/