[ALICE3] Refine pdgCodes and update histogram logic (#18205)

Author: njacazio

ALICE3/Tasks/alice3TrackingPerformance.cxx

@@ -34,43 +34,54 @@
 #include <TProfile2D.h>
 #include <TString.h>
 
+#include <algorithm>
 #include <map>
 #include <memory>
 #include <string>
 #include <vector>
 
 using namespace o2;
 using namespace o2::framework;
-std::shared_ptr<TH1> particlePdgCodes;
-std::map<int, std::shared_ptr<TH1>> particlePtDistribution;
-std::map<int, std::shared_ptr<TH1>> particleEtaDistribution;
-std::map<int, std::shared_ptr<TH1>> ptDistribution;
-std::map<int, std::shared_ptr<TH2>> ptResolutionVsPt;
-std::map<int, std::shared_ptr<TProfile2D>> ptResolutionVsEta;
-std::map<int, std::shared_ptr<TH2>> invPtResolutionVsPt;
-std::map<int, std::shared_ptr<TProfile2D>> invPtResolutionVsEta;
-std::map<int, std::shared_ptr<TH2>> dcaXyResolutionVsPt;
-std::map<int, std::shared_ptr<TH2>> dcaZResolutionVsPt;
-std::map<int, std::shared_ptr<TH2>> covariancePtPtVsPt;
-std::map<int, std::shared_ptr<TH2>> covarianceDcaXyDcaXyVsPt;
-std::map<int, std::shared_ptr<TH2>> covarianceDcaZDcaZVsPt;
 
 struct Alice3TrackingPerformance {
-  Configurable<std::vector<int>> pdgCodes{"pdgCodes", {0, 211}, "List of PDG codes to consider for efficiency calculation. (0 means all)"};
+  Configurable<std::vector<int>> pdgCodes{"pdgCodes", {211}, "List of signed PDG codes to consider for performance evaluation"};
   HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
   ConfigurableAxis ptAxis{"ptAxis", {500, 0, 100}, "#it{p}_{T} (GeV/#it{c})"};
   ConfigurableAxis etaAxis{"etaAxis", {100, -5.f, 5.f}, "#eta"};
   ConfigurableAxis invPtDeltaAxis{"invPtDeltaAxis", {1000, -0.1f, 0.1f}, "1./#it{p}_{T}^{gen} - 1./#it{p}_{T}^{reco} (GeV/#it{c})^{-1}"};
 
+  // Histograms
+  std::shared_ptr<TH1> particlePdgCodes;
+  std::map<int, std::shared_ptr<TH1>> particlePtDistribution;
+  std::map<int, std::shared_ptr<TH1>> particleEtaDistribution;
+  std::map<int, std::shared_ptr<TH1>> ptDistribution;
+  std::map<int, std::shared_ptr<TH2>> ptResolutionVsPt;
+  std::map<int, std::shared_ptr<TProfile2D>> ptResolutionVsEta;
+  std::map<int, std::shared_ptr<TH2>> invPtResolutionVsPt;
+  std::map<int, std::shared_ptr<TProfile2D>> invPtResolutionVsEta;
+  std::map<int, std::shared_ptr<TH2>> dcaXyResolutionVsPt;
+  std::map<int, std::shared_ptr<TH2>> dcaZResolutionVsPt;
+  std::map<int, std::shared_ptr<TH2>> covariancePtPtVsPt;
+  std::map<int, std::shared_ptr<TH2>> covarianceDcaXyDcaXyVsPt;
+  std::map<int, std::shared_ptr<TH2>> covarianceDcaZDcaZVsPt;
+
   void init(o2::framework::InitContext&)
   {
+    // Check if pdgCodes has duplicates and remove them
+    std::sort(pdgCodes.value.begin(), pdgCodes.value.end());
+    pdgCodes.value.erase(std::unique(pdgCodes.value.begin(), pdgCodes.value.end()), pdgCodes.value.end());
+    // Add 0 to the list of pdgCodes if it is not already present
+    if (std::find(pdgCodes.value.begin(), pdgCodes.value.end(), 0) == pdgCodes.value.end()) {
+      pdgCodes.value.push_back(0);
+    }
     const AxisSpec axisPtDelta{100, -1, 1, "(#it{p}_{T}^{reco} - #it{p}_{T}^{gen}) / #it{p}_{T}^{gen}"};
     const AxisSpec axisDcaXy{100, -1, 1, "DCA_{xy} (cm)"};
     const AxisSpec axisDcaZ{100, -1, 1, "DCA_{z} (cm)"};
     const AxisSpec axisCovariancePtPt{100, 0, 10, "cov(#it{p}_{T}, #it{p}_{T}) ((GeV/#it{c})^{2})"};
     const AxisSpec axisCovarianceDcaXyDcaXy{100, 0, 1, "cov(DCA_{xy}, DCA_{xy}) (cm^{2})"};
     const AxisSpec axisCovarianceDcaZDcaZ{100, 0, 1, "cov(DCA_{z}, DCA_{z}) (cm^{2})"};
-    particlePdgCodes = histos.add<TH1>("particlePdgCodes", "", kTH1D, {AxisSpec{100, -0.5, 99.5, "PDG Code"}});
+    const int nBinsPdgCodes = static_cast<int>(pdgCodes.value.size());
+    particlePdgCodes = histos.add<TH1>("particlePdgCodes", "", kTH1D, {AxisSpec{nBinsPdgCodes, -0.5, -0.5 + nBinsPdgCodes, "PDG Code"}});
     for (const int& pdg : pdgCodes.value) {
       std::string prefix = Form("%i", pdg);
       if (pdg < 0) {
@@ -80,6 +91,7 @@ struct Alice3TrackingPerformance {
       prefix += "/";
       auto histoName = [&](const std::string& name) { return Form("%s%s%s", prefix.c_str(), name.c_str(), tag.c_str()); };
       particlePtDistribution[pdg] = histos.add<TH1>(histoName("particlePtDistribution"), "", kTH1D, {ptAxis});
+      LOG(info) << "Added histogram: " << histoName("particlePtDistribution");
       particleEtaDistribution[pdg] = histos.add<TH1>(histoName("particleEtaDistribution"), "", kTH1D, {etaAxis});
 
       ptDistribution[pdg] = histos.add<TH1>(histoName("ptDistribution"), "", kTH1D, {ptAxis});
@@ -125,7 +137,7 @@ struct Alice3TrackingPerformance {
       }
       const auto& mcParticle = track.mcParticle();
       const float ptResolution = (track.pt() - mcParticle.pt()) / mcParticle.pt();
-      const float invptResolution = 1.f / track.pt() - 1.f / mcParticle.pt();
+      const float invptResolution = (track.pt() > 0.f ? 1.f / track.pt() : 0.f) - (mcParticle.pt() > 0.f ? 1.f / mcParticle.pt() : 0.f);
       const float covariancePtPt = track.sigma1Pt() * track.sigma1Pt() * track.pt() * track.pt() * track.pt() * track.pt();
 
       auto fillResolutionHistograms = [&](const int p) {
@@ -137,8 +149,8 @@ struct Alice3TrackingPerformance {
         dcaXyResolutionVsPt[p]->Fill(mcParticle.pt(), track.dcaXY());
         dcaZResolutionVsPt[p]->Fill(mcParticle.pt(), track.dcaZ());
         covariancePtPtVsPt[p]->Fill(mcParticle.pt(), covariancePtPt);
-        covarianceDcaXyDcaXyVsPt[p]->Fill(mcParticle.pt(), track.sigmaDcaXY2());
-        covarianceDcaZDcaZVsPt[p]->Fill(mcParticle.pt(), track.sigmaDcaZ2());
+        // covarianceDcaXyDcaXyVsPt[p]->Fill(mcParticle.pt(), track.sigmaDcaXY2());
+        // covarianceDcaZDcaZVsPt[p]->Fill(mcParticle.pt(), track.sigmaDcaZ2());
       };
 
       fillResolutionHistograms(0);