diff --git a/Detectors/ITSMFT/ITS/macros/CMakeLists.txt b/Detectors/ITSMFT/ITS/macros/CMakeLists.txt index d651d94fbc07a..75dd430cb7428 100644 --- a/Detectors/ITSMFT/ITS/macros/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/macros/CMakeLists.txt @@ -11,4 +11,5 @@ add_subdirectory(EVE) add_subdirectory(test) +add_subdirectory(tune) add_subdirectory(DCS) diff --git a/Detectors/ITSMFT/ITS/macros/tune/CMakeLists.txt b/Detectors/ITSMFT/ITS/macros/tune/CMakeLists.txt new file mode 100644 index 0000000000000..80fc8bc5e36c2 --- /dev/null +++ b/Detectors/ITSMFT/ITS/macros/tune/CMakeLists.txt @@ -0,0 +1,14 @@ +# Copyright 2019-2026 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +o2_add_test_root_macro(CheckTracklets.C + PUBLIC_LINK_LIBRARIES O2::ITStracking + LABELS its COMPILE_ONLY) diff --git a/Detectors/ITSMFT/ITS/macros/tune/CheckTracklets.C b/Detectors/ITSMFT/ITS/macros/tune/CheckTracklets.C new file mode 100644 index 0000000000000..004c249c45079 --- /dev/null +++ b/Detectors/ITSMFT/ITS/macros/tune/CheckTracklets.C @@ -0,0 +1,714 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// \file CheckTracklets.C +/// \brief Inspect ITS tune-mode tracklet diagnostics. + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +namespace +{ +struct TransitionStats { + long long all{0}; + long long ok{0}; + long long prim{0}; +}; + +struct TransitionPlots { + TH1D* hDPhiAll{nullptr}; + TH1D* hDPhiOk{nullptr}; + TH1D* hDeltaZEventOk{nullptr}; + TH1D* hDeltaZEventPrim{nullptr}; + TH1D* hDPhiTruthFraction{nullptr}; +}; + +struct TransitionIterationPlots { + TH1D* hPhiAll{nullptr}; + TH1D* hPhiOk{nullptr}; + TH1D* hTglAll{nullptr}; + TH1D* hTglOk{nullptr}; + TH1D* hDXYOk{nullptr}; + TH1D* hDXYPrim{nullptr}; + TH1D* hDZOk{nullptr}; + TH1D* hDZPrim{nullptr}; + TH1D* hDRAll{nullptr}; + TH1D* hDROk{nullptr}; + TH1D* hDZGeomAll{nullptr}; + TH1D* hDZGeomOk{nullptr}; + TH1D* hDPhiAll{nullptr}; + TH1D* hDPhiOk{nullptr}; + TH1D* hTglEventOk{nullptr}; + TH1D* hDeltaZEventOk{nullptr}; + TH1D* hDeltaZEventPrim{nullptr}; +}; + +struct IterationData { + int iteration{-1}; + TTree* tree{nullptr}; + std::map, TransitionStats> transitionStats; + TH1D* hPhiAll{nullptr}; + TH1D* hPhiOk{nullptr}; + TH1D* hTglAll{nullptr}; + TH1D* hTglOk{nullptr}; + TH1D* hDXYOk{nullptr}; + TH1D* hDXYPrim{nullptr}; + TH1D* hDZOk{nullptr}; + TH1D* hDZPrim{nullptr}; + TH1D* hDRAll{nullptr}; + TH1D* hDROk{nullptr}; + TH1D* hDZGeomAll{nullptr}; + TH1D* hDZGeomOk{nullptr}; + TH1D* hDPhiAll{nullptr}; + TH1D* hDPhiOk{nullptr}; + TH1D* hTglEventOk{nullptr}; + TH1D* hDeltaZEventOk{nullptr}; + TH1D* hDeltaZEventPrim{nullptr}; + TH2D* hDXYVsPhi{nullptr}; + TH2D* hDZVsTgl{nullptr}; +}; + +struct FormulaSet { + std::unique_ptr from; + std::unique_ptr to; + std::unique_ptr phi; + std::unique_ptr tgl; + std::unique_ptr dr; + std::unique_ptr dz; + std::unique_ptr dPhi; + std::unique_ptr ok; + std::unique_ptr prim; + std::unique_ptr dXY; + std::unique_ptr dZ; + std::unique_ptr tglEvent; + std::unique_ptr deltaZEvent; +}; + +bool parseIteration(const TString& name, int& iteration) +{ + if (!name.BeginsWith("trklt_")) { + return false; + } + const TString suffix = name(6, name.Length() - 6); + if (suffix.IsNull()) { + return false; + } + for (int i = 0; i < suffix.Length(); ++i) { + if (!std::isdigit(static_cast(suffix[i]))) { + return false; + } + } + iteration = suffix.Atoi(); + return true; +} + +std::unique_ptr makeFormula(TTree* tree, const char* name, const char* expression) +{ + auto formula = std::make_unique(name, expression, tree); + if (formula->GetNdim() <= 0) { + std::cerr << "CheckTracklets: missing or unreadable expression '" << expression << "' in tree '" << tree->GetName() << "'" << std::endl; + return nullptr; + } + return formula; +} + +bool makeFormulas(TTree* tree, FormulaSet& formulas) +{ + formulas.from = makeFormula(tree, "fromFormula", "from"); + formulas.to = makeFormula(tree, "toFormula", "to"); + formulas.phi = makeFormula(tree, "phiFormula", "trklt.phi"); + formulas.tgl = makeFormula(tree, "tglFormula", "trklt.tgl"); + formulas.dr = makeFormula(tree, "drFormula", "trklt.dr"); + formulas.dz = makeFormula(tree, "dzFormula", "trklt.dz"); + formulas.dPhi = makeFormula(tree, "dPhiFormula", "trklt.dPhi"); + formulas.ok = makeFormula(tree, "okFormula", "trklt.ok"); + formulas.prim = makeFormula(tree, "primFormula", "trklt.prim"); + formulas.dXY = makeFormula(tree, "dXYFormula", "trklt.dXY"); + formulas.dZ = makeFormula(tree, "dZFormula", "trklt.dZ"); + formulas.tglEvent = makeFormula(tree, "tglEventFormula", "trklt.tglEvent"); + formulas.deltaZEvent = makeFormula(tree, "deltaZEventFormula", "trklt.deltaZEvent"); + return formulas.from && formulas.to && formulas.phi && formulas.tgl && formulas.dr && formulas.dz && formulas.dPhi && formulas.ok && formulas.prim && formulas.dXY && formulas.dZ && formulas.tglEvent && formulas.deltaZEvent; +} + +void updateFormulas(FormulaSet& formulas) +{ + formulas.from->UpdateFormulaLeaves(); + formulas.to->UpdateFormulaLeaves(); + formulas.phi->UpdateFormulaLeaves(); + formulas.tgl->UpdateFormulaLeaves(); + formulas.dr->UpdateFormulaLeaves(); + formulas.dz->UpdateFormulaLeaves(); + formulas.dPhi->UpdateFormulaLeaves(); + formulas.ok->UpdateFormulaLeaves(); + formulas.prim->UpdateFormulaLeaves(); + formulas.dXY->UpdateFormulaLeaves(); + formulas.dZ->UpdateFormulaLeaves(); + formulas.tglEvent->UpdateFormulaLeaves(); + formulas.deltaZEvent->UpdateFormulaLeaves(); +} + +void bookIterationHistograms(IterationData& data) +{ + const TString prefix = Form("iter%d", data.iteration); + data.hPhiAll = new TH1D(Form("%s_phi_all", prefix.Data()), Form("Iteration %d;#phi;tracklets", data.iteration), 160, -TMath::Pi(), TMath::Pi()); + data.hPhiOk = new TH1D(Form("%s_phi_ok", prefix.Data()), Form("Iteration %d truth matched;#phi;tracklets", data.iteration), 160, -TMath::Pi(), TMath::Pi()); + data.hTglAll = new TH1D(Form("%s_tgl_all", prefix.Data()), Form("Iteration %d;tan(#lambda);tracklets", data.iteration), 200, -20., 20.); + data.hTglOk = new TH1D(Form("%s_tgl_ok", prefix.Data()), Form("Iteration %d truth matched;tan(#lambda);tracklets", data.iteration), 200, -20., 20.); + data.hDXYOk = new TH1D(Form("%s_dxy_ok", prefix.Data()), Form("Iteration %d truth matched;d_{XY} to true event (cm);tracklets", data.iteration), 200, 0., 3.); + data.hDXYPrim = new TH1D(Form("%s_dxy_prim", prefix.Data()), Form("Iteration %d primary truth;d_{XY} to true event (cm);tracklets", data.iteration), 200, 0., 3.); + data.hDZOk = new TH1D(Form("%s_dz_ok", prefix.Data()), Form("Iteration %d truth matched;d_{Z} to true event (cm);tracklets", data.iteration), 240, -30., 30.); + data.hDZPrim = new TH1D(Form("%s_dz_prim", prefix.Data()), Form("Iteration %d primary truth;d_{Z} to true event (cm);tracklets", data.iteration), 240, -30., 30.); + data.hDXYVsPhi = new TH2D(Form("%s_dxy_vs_phi", prefix.Data()), Form("Iteration %d truth matched;#phi;d_{XY} to true event (cm)", data.iteration), 120, -TMath::Pi(), TMath::Pi(), 120, 0., 3.); + data.hDZVsTgl = new TH2D(Form("%s_dz_vs_tgl", prefix.Data()), Form("Iteration %d truth matched;tan(#lambda);d_{Z} to true event (cm)", data.iteration), 120, -20., 20., 120, -30., 30.); + data.hDRAll = new TH1D(Form("%s_dr_all", prefix.Data()), Form("Iteration %d;#Deltar (cm);tracklets", data.iteration), 160, -5., 40.); + data.hDROk = new TH1D(Form("%s_dr_ok", prefix.Data()), Form("Iteration %d truth matched;#Deltar (cm);tracklets", data.iteration), 160, -5., 40.); + data.hDZGeomAll = new TH1D(Form("%s_dz_geom_all", prefix.Data()), Form("Iteration %d;#Deltaz (cm);tracklets", data.iteration), 240, -30., 30.); + data.hDZGeomOk = new TH1D(Form("%s_dz_geom_ok", prefix.Data()), Form("Iteration %d truth matched;#Deltaz (cm);tracklets", data.iteration), 240, -30., 30.); + data.hDPhiAll = new TH1D(Form("%s_dphi_all", prefix.Data()), Form("Iteration %d;#Delta#phi;tracklets", data.iteration), 160, 0., 0.2); + data.hDPhiOk = new TH1D(Form("%s_dphi_ok", prefix.Data()), Form("Iteration %d truth matched;#Delta#phi;tracklets", data.iteration), 160, 0., 0.2); + data.hTglEventOk = new TH1D(Form("%s_tgl_event_ok", prefix.Data()), Form("Iteration %d truth matched;tan(#lambda) from true event;tracklets", data.iteration), 200, -20., 20.); + data.hDeltaZEventOk = new TH1D(Form("%s_delta_z_event_ok", prefix.Data()), Form("Iteration %d truth matched;#Deltaz from true event (cm);tracklets", data.iteration), 240, 0., 30.); + data.hDeltaZEventPrim = new TH1D(Form("%s_delta_z_event_prim", prefix.Data()), Form("Iteration %d primary truth;#Deltaz from true event (cm);tracklets", data.iteration), 240, 0., 30.); + + for (auto* hist : {data.hPhiAll, data.hPhiOk, data.hTglAll, data.hTglOk, data.hDXYOk, data.hDXYPrim, data.hDZOk, data.hDZPrim, data.hDRAll, data.hDROk, data.hDZGeomAll, data.hDZGeomOk, data.hDPhiAll, data.hDPhiOk, data.hTglEventOk, data.hDeltaZEventOk, data.hDeltaZEventPrim}) { + hist->Sumw2(); + } +} + +TransitionPlots& getTransitionPlots(std::map, TransitionPlots>& plots, const std::pair& transition) +{ + auto& transitionPlots = plots[transition]; + if (transitionPlots.hDPhiAll) { + return transitionPlots; + } + + const TString prefix = Form("L%d_L%d", transition.first, transition.second); + transitionPlots.hDPhiAll = new TH1D(Form("%s_dphi_all", prefix.Data()), Form("L%d-L%d;#Delta#phi;tracklets", transition.first, transition.second), 160, 0., 0.2); + transitionPlots.hDPhiOk = new TH1D(Form("%s_dphi_ok", prefix.Data()), Form("L%d-L%d truth matched;#Delta#phi;tracklets", transition.first, transition.second), 160, 0., 0.2); + transitionPlots.hDeltaZEventOk = new TH1D(Form("%s_delta_z_event_ok", prefix.Data()), Form("L%d-L%d truth matched;#Deltaz from true event (cm);tracklets", transition.first, transition.second), 240, 0., 30.); + transitionPlots.hDeltaZEventPrim = new TH1D(Form("%s_delta_z_event_prim", prefix.Data()), Form("L%d-L%d primary truth;#Deltaz from true event (cm);tracklets", transition.first, transition.second), 240, 0., 30.); + + for (auto* hist : {transitionPlots.hDPhiAll, transitionPlots.hDPhiOk, transitionPlots.hDeltaZEventOk, transitionPlots.hDeltaZEventPrim}) { + hist->Sumw2(); + } + return transitionPlots; +} + +TransitionIterationPlots& getTransitionIterationPlots(std::map, std::map>& plots, const std::pair& transition, const int iteration) +{ + auto& iterationPlots = plots[transition][iteration]; + if (iterationPlots.hPhiAll) { + return iterationPlots; + } + + const TString prefix = Form("L%d_L%d_iter%d", transition.first, transition.second, iteration); + const TString title = Form("L%d-L%d iteration %d", transition.first, transition.second, iteration); + iterationPlots.hPhiAll = new TH1D(Form("%s_phi_all", prefix.Data()), Form("%s;#phi;tracklets", title.Data()), 160, -TMath::Pi(), TMath::Pi()); + iterationPlots.hPhiOk = new TH1D(Form("%s_phi_ok", prefix.Data()), Form("%s truth matched;#phi;tracklets", title.Data()), 160, -TMath::Pi(), TMath::Pi()); + iterationPlots.hTglAll = new TH1D(Form("%s_tgl_all", prefix.Data()), Form("%s;tan(#lambda);tracklets", title.Data()), 200, -20., 20.); + iterationPlots.hTglOk = new TH1D(Form("%s_tgl_ok", prefix.Data()), Form("%s truth matched;tan(#lambda);tracklets", title.Data()), 200, -20., 20.); + iterationPlots.hDXYOk = new TH1D(Form("%s_dxy_ok", prefix.Data()), Form("%s truth matched;d_{XY} to true event (cm);tracklets", title.Data()), 200, 0., 3.); + iterationPlots.hDXYPrim = new TH1D(Form("%s_dxy_prim", prefix.Data()), Form("%s primary truth;d_{XY} to true event (cm);tracklets", title.Data()), 200, 0., 3.); + iterationPlots.hDZOk = new TH1D(Form("%s_dz_ok", prefix.Data()), Form("%s truth matched;d_{Z} to true event (cm);tracklets", title.Data()), 240, -30., 30.); + iterationPlots.hDZPrim = new TH1D(Form("%s_dz_prim", prefix.Data()), Form("%s primary truth;d_{Z} to true event (cm);tracklets", title.Data()), 240, -30., 30.); + iterationPlots.hDRAll = new TH1D(Form("%s_dr_all", prefix.Data()), Form("%s;#Deltar (cm);tracklets", title.Data()), 160, -5., 40.); + iterationPlots.hDROk = new TH1D(Form("%s_dr_ok", prefix.Data()), Form("%s truth matched;#Deltar (cm);tracklets", title.Data()), 160, -5., 40.); + iterationPlots.hDZGeomAll = new TH1D(Form("%s_dz_geom_all", prefix.Data()), Form("%s;#Deltaz (cm);tracklets", title.Data()), 240, -30., 30.); + iterationPlots.hDZGeomOk = new TH1D(Form("%s_dz_geom_ok", prefix.Data()), Form("%s truth matched;#Deltaz (cm);tracklets", title.Data()), 240, -30., 30.); + iterationPlots.hDPhiAll = new TH1D(Form("%s_dphi_all", prefix.Data()), Form("%s;#Delta#phi;tracklets", title.Data()), 160, 0., 0.2); + iterationPlots.hDPhiOk = new TH1D(Form("%s_dphi_ok", prefix.Data()), Form("%s truth matched;#Delta#phi;tracklets", title.Data()), 160, 0., 0.2); + iterationPlots.hTglEventOk = new TH1D(Form("%s_tgl_event_ok", prefix.Data()), Form("%s truth matched;tan(#lambda) from true event;tracklets", title.Data()), 200, -20., 20.); + iterationPlots.hDeltaZEventOk = new TH1D(Form("%s_delta_z_event_ok", prefix.Data()), Form("%s truth matched;#Deltaz from true event (cm);tracklets", title.Data()), 240, 0., 30.); + iterationPlots.hDeltaZEventPrim = new TH1D(Form("%s_delta_z_event_prim", prefix.Data()), Form("%s primary truth;#Deltaz from true event (cm);tracklets", title.Data()), 240, 0., 30.); + + for (auto* hist : {iterationPlots.hPhiAll, iterationPlots.hPhiOk, iterationPlots.hTglAll, iterationPlots.hTglOk, iterationPlots.hDXYOk, iterationPlots.hDXYPrim, iterationPlots.hDZOk, iterationPlots.hDZPrim, iterationPlots.hDRAll, iterationPlots.hDROk, iterationPlots.hDZGeomAll, iterationPlots.hDZGeomOk, iterationPlots.hDPhiAll, iterationPlots.hDPhiOk, iterationPlots.hTglEventOk, iterationPlots.hDeltaZEventOk, iterationPlots.hDeltaZEventPrim}) { + hist->Sumw2(); + } + return iterationPlots; +} + +void setLine(TH1* hist, int color, int style = 1) +{ + hist->SetLineColor(color); + hist->SetMarkerColor(color); + hist->SetLineStyle(style); + hist->SetLineWidth(2); +} + +TH1D* normalizedClone(const TH1D* source, const char* name) +{ + auto* clone = static_cast(source->Clone(name)); + clone->SetDirectory(nullptr); + if (clone->Integral() > 0.) { + clone->Scale(1. / clone->Integral()); + } + return clone; +} + +void drawPair(TH1D* all, TH1D* selected, const char* allLabel, const char* selectedLabel) +{ + setLine(all, kBlack); + setLine(selected, kRed + 1); + all->Draw("hist"); + selected->Draw("hist same"); + auto* legend = new TLegend(0.58, 0.74, 0.88, 0.88); + legend->SetBorderSize(0); + legend->SetFillStyle(0); + legend->AddEntry(all, Form("%s: %.0f", allLabel, all->GetEntries()), "l"); + legend->AddEntry(selected, Form("%s: %.0f", selectedLabel, selected->GetEntries()), "l"); + legend->Draw(); +} + +void drawSingle(TH1D* hist, int color) +{ + setLine(hist, color); + hist->Draw("hist"); +} + +void drawIterationBreakdown(const std::map& plots, TH1D* TransitionIterationPlots::* member, const char* title) +{ + static const std::array colors = {kBlack, kRed + 1, kBlue + 1, kGreen + 2, kMagenta + 1, kOrange + 7, kCyan + 2, kViolet + 5, kAzure + 7, kGray + 2}; + bool first = true; + auto* legend = new TLegend(0.62, 0.62, 0.88, 0.88); + legend->SetBorderSize(0); + legend->SetFillStyle(0); + + size_t i = 0; + for (const auto& [iteration, iterationPlots] : plots) { + auto* hist = iterationPlots.*member; + setLine(hist, colors[i % colors.size()]); + hist->SetTitle(title); + hist->Draw(first ? "hist" : "hist same"); + first = false; + legend->AddEntry(hist, Form("iteration %d: %.0f", iteration, hist->GetEntries()), "l"); + ++i; + } + legend->Draw(); +} + +void drawOverlay(const std::vector& iterations, bool drawDXY) +{ + static const std::array colors = {kBlack, kRed + 1, kBlue + 1, kGreen + 2, kMagenta + 1, kOrange + 7, kCyan + 2, kViolet + 5, kAzure + 7, kGray + 2}; + bool first = true; + auto* legend = new TLegend(0.62, 0.62, 0.88, 0.88); + legend->SetBorderSize(0); + legend->SetFillStyle(0); + + for (size_t i = 0; i < iterations.size(); ++i) { + const auto* source = drawDXY ? iterations[i].hDXYOk : iterations[i].hDZOk; + auto* hist = normalizedClone(source, Form("%s_norm", source->GetName())); + setLine(hist, colors[i % colors.size()]); + hist->SetTitle(drawDXY ? "Truth matched d_{XY};d_{XY} to true event (cm);normalized tracklets" : "Truth matched d_{Z};d_{Z} to true event (cm);normalized tracklets"); + hist->Draw(first ? "hist" : "hist same"); + first = false; + legend->AddEntry(hist, Form("iteration %d", iterations[i].iteration), "l"); + } + legend->Draw(); +} + +void writeIterationHistograms(TDirectory* directory, const IterationData& data) +{ + directory->cd(); + data.hPhiAll->Write(); + data.hPhiOk->Write(); + data.hTglAll->Write(); + data.hTglOk->Write(); + data.hDXYOk->Write(); + data.hDXYPrim->Write(); + data.hDZOk->Write(); + data.hDZPrim->Write(); + data.hDRAll->Write(); + data.hDROk->Write(); + data.hDZGeomAll->Write(); + data.hDZGeomOk->Write(); + data.hDPhiAll->Write(); + data.hDPhiOk->Write(); + data.hTglEventOk->Write(); + data.hDeltaZEventOk->Write(); + data.hDeltaZEventPrim->Write(); + data.hDXYVsPhi->Write(); + data.hDZVsTgl->Write(); +} + +void writeTransitionIterationHistograms(TDirectory* directory, const TransitionIterationPlots& plots) +{ + directory->cd(); + plots.hPhiAll->Write(); + plots.hPhiOk->Write(); + plots.hTglAll->Write(); + plots.hTglOk->Write(); + plots.hDXYOk->Write(); + plots.hDXYPrim->Write(); + plots.hDZOk->Write(); + plots.hDZPrim->Write(); + plots.hDRAll->Write(); + plots.hDROk->Write(); + plots.hDZGeomAll->Write(); + plots.hDZGeomOk->Write(); + plots.hDPhiAll->Write(); + plots.hDPhiOk->Write(); + plots.hTglEventOk->Write(); + plots.hDeltaZEventOk->Write(); + plots.hDeltaZEventPrim->Write(); +} + +} // namespace + +void CheckTracklets(std::string inputFile = "its_tune.root", + std::string outputRoot = "CheckTracklets.root", + std::string outputPdf = "CheckTracklets.pdf") +{ + TH1::AddDirectory(kFALSE); + gStyle->SetOptStat(10); + + std::unique_ptr input(TFile::Open(inputFile.c_str(), "READ")); + if (!input || input->IsZombie()) { + std::cerr << "CheckTracklets: cannot open input file '" << inputFile << "'" << std::endl; + return; + } + + std::vector iterations; + TIter nextKey(input->GetListOfKeys()); + while (auto* key = static_cast(nextKey())) { + const TString name = key->GetName(); + int iteration = -1; + if (!parseIteration(name, iteration)) { + continue; + } + auto* klass = gROOT->GetClass(key->GetClassName()); + if (!klass || !klass->InheritsFrom(TTree::Class())) { + continue; + } + auto* tree = input->Get(name); + if (!tree) { + continue; + } + iterations.push_back({iteration, tree}); + } + + std::sort(iterations.begin(), iterations.end(), [](const auto& lhs, const auto& rhs) { + return lhs.iteration < rhs.iteration; + }); + + if (iterations.empty()) { + std::cerr << "CheckTracklets: no trklt_ trees found in '" << inputFile << "'" << std::endl; + return; + } + + std::set> transitions; + std::map, TransitionPlots> transitionPlots; + std::map, std::map> transitionIterationPlots; + for (auto& data : iterations) { + FormulaSet formulas; + if (!makeFormulas(data.tree, formulas)) { + return; + } + bookIterationHistograms(data); + const auto entries = data.tree->GetEntries(); + for (Long64_t entry = 0; entry < entries; ++entry) { + data.tree->LoadTree(entry); + data.tree->GetEntry(entry); + updateFormulas(formulas); + + const int from = static_cast(formulas.from->EvalInstance()); + const int to = static_cast(formulas.to->EvalInstance()); + const float phi = static_cast(formulas.phi->EvalInstance()); + const float tgl = static_cast(formulas.tgl->EvalInstance()); + const float dr = static_cast(formulas.dr->EvalInstance()); + const float dz = static_cast(formulas.dz->EvalInstance()); + const float dPhi = static_cast(formulas.dPhi->EvalInstance()); + const bool ok = formulas.ok->EvalInstance() != 0.; + const bool prim = formulas.prim->EvalInstance() != 0.; + const float dXY = static_cast(formulas.dXY->EvalInstance()); + const float dZ = static_cast(formulas.dZ->EvalInstance()); + const float tglEvent = static_cast(formulas.tglEvent->EvalInstance()); + const float deltaZEvent = static_cast(formulas.deltaZEvent->EvalInstance()); + + const auto transition = std::make_pair(from, to); + auto& stats = data.transitionStats[transition]; + auto& plots = getTransitionPlots(transitionPlots, transition); + auto& perIterationPlots = getTransitionIterationPlots(transitionIterationPlots, transition, data.iteration); + ++stats.all; + transitions.insert(transition); + + data.hPhiAll->Fill(phi); + data.hTglAll->Fill(tgl); + data.hDRAll->Fill(dr); + data.hDZGeomAll->Fill(dz); + data.hDPhiAll->Fill(dPhi); + plots.hDPhiAll->Fill(dPhi); + perIterationPlots.hPhiAll->Fill(phi); + perIterationPlots.hTglAll->Fill(tgl); + perIterationPlots.hDRAll->Fill(dr); + perIterationPlots.hDZGeomAll->Fill(dz); + perIterationPlots.hDPhiAll->Fill(dPhi); + + if (!ok) { + continue; + } + ++stats.ok; + data.hPhiOk->Fill(phi); + data.hTglOk->Fill(tgl); + data.hDXYOk->Fill(dXY); + data.hDZOk->Fill(dZ); + data.hDROk->Fill(dr); + data.hDZGeomOk->Fill(dz); + data.hDPhiOk->Fill(dPhi); + data.hTglEventOk->Fill(tglEvent); + data.hDeltaZEventOk->Fill(deltaZEvent); + data.hDXYVsPhi->Fill(phi, dXY); + data.hDZVsTgl->Fill(tgl, dZ); + plots.hDPhiOk->Fill(dPhi); + plots.hDeltaZEventOk->Fill(deltaZEvent); + perIterationPlots.hPhiOk->Fill(phi); + perIterationPlots.hTglOk->Fill(tgl); + perIterationPlots.hDXYOk->Fill(dXY); + perIterationPlots.hDZOk->Fill(dZ); + perIterationPlots.hDROk->Fill(dr); + perIterationPlots.hDZGeomOk->Fill(dz); + perIterationPlots.hDPhiOk->Fill(dPhi); + perIterationPlots.hTglEventOk->Fill(tglEvent); + perIterationPlots.hDeltaZEventOk->Fill(deltaZEvent); + if (prim) { + ++stats.prim; + data.hDXYPrim->Fill(dXY); + data.hDZPrim->Fill(dZ); + data.hDeltaZEventPrim->Fill(deltaZEvent); + plots.hDeltaZEventPrim->Fill(deltaZEvent); + perIterationPlots.hDXYPrim->Fill(dXY); + perIterationPlots.hDZPrim->Fill(dZ); + perIterationPlots.hDeltaZEventPrim->Fill(deltaZEvent); + } + } + } + + if (transitions.empty()) { + std::cerr << "CheckTracklets: no tracklet entries found in '" << inputFile << "'" << std::endl; + return; + } + + auto* hCounts = new TH2D("hTrackletCounts", "Tracklet counts;iteration;transition", iterations.size(), 0, iterations.size(), transitions.size(), 0, transitions.size()); + auto* hTruthFraction = new TH2D("hTruthMatchedFraction", "Truth matched fraction;iteration;transition", iterations.size(), 0, iterations.size(), transitions.size(), 0, transitions.size()); + auto* hPrimaryFraction = new TH2D("hPrimaryFraction", "Primary fraction among truth matched;iteration;transition", iterations.size(), 0, iterations.size(), transitions.size(), 0, transitions.size()); + + int ybin = 1; + for (const auto& transition : transitions) { + const auto label = Form("L%d-L%d", transition.first, transition.second); + hCounts->GetYaxis()->SetBinLabel(ybin, label); + hTruthFraction->GetYaxis()->SetBinLabel(ybin, label); + hPrimaryFraction->GetYaxis()->SetBinLabel(ybin, label); + ++ybin; + } + + for (size_t i = 0; i < iterations.size(); ++i) { + const int xbin = static_cast(i) + 1; + const auto label = Form("%d", iterations[i].iteration); + hCounts->GetXaxis()->SetBinLabel(xbin, label); + hTruthFraction->GetXaxis()->SetBinLabel(xbin, label); + hPrimaryFraction->GetXaxis()->SetBinLabel(xbin, label); + + ybin = 1; + for (const auto& transition : transitions) { + const auto found = iterations[i].transitionStats.find(transition); + if (found != iterations[i].transitionStats.end()) { + const auto& stats = found->second; + hCounts->SetBinContent(xbin, ybin, stats.all); + if (stats.all > 0) { + hTruthFraction->SetBinContent(xbin, ybin, static_cast(stats.ok) / static_cast(stats.all)); + } + if (stats.ok > 0) { + hPrimaryFraction->SetBinContent(xbin, ybin, static_cast(stats.prim) / static_cast(stats.ok)); + } + } + ++ybin; + } + } + + for (auto& [transition, plots] : transitionPlots) { + plots.hDPhiTruthFraction = static_cast(plots.hDPhiOk->Clone(Form("L%d_L%d_dphi_truth_fraction", transition.first, transition.second))); + plots.hDPhiTruthFraction->SetTitle(Form("L%d-L%d truth matched fraction;#Delta#phi;ok / all", transition.first, transition.second)); + plots.hDPhiTruthFraction->Divide(plots.hDPhiOk, plots.hDPhiAll, 1., 1., "B"); + plots.hDPhiTruthFraction->SetMinimum(0.); + plots.hDPhiTruthFraction->SetMaximum(1.05); + } + + std::unique_ptr output(TFile::Open(outputRoot.c_str(), "RECREATE")); + if (!output || output->IsZombie()) { + std::cerr << "CheckTracklets: cannot create output file '" << outputRoot << "'" << std::endl; + return; + } + + auto* summaryDir = output->mkdir("summary"); + summaryDir->cd(); + hCounts->Write(); + hTruthFraction->Write(); + hPrimaryFraction->Write(); + for (const auto& [transition, plots] : transitionPlots) { + auto* transitionDir = summaryDir->mkdir(Form("transition_L%d_L%d", transition.first, transition.second)); + transitionDir->cd(); + plots.hDPhiAll->Write(); + plots.hDPhiOk->Write(); + plots.hDPhiTruthFraction->Write(); + plots.hDeltaZEventOk->Write(); + plots.hDeltaZEventPrim->Write(); + const auto transitionIterations = transitionIterationPlots.find(transition); + if (transitionIterations != transitionIterationPlots.end()) { + auto* iterationsDir = transitionDir->mkdir("iterations"); + for (const auto& [iteration, iterationPlots] : transitionIterations->second) { + auto* iterationDir = iterationsDir->mkdir(Form("iteration_%d", iteration)); + writeTransitionIterationHistograms(iterationDir, iterationPlots); + } + } + } + for (const auto& data : iterations) { + auto* iterDir = output->mkdir(Form("iteration_%d", data.iteration)); + writeIterationHistograms(iterDir, data); + } + output->Close(); + + TCanvas canvas("cCheckTracklets", "ITS tune tracklets", 1800, 1100); + canvas.Print((outputPdf + "[").c_str()); + + canvas.Clear(); + canvas.Divide(2, 2); + canvas.cd(1); + gPad->SetRightMargin(0.16); + hCounts->Draw("colz text"); + canvas.cd(2); + gPad->SetRightMargin(0.16); + hTruthFraction->SetMinimum(0.); + hTruthFraction->SetMaximum(1.); + hTruthFraction->Draw("colz text"); + canvas.cd(3); + gPad->SetRightMargin(0.16); + hPrimaryFraction->SetMinimum(0.); + hPrimaryFraction->SetMaximum(1.); + hPrimaryFraction->Draw("colz text"); + canvas.cd(4); + drawOverlay(iterations, true); + canvas.Print(outputPdf.c_str()); + + canvas.Clear(); + canvas.Divide(2, 1); + canvas.cd(1); + drawOverlay(iterations, true); + canvas.cd(2); + drawOverlay(iterations, false); + canvas.Print(outputPdf.c_str()); + + for (const auto& [transition, plots] : transitionPlots) { + canvas.Clear(); + canvas.Divide(2, 1); + canvas.cd(1); + drawSingle(plots.hDPhiTruthFraction, kBlue + 1); + canvas.cd(2); + drawPair(plots.hDeltaZEventOk, plots.hDeltaZEventPrim, "truth", "primary"); + canvas.Print(outputPdf.c_str()); + } + + for (const auto& [transition, plots] : transitionIterationPlots) { + canvas.Clear(); + canvas.Divide(3, 2); + canvas.cd(1); + drawIterationBreakdown(plots, &TransitionIterationPlots::hPhiOk, Form("L%d-L%d truth matched by iteration;#phi;tracklets", transition.first, transition.second)); + canvas.cd(2); + drawIterationBreakdown(plots, &TransitionIterationPlots::hTglOk, Form("L%d-L%d truth matched by iteration;tan(#lambda);tracklets", transition.first, transition.second)); + canvas.cd(3); + drawIterationBreakdown(plots, &TransitionIterationPlots::hDXYOk, Form("L%d-L%d truth matched by iteration;d_{XY} to true event (cm);tracklets", transition.first, transition.second)); + canvas.cd(4); + drawIterationBreakdown(plots, &TransitionIterationPlots::hDZOk, Form("L%d-L%d truth matched by iteration;d_{Z} to true event (cm);tracklets", transition.first, transition.second)); + canvas.cd(5); + drawIterationBreakdown(plots, &TransitionIterationPlots::hDPhiOk, Form("L%d-L%d truth matched by iteration;#Delta#phi;tracklets", transition.first, transition.second)); + canvas.cd(6); + drawIterationBreakdown(plots, &TransitionIterationPlots::hDeltaZEventOk, Form("L%d-L%d truth matched by iteration;#Deltaz from true event (cm);tracklets", transition.first, transition.second)); + canvas.Print(outputPdf.c_str()); + + canvas.Clear(); + canvas.Divide(3, 2); + canvas.cd(1); + drawIterationBreakdown(plots, &TransitionIterationPlots::hDRAll, Form("L%d-L%d all tracklets by iteration;#Deltar (cm);tracklets", transition.first, transition.second)); + canvas.cd(2); + drawIterationBreakdown(plots, &TransitionIterationPlots::hDROk, Form("L%d-L%d truth matched by iteration;#Deltar (cm);tracklets", transition.first, transition.second)); + canvas.cd(3); + drawIterationBreakdown(plots, &TransitionIterationPlots::hDZGeomAll, Form("L%d-L%d all tracklets by iteration;#Deltaz (cm);tracklets", transition.first, transition.second)); + canvas.cd(4); + drawIterationBreakdown(plots, &TransitionIterationPlots::hDZGeomOk, Form("L%d-L%d truth matched by iteration;#Deltaz (cm);tracklets", transition.first, transition.second)); + canvas.cd(5); + drawIterationBreakdown(plots, &TransitionIterationPlots::hDPhiAll, Form("L%d-L%d all tracklets by iteration;#Delta#phi;tracklets", transition.first, transition.second)); + canvas.cd(6); + drawIterationBreakdown(plots, &TransitionIterationPlots::hTglEventOk, Form("L%d-L%d truth matched by iteration;tan(#lambda) from true event;tracklets", transition.first, transition.second)); + canvas.Print(outputPdf.c_str()); + } + + for (auto& data : iterations) { + canvas.Clear(); + canvas.Divide(3, 2); + canvas.cd(1); + drawPair(data.hPhiAll, data.hPhiOk, "all", "truth"); + canvas.cd(2); + drawPair(data.hTglAll, data.hTglOk, "all", "truth"); + canvas.cd(3); + drawPair(data.hDXYOk, data.hDXYPrim, "truth", "primary"); + canvas.cd(4); + drawPair(data.hDZOk, data.hDZPrim, "truth", "primary"); + canvas.cd(5); + gPad->SetRightMargin(0.16); + data.hDXYVsPhi->Draw("colz"); + canvas.cd(6); + gPad->SetRightMargin(0.16); + data.hDZVsTgl->Draw("colz"); + canvas.Print(outputPdf.c_str()); + + canvas.Clear(); + canvas.Divide(3, 2); + canvas.cd(1); + drawPair(data.hDRAll, data.hDROk, "all", "truth"); + canvas.cd(2); + drawPair(data.hDZGeomAll, data.hDZGeomOk, "all", "truth"); + canvas.cd(3); + drawPair(data.hDPhiAll, data.hDPhiOk, "all", "truth"); + canvas.cd(4); + drawPair(data.hDeltaZEventOk, data.hDeltaZEventPrim, "truth", "primary"); + canvas.cd(5); + drawSingle(data.hTglEventOk, kBlue + 1); + canvas.Print(outputPdf.c_str()); + } + + canvas.Print((outputPdf + "]").c_str()); + + std::cout << "CheckTracklets: wrote " << outputRoot << " and " << outputPdf << std::endl; +} diff --git a/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt b/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt index 1dd64b6f1874b..0029eb20624a9 100644 --- a/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt @@ -21,9 +21,11 @@ o2_add_library(ITStracking src/IOUtils.cxx src/Tracker.cxx src/TrackerTraits.cxx + src/TrackerTraitsMC.cxx src/TrackingConfigParam.cxx src/Vertexer.cxx src/VertexerTraits.cxx + src/TuneExt.cxx PUBLIC_LINK_LIBRARIES O2::GPUCommon Microsoft.GSL::GSL @@ -55,6 +57,7 @@ o2_target_root_dictionary(ITStracking include/ITStracking/Definitions.h include/ITStracking/FastMultEstConfig.h include/ITStracking/TrackingConfigParam.h + include/ITStracking/TuneExt.h LINKDEF src/TrackingLinkDef.h) if(CUDA_ENABLED OR HIP_ENABLED) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h index f536e86fe95d5..4f26c545d6d83 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h @@ -85,11 +85,12 @@ class TrackerTraits std::shared_ptr mTaskArena; protected: - o2::gpu::GPUChainITS* mChain = nullptr; - TimeFrame* mTimeFrame; - std::vector mTrkParams; + void createTrackletMC(); - float mBz{-999.f}; + o2::gpu::GPUChainITS* mChain{nullptr}; + TimeFrame* mTimeFrame{nullptr}; + std::vector mTrkParams; + float mBz{constants::UnsetValue}; }; } // namespace its diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraitsMC.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraitsMC.h new file mode 100644 index 0000000000000..92e236768fb9a --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraitsMC.h @@ -0,0 +1,44 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// + +#ifndef TRACKINGITSU_INCLUDE_TRACKERTRAITSMC_H_ +#define TRACKINGITSU_INCLUDE_TRACKERTRAITSMC_H_ + +#include "ITStracking/TrackerTraits.h" + +namespace o2::its +{ + +// This class faciliates the tuning of the reconstruction parameters. +// Two modes are foreseen: +// 1. Inspect and dump the artefacts the reconstrcution produces +// 2. Run the same over MC truth +// Both should faciliate finding sources of in-efficiency (c.f., missing links) +// and allow to tune overall the imposed parameters in a consistent way. +template +class TrackerTraitsMC : public TrackerTraits +{ + public: + TrackerTraitsMC() = default; + ~TrackerTraitsMC() override = default; + + void computeLayerTracklets(const int iteration, int iVertex) override; + void computeLayerCells(const int iteration) override; + void findCellsNeighbours(const int iteration) override; + void findRoads(const int iteration) override; + + const char* getName() const noexcept override { return "TUNE"; } +}; + +} // namespace o2::its + +#endif /* TRACKINGITSU_INCLUDE_TRACKERTRAITSMC_H_ */ diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TuneExt.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TuneExt.h new file mode 100644 index 0000000000000..b30bd4ba9320c --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TuneExt.h @@ -0,0 +1,42 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// Holder of auxiallary information used for tuning +/// + +#include "ITStracking/Constants.h" +#include "GPUCommonRtypes.h" + +namespace o2::its +{ +struct TrackletMC final { + float tgl{constants::UnsetValue}; // tanLambda + float phi{constants::UnsetValue}; // phi + float rIn{constants::UnsetValue}; + float zIn{constants::UnsetValue}; + float phiIn{constants::UnsetValue}; + float rOut{constants::UnsetValue}; + float zOut{constants::UnsetValue}; + float phiOut{constants::UnsetValue}; + float dr{constants::UnsetValue}; + float dz{constants::UnsetValue}; + float dPhi{constants::UnsetValue}; + bool ok{false}; // truth + /// below only metrics valid if ok + bool prim{false}; // primary + float dXY{constants::UnsetValue}; // transverse distance to event + float dZ{constants::UnsetValue}; // longitudinal distance to event + float deltaZEvent{constants::UnsetValue}; + float tglEvent{constants::UnsetValue}; + ClassDefNV(TrackletMC, 2); +}; + +} // namespace o2::its diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx index c4439dc74d29e..f1e5de27595bf 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx @@ -223,26 +223,7 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iVer /// Create tracklets labels if (mTimeFrame->hasMCinformation() && mTrkParams[iteration].CreateArtefactLabels) { - tbb::parallel_for(0, static_cast(topology.nTransitions), [&](const int transitionId) { - const auto& transition = topology.getTransition(transitionId); - for (auto& trk : mTimeFrame->getTracklets()[transitionId]) { - MCCompLabel label; - int currentId{mTimeFrame->getClusters()[transition.fromLayer][trk.firstClusterIndex].clusterId}; - int nextId{mTimeFrame->getClusters()[transition.toLayer][trk.secondClusterIndex].clusterId}; - for (const auto& lab1 : mTimeFrame->getClusterLabels(transition.fromLayer, currentId)) { - for (const auto& lab2 : mTimeFrame->getClusterLabels(transition.toLayer, nextId)) { - if (lab1 == lab2 && lab1.isValid()) { - label = lab1; - break; - } - } - if (label.isValid()) { - break; - } - } - mTimeFrame->getTrackletsLabel(transitionId).emplace_back(label); - } - }); + createTrackletMC(); } }); } @@ -907,6 +888,32 @@ void TrackerTraits::markTracks(int iteration) } } +template +void TrackerTraits::createTrackletMC() +{ + const auto topology = mTimeFrame->getTrackingTopologyView(); + tbb::parallel_for(0, static_cast(topology.nTransitions), [&](const int transitionId) { + const auto& transition = topology.getTransition(transitionId); + for (auto& trk : mTimeFrame->getTracklets()[transitionId]) { + MCCompLabel label; + int currentId{mTimeFrame->getClusters()[transition.fromLayer][trk.firstClusterIndex].clusterId}; + int nextId{mTimeFrame->getClusters()[transition.toLayer][trk.secondClusterIndex].clusterId}; + for (const auto& lab1 : mTimeFrame->getClusterLabels(transition.fromLayer, currentId)) { + for (const auto& lab2 : mTimeFrame->getClusterLabels(transition.toLayer, nextId)) { + if (lab1 == lab2 && lab1.isValid()) { + label = lab1; + break; + } + } + if (label.isValid()) { + break; + } + } + mTimeFrame->getTrackletsLabel(transitionId).emplace_back(label); + } + }); +} + template void TrackerTraits::setBz(float bz) { diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraitsMC.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraitsMC.cxx new file mode 100644 index 0000000000000..e502566bfa5a7 --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraitsMC.cxx @@ -0,0 +1,128 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// + +#include +#include +#include +#include + +#include "ITStracking/TrackerTraitsMC.h" +#include "ITStracking/TuneExt.h" +#include "CommonUtils/TreeStreamRedirector.h" +#include "CommonConstants/MathConstants.h" +#include "Steer/MCKinematicsReader.h" + +#include "Framework/Logger.h" + +namespace o2::its +{ + +static o2::utils::TreeStreamRedirector* gDBG{nullptr}; +static o2::steer::MCKinematicsReader* gMCReader{nullptr}; + +template +void TrackerTraitsMC::computeLayerTracklets(const int iteration, int iVertex) +{ + if (!gDBG) { + gDBG = new o2::utils::TreeStreamRedirector("its_tune.root"); + } + if (!gMCReader) { + gMCReader = new o2::steer::MCKinematicsReader("collisioncontext.root"); + } + + // Create all tracklets we find in this iteration and dump them. + const std::string treeName = std::format("trklt_{}", iteration); + TrackerTraits::computeLayerTracklets(iteration, iVertex); + this->createTrackletMC(); + const auto topology = this->mTimeFrame->getTrackingTopologyView(); + for (int transitionId{0}; transitionId < topology.nTransitions; ++transitionId) { + const auto& transition = topology.getTransition(transitionId); + for (int iTrklt{0}; iTrklt < this->mTimeFrame->getTracklets()[transitionId].size(); ++iTrklt) { + const auto& lbl = this->mTimeFrame->getTrackletsLabel(transitionId)[iTrklt]; + const auto trklt = this->mTimeFrame->getTracklets()[transitionId][iTrklt]; + const auto& firstCluster = this->mTimeFrame->getClusters()[transition.fromLayer][trklt.firstClusterIndex]; + const auto& secondCluster = this->mTimeFrame->getClusters()[transition.toLayer][trklt.secondClusterIndex]; + const float deltaPhi = std::abs(firstCluster.phi - secondCluster.phi); + TrackletMC trkltMC{ + .tgl = trklt.tanLambda, + .phi = trklt.phi, + .rIn = firstCluster.radius, + .zIn = firstCluster.zCoordinate, + .phiIn = firstCluster.phi, + .rOut = secondCluster.radius, + .zOut = secondCluster.zCoordinate, + .phiOut = secondCluster.phi, + .dr = secondCluster.radius - firstCluster.radius, + .dz = secondCluster.zCoordinate - firstCluster.zCoordinate, + .dPhi = std::min(deltaPhi, static_cast(o2::constants::math::TwoPI) - deltaPhi), + .ok = lbl.isValid(), + }; + float dcaXY = std::numeric_limits::max(), dcaZ = std::numeric_limits::max(); + if (lbl.isValid()) { + const auto& eve = gMCReader->getMCEventHeader(lbl.getSourceID(), lbl.getEventID()); + const float dx = secondCluster.xCoordinate - firstCluster.xCoordinate; + const float dy = secondCluster.yCoordinate - firstCluster.yCoordinate; + const float dz = secondCluster.zCoordinate - firstCluster.zCoordinate; + trkltMC.tglEvent = (firstCluster.zCoordinate - eve.GetZ()) / firstCluster.radius; + trkltMC.deltaZEvent = std::abs((trkltMC.tglEvent * (secondCluster.radius - firstCluster.radius)) + firstCluster.zCoordinate - secondCluster.zCoordinate); + const float dxy2 = math_utils::hypot(dx, dy); + if (dxy2 > constants::Tolerance) { + const float t = ((eve.GetX() - firstCluster.xCoordinate) * dx + (eve.GetY() - firstCluster.yCoordinate) * dy) / dxy2; + const float xAtDCA = firstCluster.xCoordinate + t * dx; + const float yAtDCA = firstCluster.yCoordinate + t * dy; + const float zAtDCA = firstCluster.zCoordinate + t * dz; + const float curDCAx = xAtDCA - eve.GetX(); + const float curDCAy = yAtDCA - eve.GetY(); + dcaXY = math_utils::hypot(curDCAx, curDCAy); + dcaZ = zAtDCA - eve.GetZ(); + trkltMC.dXY = dcaXY; + trkltMC.dZ = dcaZ; + } + const auto* mcTrk = gMCReader->getTrack(lbl); + if (mcTrk) { + trkltMC.prim = mcTrk->isPrimary(); + } + } + (*gDBG) << treeName.c_str() + << "from=" << transition.fromLayer + << "to=" << transition.toLayer + << "trklt=" << trkltMC + << "\n"; + } + } +} + +template +void TrackerTraitsMC::computeLayerCells(const int iteration) +{ + TrackerTraits::computeLayerCells(iteration); +} + +template +void TrackerTraitsMC::findCellsNeighbours(const int iteration) +{ + TrackerTraits::findCellsNeighbours(iteration); +} + +template +void TrackerTraitsMC::findRoads(const int iteration) +{ + TrackerTraits::findRoads(iteration); + + if (this->mTrkParams.size() - 1 == iteration) { + gDBG->Close(); + } +} + +template class TrackerTraitsMC<7>; + +} // namespace o2::its diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackingLinkDef.h b/Detectors/ITSMFT/ITS/tracking/src/TrackingLinkDef.h index 46af692fe0c15..28423ddcf2154 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingLinkDef.h +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingLinkDef.h @@ -18,6 +18,9 @@ #pragma link C++ class o2::its::Tracklet + ; #pragma link C++ class std::vector < o2::its::Tracklet> + ; +#pragma link C++ class o2::its::TrackletMC + ; +#pragma link C++ class std::vector < o2::its::TrackletMC> + ; + #pragma link C++ class o2::its::Cluster + ; #pragma link C++ class std::vector < o2::its::Cluster> + ; diff --git a/Detectors/ITSMFT/ITS/tracking/src/TuneExt.cxx b/Detectors/ITSMFT/ITS/tracking/src/TuneExt.cxx new file mode 100644 index 0000000000000..3a7adcfdabc5c --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/src/TuneExt.cxx @@ -0,0 +1,13 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "ITStracking/TuneExt.h" +ClassImp(o2::its::TrackletMC); diff --git a/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/RecoWorkflow.h b/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/RecoWorkflow.h index 3068954c92003..1f9f6f69609b9 100644 --- a/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/RecoWorkflow.h +++ b/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/RecoWorkflow.h @@ -28,7 +28,7 @@ namespace reco_workflow framework::WorkflowSpec getWorkflow(bool useMC, bool doStag, TrackingMode::Type trmode, const bool overrideBeamPosition = false, bool upstreamDigits = false, bool upstreamClusters = false, bool clrofOnly = false, bool disableRootOutput = false, bool useGeom = false, int useTrig = 0, - bool useGPUWF = false, o2::gpu::gpudatatypes::DeviceType dType = o2::gpu::gpudatatypes::DeviceType::CPU); + bool useGPUWF = false, o2::gpu::gpudatatypes::DeviceType dType = o2::gpu::gpudatatypes::DeviceType::CPU, bool enableTuning = false); } } // namespace its diff --git a/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/TrackerSpec.h b/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/TrackerSpec.h index 8ce63efcb7a3b..4df986abaea7b 100644 --- a/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/TrackerSpec.h +++ b/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/TrackerSpec.h @@ -46,7 +46,8 @@ class TrackerDPL : public framework::Task int trgType, const TrackingMode::Type trMode = TrackingMode::Unset, const bool overrBeamEst = false, - o2::gpu::gpudatatypes::DeviceType dType = o2::gpu::gpudatatypes::DeviceType::CPU); + o2::gpu::gpudatatypes::DeviceType dType = o2::gpu::gpudatatypes::DeviceType::CPU, + bool enableTuning = false); ~TrackerDPL() override = default; void init(framework::InitContext& ic) final; void run(framework::ProcessingContext& pc) final; @@ -61,10 +62,11 @@ class TrackerDPL : public framework::Task std::unique_ptr mChainITS = nullptr; std::shared_ptr mGGCCDBRequest; ITSTrackingInterface mITSTrackingInterface; + bool mEnableTuning = false; TStopwatch mTimer; }; -framework::DataProcessorSpec getTrackerSpec(bool useMC, bool doStag, bool useGeom, int useTrig, TrackingMode::Type trMode, const bool overrBeamEst = false, o2::gpu::gpudatatypes::DeviceType dType = o2::gpu::gpudatatypes::DeviceType::CPU); +framework::DataProcessorSpec getTrackerSpec(bool useMC, bool doStag, bool useGeom, int useTrig, TrackingMode::Type trMode, const bool overrBeamEst = false, o2::gpu::gpudatatypes::DeviceType dType = o2::gpu::gpudatatypes::DeviceType::CPU, bool enableTuning = false); } // namespace o2::its diff --git a/Detectors/ITSMFT/ITS/workflow/src/RecoWorkflow.cxx b/Detectors/ITSMFT/ITS/workflow/src/RecoWorkflow.cxx index 06b3f019a6be7..f582003252fe5 100644 --- a/Detectors/ITSMFT/ITS/workflow/src/RecoWorkflow.cxx +++ b/Detectors/ITSMFT/ITS/workflow/src/RecoWorkflow.cxx @@ -37,7 +37,8 @@ framework::WorkflowSpec getWorkflow(bool useMC, bool doStag, bool useGeom, int useTrig, bool useGPUWF, - o2::gpu::gpudatatypes::DeviceType dtype) + o2::gpu::gpudatatypes::DeviceType dtype, + bool enableTuning) { framework::WorkflowSpec specs; if (!(upstreamDigits || upstreamClusters)) { @@ -80,7 +81,7 @@ framework::WorkflowSpec getWorkflow(bool useMC, bool doStag, .algorithm = AlgorithmSpec{adoptTask(task)}, .options = taskOptions}); } else { - specs.emplace_back(o2::its::getTrackerSpec(useMC, doStag, useGeom, useTrig, trmode, overrideBeamPosition, dtype)); + specs.emplace_back(o2::its::getTrackerSpec(useMC, doStag, useGeom, useTrig, trmode, overrideBeamPosition, dtype, enableTuning)); } if (!disableRootOutput) { specs.emplace_back(o2::its::getTrackWriterSpec(useMC)); diff --git a/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx b/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx index bbafc48e931ed..5680333e27bf9 100644 --- a/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx +++ b/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx @@ -11,14 +11,13 @@ #include -#include "Framework/ControlService.h" -#include "Framework/ConfigParamRegistry.h" #include "Framework/CCDBParamSpec.h" #include "Framework/DeviceSpec.h" +#include "Framework/Logger.h" #include "DataFormatsITSMFT/DPLAlpideParam.h" #include "ITSWorkflow/TrackerSpec.h" -#include "ITStracking/Definitions.h" #include "ITStracking/TrackingConfigParam.h" +#include "ITStracking/TrackerTraitsMC.h" namespace o2 { @@ -31,19 +30,24 @@ TrackerDPL::TrackerDPL(std::shared_ptr gr, int trgType, const TrackingMode::Type trMode, const bool overrBeamEst, - o2::gpu::gpudatatypes::DeviceType dType) : mGGCCDBRequest(gr), - mRecChain{o2::gpu::GPUReconstruction::CreateInstance(dType, true)}, - mITSTrackingInterface{isMC, doStag, trgType, overrBeamEst} + o2::gpu::gpudatatypes::DeviceType dType, + bool enableTuning) : mRecChain{o2::gpu::GPUReconstruction::CreateInstance(dType, true)}, + mGGCCDBRequest(gr), + mITSTrackingInterface{isMC, doStag, trgType, overrBeamEst}, + mEnableTuning{enableTuning} { mITSTrackingInterface.setTrackingMode(trMode); } -void TrackerDPL::init(InitContext& ic) +void TrackerDPL::init(InitContext& /*ic*/) { mTimer.Stop(); mTimer.Reset(); o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); mChainITS.reset(mRecChain->AddChain()); + if (mEnableTuning) { + mChainITS->SetITSTrackerTraits(std::make_unique>()); + } mITSTrackingInterface.setTraitsFromProvider(mChainITS->GetITSVertexerTraits(), mChainITS->GetITSTrackerTraits(), mChainITS->GetITSTimeframe()); @@ -78,7 +82,7 @@ void TrackerDPL::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) mITSTrackingInterface.finaliseCCDB(matcher, obj); } -void TrackerDPL::endOfStream(EndOfStreamContext& ec) +void TrackerDPL::endOfStream(EndOfStreamContext& /*ec*/) { end(); } @@ -92,7 +96,7 @@ void TrackerDPL::end() } } -DataProcessorSpec getTrackerSpec(bool useMC, bool doStag, bool useGeom, int trgType, TrackingMode::Type trMode, const bool overrBeamEst, o2::gpu::gpudatatypes::DeviceType dType) +DataProcessorSpec getTrackerSpec(bool useMC, bool doStag, bool useGeom, int trgType, TrackingMode::Type trMode, const bool overrBeamEst, o2::gpu::gpudatatypes::DeviceType dType, bool enableTuning) { const int mLayers = doStag ? o2::itsmft::DPLAlpideParam::getNLayers() : 1; std::vector inputs; @@ -149,7 +153,8 @@ DataProcessorSpec getTrackerSpec(bool useMC, bool doStag, bool useGeom, int trgT trgType, trMode, overrBeamEst, - dType)}, + dType, + enableTuning)}, .options = Options{}}; } diff --git a/Detectors/ITSMFT/ITS/workflow/src/its-reco-workflow.cxx b/Detectors/ITSMFT/ITS/workflow/src/its-reco-workflow.cxx index f1d60b8ac2c9b..001263e3d917e 100644 --- a/Detectors/ITSMFT/ITS/workflow/src/its-reco-workflow.cxx +++ b/Detectors/ITSMFT/ITS/workflow/src/its-reco-workflow.cxx @@ -13,6 +13,7 @@ #include "DataFormatsITSMFT/DPLAlpideParamInitializer.h" #include "CommonUtils/ConfigurableParam.h" #include "ITStracking/Configuration.h" +#include "ITStracking/TrackingConfigParam.h" #include "DetectorsRaw/HBFUtilsInitializer.h" #include "Framework/CallbacksPolicy.h" #include "Framework/ConfigContext.h" @@ -46,6 +47,7 @@ void customize(std::vector& workflowOptions) {"ccdb-meanvertex-seed", o2::framework::VariantType::Bool, false, {"use MeanVertex from CCDB if available to provide beam position seed (default: false)"}}, {"select-with-triggers", o2::framework::VariantType::String, "none", {"use triggers to prescale processed ROFs: phys, trd, none"}}, {"tracking-mode", o2::framework::VariantType::String, "sync", {"sync,async,cosmics,unset,off"}}, + {"enable-tuning", o2::framework::VariantType::Bool, false, {"enable MC-based CPU tracker tuning"}}, {"disable-tracking", o2::framework::VariantType::Bool, false, {"disable tracking step"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}, {"use-full-geometry", o2::framework::VariantType::Bool, false, {"use full geometry instead of the light-weight ITS part"}}, @@ -70,6 +72,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) auto trmode = configcontext.options().get("tracking-mode"); auto selTrig = configcontext.options().get("select-with-triggers"); auto useGpuWF = configcontext.options().get("use-gpu-workflow"); + auto enableTuning = configcontext.options().get("enable-tuning"); auto gpuDevice = static_cast(configcontext.options().get("gpu-device")); auto extDigits = configcontext.options().get("digits-from-upstream"); auto extClusters = configcontext.options().get("clusters-from-upstream"); @@ -92,10 +95,20 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) LOG(fatal) << "Unknown trigger type requested for events prescaling: " << selTrig; } } + const auto trackingMode = o2::its::TrackingMode::fromString(trmode); + if (enableTuning && !useMC) { + LOGP(fatal, "ITS tracker tuning requires MC labels"); + } + if (enableTuning && useGpuWF) { + LOGP(fatal, "ITS tracker tuning is CPU-only and cannot use the GPU workflow"); + } + if (enableTuning && gpuDevice != o2::gpu::gpudatatypes::DeviceType::CPU) { + LOGP(fatal, "ITS tracker tuning is CPU-only and cannot use GPU tracker traits"); + } auto wf = o2::its::reco_workflow::getWorkflow( useMC, doStag, - o2::its::TrackingMode::fromString(trmode), + trackingMode, beamPosOVerride, extDigits, extClusters, @@ -104,7 +117,8 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) useGeom, trType, useGpuWF, - gpuDevice); + gpuDevice, + enableTuning); // configure dpl timer to inject correct firstTForbit: start from the 1st orbit of TF containing 1st sampled orbit o2::raw::HBFUtilsInitializer hbfIni(configcontext, wf); diff --git a/GPU/GPUTracking/Global/GPUChainITS.cxx b/GPU/GPUTracking/Global/GPUChainITS.cxx index 598f7a61cac1a..2631df39a999d 100644 --- a/GPU/GPUTracking/Global/GPUChainITS.cxx +++ b/GPU/GPUTracking/Global/GPUChainITS.cxx @@ -59,6 +59,11 @@ o2::its::TrackerTraits<7>* GPUChainITS::GetITSTrackerTraits() return mITSTrackerTraits.get(); } +void GPUChainITS::SetITSTrackerTraits(std::unique_ptr> trackerTraits) +{ + mITSTrackerTraits = std::move(trackerTraits); +} + o2::its::VertexerTraits<7>* GPUChainITS::GetITSVertexerTraits() { if (mITSVertexerTraits == nullptr) { diff --git a/GPU/GPUTracking/Global/GPUChainITS.h b/GPU/GPUTracking/Global/GPUChainITS.h index ee466365a157d..3ef7c5b3ab83e 100644 --- a/GPU/GPUTracking/Global/GPUChainITS.h +++ b/GPU/GPUTracking/Global/GPUChainITS.h @@ -22,6 +22,8 @@ struct Cluster; struct TrackingFrameInfo; class TrackITSExt; class GPUFrameworkExternalAllocator; +template +class TrackerTraits; } // namespace o2::its namespace o2::gpu @@ -42,6 +44,7 @@ class GPUChainITS final : public GPUChain void MemorySize(size_t&, size_t&) final {}; o2::its::TrackerTraits<7>* GetITSTrackerTraits(); + void SetITSTrackerTraits(std::unique_ptr> trackerTraits); o2::its::VertexerTraits<7>* GetITSVertexerTraits(); o2::its::TimeFrame<7>* GetITSTimeframe();