[PWGCF/FemtoUniverse] Add EfficiencyCorrection class & fix pT binning in FemtoUniverseParticleHisto

PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCalculator.h

@@ -19,6 +19,7 @@
 #include <vector>
 #include <map>
 #include <string>
+#include <algorithm>
 
 #include "Framework/Configurable.h"
 #include "CCDB/BasicCCDBManager.h"
@@ -29,11 +30,11 @@
 {
 enum ParticleNo : size_t {
   ONE = 1,
-  TWO
+  TWO,
 };
 
 template <size_t T>
 concept isOneOrTwo = T == ParticleNo::ONE || T == ParticleNo::TWO;
 
 template <typename T>
 consteval auto getHistDim() -> int
@@ -50,8 +51,8 @@
 
 struct EfficiencyConfigurableGroup : ConfigurableGroup {
   Configurable<bool> confEfficiencyApplyCorrections{"confEfficiencyApplyCorrections", false, "Should apply corrections from efficiency"};
-  Configurable<int> confEfficiencyCCDBTrainNumber{"confEfficiencyCCDBTrainNumber", -1, "Train number for which to query CCDB objects (set to -1 to ignore)"};
-  Configurable<std::vector<std::string>> confEfficiencyCCDBTimestamps{"confEfficiencyCCDBTimestamps", {}, "Timestamps of efficiency histograms in CCDB, to query for specific objects (default: ['-1', '-1'], gets the latest valid objects for both)"};
+  Configurable<int> confEfficiencyCCDBTrainNumber{"confEfficiencyCCDBTrainNumber", 0, "Train number for which to query CCDB objects (set 0 to ignore)"};
+  Configurable<std::vector<std::string>> confEfficiencyCCDBTimestamps{"confEfficiencyCCDBTimestamps", {}, "Timestamps of efficiency histograms in CCDB, to query for specific objects (default: [], set 0 to ignore, useful when running subwagons)"};
 
   // NOTE: in the future we might move the below configurables to a separate struct, eg. CCDBConfigurableGroup
   Configurable<std::string> confCCDBUrl{"confCCDBUrl", "http://alice-ccdb.cern.ch", "CCDB URL to be used"};
@@ -73,17 +74,23 @@
     ccdb.setLocalObjectValidityChecking();
     ccdb.setFatalWhenNull(false);
 
-    int64_t now = duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
+    auto now = duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
     ccdb.setCreatedNotAfter(now);
 
     shouldApplyCorrections = config->confEfficiencyApplyCorrections;
 
     if (config->confEfficiencyApplyCorrections && !config->confEfficiencyCCDBTimestamps.value.empty()) {
-      for (const auto& timestamp : config->confEfficiencyCCDBTimestamps.value) {
-        hLoaded.push_back(loadEfficiencyFromCCDB(std::stol(timestamp)));
+      for (auto idx = 0UL; idx < config->confEfficiencyCCDBTimestamps.value.size(); idx++) {
+        auto timestamp = 0L;
+        try {
+          timestamp = std::max(0L, std::stol(config->confEfficiencyCCDBTimestamps.value[idx]));
+        } catch (const std::exception&) {
+          LOGF(error, notify("Could not parse CCDB timestamp \"%s\""), config->confEfficiencyCCDBTimestamps.value[idx]);
+          continue;
+        }
+
+        hLoaded[idx] = timestamp > 0 ? loadEfficiencyFromCCDB(timestamp) : nullptr;
       }
-
-      LOGF(info, notify("Successfully loaded %d efficiency histogram(s)"), hLoaded.size());
     }
   }
 
@@ -92,15 +99,12 @@
   auto getWeight(ParticleNo partNo, const BinVars&... binVars) const -> float
   {
     auto weight = 1.0f;
+    auto hEff = hLoaded[partNo - 1];
 
-    if (partNo - 1 < config->confEfficiencyCCDBTimestamps.value.size()) {
-      auto hEff = hLoaded[partNo - 1];
-
-      if (shouldApplyCorrections && hEff) {
-        auto bin = hEff->FindBin(binVars...);
-        auto eff = hEff->GetBinContent(bin);
-        weight /= eff > 0 ? eff : 1.0f;
-      }
+    if (shouldApplyCorrections && hEff) {
+      auto bin = hEff->FindBin(binVars...);
+      auto eff = hEff->GetBinContent(bin);
+      weight /= eff > 0 ? eff : 1.0f;
     }
 
     return weight;
@@ -143,6 +147,7 @@
       LOGF(warn, notify("Histogram \"%s/%ld\" has been loaded, but it is empty"), config->confCCDBPath.value, timestamp);
     }
 
+    LOGF(info, notify("Successfully loaded %ld"), timestamp);
     return hEff;
   }
 
@@ -151,7 +156,7 @@
   bool shouldApplyCorrections = false;
 
   o2::ccdb::BasicCCDBManager& ccdb{o2::ccdb::BasicCCDBManager::instance()};
-  std::vector<HistType*> hLoaded{};
+  std::array<HistType*, 2> hLoaded{};
 };
 
 } // namespace o2::analysis::femto_universe::efficiency

PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h

@@ -0,0 +1,154 @@
+// Copyright 2019-2022 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 FemtoUniverseEfficiencyCorrection.h
+/// \brief Abstraction for applying efficiency corrections based on weights from CCDB
+/// \author Dawid Karpiński, WUT Warsaw, dawid.karpinski@cern.ch
+
+#ifndef PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSEEFFICIENCYCORRECTION_H_
+#define PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSEEFFICIENCYCORRECTION_H_
+
+#include <vector>
+#include <string>
+#include <algorithm>
+
+#include "Framework/Configurable.h"
+#include "CCDB/BasicCCDBManager.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+
+namespace o2::analysis::femto_universe::efficiency_correction
+{
+enum ParticleNo : size_t {
+  ONE = 1,
+  TWO,
+};
+
+template <size_t T>
+concept isOneOrTwo = T == ParticleNo::ONE || T == ParticleNo::TWO;
+
+template <typename T>
+consteval auto getHistDim() -> int
+{
+  if (std::is_same_v<T, TH1>)
+    return 1;
+  else if (std::is_same_v<T, TH2>)
+    return 2;
+  else if (std::is_same_v<T, TH3>)
+    return 3;
+  else
+    return -1;
+}
+
+struct EffCorConfigurableGroup : framework::ConfigurableGroup {
+  framework::Configurable<bool> confEffCorApply{"confEffCorApply", false, "[Efficiency Correction] Should apply efficiency corrections"};
+  framework::Configurable<std::string> confEffCorCCDBUrl{"confEffCorCCDBUrl", "http://alice-ccdb.cern.ch", "[Efficiency Correction] CCDB URL to use"};
+  framework::Configurable<std::string> confEffCorCCDBPath{"confEffCorCCDBPath", "", "[Efficiency Correction] CCDB path to histograms"};
+  framework::Configurable<std::vector<std::string>> confEffCorCCDBTimestamps{"confEffCorCCDBTimestamps", {}, "[Efficiency Correction] Timestamps of histograms in CCDB (0 can be used as a placeholder, e.g. when running subwagons)"};
+};
+
+template <typename HistType>
+  requires std::is_base_of_v<TH1, HistType>
+class EfficiencyCorrection
+{
+ public:
+  explicit EfficiencyCorrection(EffCorConfigurableGroup* config) : config(config) // o2-linter: disable=name/function-variable
+  {
+  }
+
+  auto init() -> void
+  {
+    ccdb.setURL(config->confEffCorCCDBUrl);
+    ccdb.setLocalObjectValidityChecking();
+    ccdb.setFatalWhenNull(false);
+
+    auto now = duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
+    ccdb.setCreatedNotAfter(now);
+
+    shouldApplyCorrection = config->confEffCorApply;
+
+    if (shouldApplyCorrection && !config->confEffCorCCDBTimestamps.value.empty()) {
+      for (auto idx = 0UL; idx < config->confEffCorCCDBTimestamps.value.size(); idx++) {
+        auto timestamp = 0L;
+        try {
+          timestamp = std::max(0L, std::stol(config->confEffCorCCDBTimestamps.value[idx]));
+        } catch (const std::exception&) {
+          LOGF(error, notify("Could not parse CCDB timestamp \"%s\""), config->confEffCorCCDBTimestamps.value[idx]);
+          continue;
+        }
+
+        hLoaded[idx] = timestamp > 0 ? loadHistFromCCDB(timestamp) : nullptr;
+      }
+    }
+  }
+
+  template <typename... BinVars>
+    requires(sizeof...(BinVars) == getHistDim<HistType>())
+  auto getWeight(ParticleNo partNo, const BinVars&... binVars) const -> float
+  {
+    auto weight = 1.0f;
+    auto hWeights = hLoaded[partNo - 1];
+
+    if (shouldApplyCorrection && hWeights) {
+      auto bin = hWeights->FindBin(binVars...);
+      weight = hWeights->GetBinContent(bin);
+    }
+
+    return weight;
+  }
+
+ private:
+  static inline auto notify(const std::string& msg) -> const std::string
+  {
+    return fmt::format("[EFFICIENCY CORRECTION] {}", msg);
+  }
+
+  static auto isHistEmpty(HistType* hist) -> bool
+  {
+    if (!hist) {
+      return true;
+    }
+    for (auto idx = 0; idx <= hist->GetNbinsX() + 1; idx++) {
+      if (hist->GetBinContent(idx) > 0) {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  auto loadHistFromCCDB(const int64_t timestamp) const -> HistType*
+  {
+    auto hWeights = ccdb.getForTimeStamp<HistType>(config->confEffCorCCDBPath, timestamp);
+    if (!hWeights || hWeights->IsZombie()) {
+      LOGF(error, notify("Could not load histogram \"%s/%ld\""), config->confEffCorCCDBPath.value, timestamp);
+      return nullptr;
+    }
+
+    if (isHistEmpty(hWeights)) {
+      LOGF(warn, notify("Histogram \"%s/%ld\" has been loaded, but it is empty"), config->confEffCorCCDBUrl.value, timestamp);
+    }
+
+    LOGF(info, notify("Successfully loaded %ld"), timestamp);
+    return hWeights;
+  }
+
+  EffCorConfigurableGroup* config{};
+
+  bool shouldApplyCorrection = false;
+
+  o2::ccdb::BasicCCDBManager& ccdb{o2::ccdb::BasicCCDBManager::instance()};
+  std::array<HistType*, 2> hLoaded{};
+};
+
+} // namespace o2::analysis::femto_universe::efficiency_correction
+
+#endif // PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSEEFFICIENCYCORRECTION_H_

PWGCF/FemtoUniverse/Core/FemtoUniverseParticleHisto.h

@@ -26,7 +26,7 @@
 #include "Framework/HistogramRegistry.h"
 #include "CommonConstants/MathConstants.h"
 
 using namespace o2::framework; // o2-linter: disable=using-directive
 
 namespace o2::analysis::femto_universe // o2-linter: disable=name/namespace
 {
@@ -51,11 +51,11 @@
   /// \param tempFitVarpTAxis axis object for the pT axis in the pT vs. tempFitVar plots
   /// \param tempFitVarAxis axis object for the tempFitVar axis
   template <o2::aod::femtouniverse_mc_particle::MCType mc, typename T>
   void init_base(std::string folderName, std::string tempFitVarAxisTitle, T& tempFitVarpTAxis, T& tempFitVarAxis) // o2-linter: disable=name/function-variable
   {
     std::string folderSuffix = static_cast<std::string>(o2::aod::femtouniverse_mc_particle::MCTypeName[mc]).c_str();
     /// Histograms of the kinematic properties
-    mHistogramRegistry->add((folderName + folderSuffix + "/hPt").c_str(), "; #it{p}_{T} (GeV/#it{c}); Entries", kTH1F, {{240, 0, 6}});
+    mHistogramRegistry->add((folderName + folderSuffix + "/hPt").c_str(), "; #it{p}_{T} (GeV/#it{c}); Entries", kTH1F, {tempFitVarpTAxis});
     mHistogramRegistry->add((folderName + folderSuffix + "/hEta").c_str(), "; #eta; Entries", kTH1F, {{200, -1.5, 1.5}});
     mHistogramRegistry->add((folderName + folderSuffix + "/hPhi").c_str(), "; #phi; Entries", kTH1F, {{200, 0, o2::constants::math::TwoPI}});
     mHistogramRegistry->add((folderName + folderSuffix + "/hPhiEta").c_str(), "; #phi; #eta", kTH2F, {{200, 0, o2::constants::math::TwoPI}, {200, -1.5, 1.5}});
@@ -68,7 +68,7 @@
 
   // comment
   template <o2::aod::femtouniverse_mc_particle::MCType mc>
   void init_debug(std::string folderName) // o2-linter: disable=name/function-variable
   {
     std::string folderSuffix = static_cast<std::string>(o2::aod::femtouniverse_mc_particle::MCTypeName[mc]).c_str();
     if constexpr (mParticleType == o2::aod::femtouniverseparticle::ParticleType::kTrack || mParticleType == o2::aod::femtouniverseparticle::ParticleType::kV0Child || mParticleType == o2::aod::femtouniverseparticle::ParticleType::kCascadeBachelor || mParticleType == o2::aod::femtouniverseparticle::ParticleType::kMCTruthTrack) {
@@ -135,7 +135,7 @@
     /// Particle-type specific histograms
     std::string folderSuffix = static_cast<std::string>(o2::aod::femtouniverse_mc_particle::MCTypeName[o2::aod::femtouniverse_mc_particle::MCType::kTruth]).c_str();
 
-    mHistogramRegistry->add((folderName + folderSuffix + "/hPt_ReconNoFake").c_str(), "; #it{p}_{T} (GeV/#it{c}); Entries", kTH1F, {{240, 0, 6}});
+    mHistogramRegistry->add((folderName + folderSuffix + "/hPt_ReconNoFake").c_str(), "; #it{p}_{T} (GeV/#it{c}); Entries", kTH1F, {tempFitVarpTAxis});
 
     if constexpr (mParticleType == o2::aod::femtouniverseparticle::ParticleType::kTrack || mParticleType == o2::aod::femtouniverseparticle::ParticleType::kV0Child || mParticleType == o2::aod::femtouniverseparticle::ParticleType::kCascadeBachelor || mParticleType == o2::aod::femtouniverseparticle::ParticleType::kMCTruthTrack) {
       /// Track histograms
@@ -375,17 +375,17 @@
         } else if (confPDG == mConfPDGCodePart[1]) {
           binPDG = 1;
         } else if (confPDG == mConfPDGCodePart[2]) {
           binPDG = 2; // o2-linter: disable=pdg/explicit-code
         } else {
           binPDG = 3; // o2-linter: disable=pdg/explicit-code
         }
         if (std::abs(pdgcode) == 211) {
           mHistogramRegistry->fill(histFolder + HIST("_MC/hMisidentification"),
                                    binPDG, 0, part.pt());
         } else if (std::abs(pdgcode) == 321) {
           mHistogramRegistry->fill(histFolder + HIST("_MC/hMisidentification"),
                                    binPDG, 1, part.pt());
         } else if (std::abs(pdgcode) == 2212) {
           mHistogramRegistry->fill(histFolder + HIST("_MC/hMisidentification"),
                                    binPDG, 2, part.pt());
         } else {

PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx

@@ -33,11 +33,11 @@
 #include "PWGCF/FemtoUniverse/Core/FemtoUniverseContainer.h"
 #include "PWGCF/FemtoUniverse/Core/FemtoUniverseDetaDphiStar.h"
 #include "PWGCF/FemtoUniverse/Core/FemtoUniverseTrackSelection.h"
-#include "PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCalculator.h"
+#include "PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h"
 
 using namespace o2;
 using namespace o2::analysis::femto_universe;
-using namespace o2::analysis::femto_universe::efficiency;
+using namespace o2::analysis::femto_universe::efficiency_correction;
 using namespace o2::framework;
 using namespace o2::framework::expressions;
 using namespace o2::soa;
@@ -177,8 +177,8 @@ struct FemtoUniversePairTaskTrackTrackExtended {
   HistogramRegistry resultRegistry{"Correlations", {}, OutputObjHandlingPolicy::AnalysisObject};
   HistogramRegistry mixQaRegistry{"mixQaRegistry", {}, OutputObjHandlingPolicy::AnalysisObject};
 
-  EfficiencyConfigurableGroup effConfGroup;
-  EfficiencyCalculator<TH1> efficiencyCalculator{&effConfGroup};
+  EffCorConfigurableGroup effCorConfGroup;
+  EfficiencyCorrection<TH1> effCorrection{&effCorConfGroup};
 
   /// @brief Counter for particle swapping
   int fNeventsProcessed = 0;
@@ -328,7 +328,7 @@ struct FemtoUniversePairTaskTrackTrackExtended {
         hMCTruth2.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarPDGBins, false, tracktwofilter.confPDGCodePartTwo, false);
       }
     }
-    efficiencyCalculator.init();
+    effCorrection.init();
 
     eventHisto.init(&qaRegistry);
     trackHistoPartOne.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, trackonefilter.confPDGCodePartOne, true); // last true = isDebug
@@ -505,9 +505,9 @@ struct FemtoUniversePairTaskTrackTrackExtended {
           continue;
         }
 
-        float weight = efficiencyCalculator.getWeight(ParticleNo::ONE, p1.pt());
+        float weight = effCorrection.getWeight(ParticleNo::ONE, p1.pt());
         if (!confIsSame) {
-          weight *= efficiencyCalculator.getWeight(ParticleNo::TWO, p2.pt());
+          weight *= effCorrection.getWeight(ParticleNo::TWO, p2.pt());
         }
 
         if (swpart)
@@ -566,9 +566,9 @@ struct FemtoUniversePairTaskTrackTrackExtended {
           continue;
         }
 
-        float weight = efficiencyCalculator.getWeight(ParticleNo::ONE, p1.pt());
+        float weight = effCorrection.getWeight(ParticleNo::ONE, p1.pt());
         if (!confIsSame) {
-          weight *= efficiencyCalculator.getWeight(ParticleNo::TWO, p2.pt());
+          weight *= effCorrection.getWeight(ParticleNo::TWO, p2.pt());
         }
 
         sameEventCont.setPair<isMC>(p1, p2, multCol, twotracksconfigs.confUse3D, weight);
@@ -674,9 +674,9 @@ struct FemtoUniversePairTaskTrackTrackExtended {
         }
       }
 
-      float weight = efficiencyCalculator.getWeight(ParticleNo::ONE, p1.pt());
+      float weight = effCorrection.getWeight(ParticleNo::ONE, p1.pt());
       if (!confIsSame) {
-        weight *= efficiencyCalculator.getWeight(ParticleNo::TWO, p2.pt());
+        weight *= effCorrection.getWeight(ParticleNo::TWO, p2.pt());
       }
 
       if (swpart)