From 1be3775d423be3addb55e268025bad4bf37ad67d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dawid=20Karpi=C5=84ski?= Date: Mon, 12 May 2025 19:02:18 +0200 Subject: [PATCH 1/9] fix: correct classification of material --- .../TableProducer/femtoUniverseProducerTask.cxx | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx index ce06e110b2c..1f532014b93 100644 --- a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx +++ b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx @@ -708,10 +708,11 @@ struct FemtoUniverseProducerTask { particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kPrimary; } else if (!motherparticlesMC.empty()) { auto motherparticleMC = motherparticlesMC.front(); - if (motherparticleMC.producedByGenerator()) + if (motherparticleMC.producedByGenerator()) { particleOrigin = checkDaughterType(fdparttype, motherparticleMC.pdgCode()); - } else { - particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kMaterial; + } else { + particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kMaterial; + } } } else { particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kFake; From cbd8db0acea137b5850a4efc47e4c2b1e58e189b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dawid=20Karpi=C5=84ski?= Date: Thu, 15 May 2025 21:21:41 +0200 Subject: [PATCH 2/9] feat: fill 3d histograms for efficiency correction --- .../Core/FemtoUniverseEfficiencyCorrection.h | 167 ++++++++++++++---- ...emtoUniversePairTaskTrackTrackExtended.cxx | 63 +++---- 2 files changed, 163 insertions(+), 67 deletions(-) diff --git a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h index b866649ab04..da751a24a52 100644 --- a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h +++ b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h @@ -16,15 +16,22 @@ #ifndef PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSEEFFICIENCYCORRECTION_H_ #define PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSEEFFICIENCYCORRECTION_H_ +#include +#include +#include +#include +#include + +#include +#include +#include +#include + #include #include #include -#include "Framework/Configurable.h" -#include "CCDB/BasicCCDBManager.h" -#include "TH1.h" -#include "TH2.h" -#include "TH3.h" +#include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h" namespace o2::analysis::femto_universe::efficiency_correction { @@ -34,30 +41,17 @@ enum ParticleNo : size_t { }; template -concept isOneOrTwo = T == ParticleNo::ONE || T == ParticleNo::TWO; - -template -consteval auto getHistDim() -> int -{ - if (std::is_same_v) - return 1; - else if (std::is_same_v) - return 2; - else if (std::is_same_v) - return 3; - else - return -1; -} +concept IsOneOrTwo = T == ParticleNo::ONE || T == ParticleNo::TWO; struct EffCorConfigurableGroup : framework::ConfigurableGroup { framework::Configurable confEffCorApply{"confEffCorApply", false, "[Efficiency Correction] Should apply efficiency corrections"}; + framework::Configurable confEffCorFillHist{"confEffCorFillHist", false, "[Efficiency Correction] Should fill histograms for efficiency corrections"}; framework::Configurable confEffCorCCDBUrl{"confEffCorCCDBUrl", "http://alice-ccdb.cern.ch", "[Efficiency Correction] CCDB URL to use"}; framework::Configurable confEffCorCCDBPath{"confEffCorCCDBPath", "", "[Efficiency Correction] CCDB path to histograms"}; framework::Configurable> confEffCorCCDBTimestamps{"confEffCorCCDBTimestamps", {}, "[Efficiency Correction] Timestamps of histograms in CCDB (0 can be used as a placeholder, e.g. when running subwagons)"}; + framework::Configurable> confEffCorVariables{"confEffCorVariables", {"pt"}, "[Efficiency Correction] Variables for efficiency correction histogram dimensions (available: pt, eta, cent-mult)"}; }; -template - requires std::is_base_of_v class EfficiencyCorrection { public: @@ -65,8 +59,25 @@ class EfficiencyCorrection { } - auto init() -> void + auto init(framework::HistogramRegistry* registry, std::vector axisSpecs) -> void { + shouldFillHistograms = config->confEffCorFillHist; + + histRegistry = registry; + if (shouldFillHistograms) { + for (const auto& suffix : histSuffix) { + auto path = std::format("{}/{}", histDirectory, suffix); + + registry->add((path + "/hMCTruth").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); + + registry->add((path + "/hPrimary").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); + registry->add((path + "/hSecondary").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); + registry->add((path + "/hMaterial").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); + registry->add((path + "/hFake").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); + registry->add((path + "/hOther").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); + } + } + ccdb.setURL(config->confEffCorCCDBUrl); ccdb.setLocalObjectValidityChecking(); ccdb.setFatalWhenNull(false); @@ -91,15 +102,104 @@ class EfficiencyCorrection } } - template - requires(sizeof...(BinVars) == getHistDim()) - auto getWeight(ParticleNo partNo, const BinVars&... binVars) const -> float + template + requires IsOneOrTwo + void fillTruthHist(auto particle) + { + if (!shouldFillHistograms) { + return; + } + + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hMCTruth"), + particle.pt(), + particle.eta(), + particle.fdCollision().multV0M()); + } + + template + requires IsOneOrTwo + void fillRecoHist(auto particle, int particlePDG) + { + if (!shouldFillHistograms) { + return; + } + + if (!particle.has_fdMCParticle()) { + return; + } + + auto mcParticle = particle.fdMCParticle(); + + if (mcParticle.pdgMCTruth() == particlePDG) { + // TODO question: fill with particle vs mcParticle, pt, eta, multV0M? + switch (mcParticle.partOriginMCTruth()) { + case (o2::aod::femtouniverse_mc_particle::kPrimary): + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hPrimary"), + mcParticle.pt(), + particle.eta(), + particle.fdCollision().multV0M()); + break; + + case (o2::aod::femtouniverse_mc_particle::kDaughter): + case (o2::aod::femtouniverse_mc_particle::kDaughterLambda): + case (o2::aod::femtouniverse_mc_particle::kDaughterSigmaplus): + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hSecondary"), + mcParticle.pt(), + particle.eta(), + particle.fdCollision().multV0M()); + break; + + case (o2::aod::femtouniverse_mc_particle::kMaterial): + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hMaterial"), + mcParticle.pt(), + particle.eta(), + particle.fdCollision().multV0M()); + break; + + case (o2::aod::femtouniverse_mc_particle::kFake): + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hFake"), + mcParticle.pt(), + particle.eta(), + particle.fdCollision().multV0M()); + break; + + default: + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hOther"), + mcParticle.pt(), + particle.eta(), + particle.fdCollision().multV0M()); + break; + } + } + } + + auto getWeight(ParticleNo partNo, auto particle) const -> float { auto weight = 1.0f; auto hWeights = hLoaded[partNo - 1]; if (shouldApplyCorrection && hWeights) { - auto bin = hWeights->FindBin(binVars...); + auto dim = hWeights->GetDimension(); + + if (dim != 3) { + LOGF(fatal, notify("Histogram \"%s\" has wrong dimension %d != 3"), config->confEffCorCCDBPath.value, dim); + return weight; + } + + auto bin = -1; + if (config->confEffCorVariables.value == std::vector{"pt"}) { + auto projected = hWeights->Project3D("x"); + bin = projected->FindBin(particle.pt()); + } else if (config->confEffCorVariables.value == std::vector{"pt", "eta"}) { + auto projected = hWeights->Project3D("xy"); + bin = projected->FindBin(particle.pt(), particle.eta()); + } else if (config->confEffCorVariables.value == std::vector{"pt", "cent-mult"}) { + auto projected = hWeights->Project3D("xz"); + bin = projected->FindBin(particle.pt(), particle.fdCollision().multV0M()); + } else { + LOGF(fatal, notify("unknown configuration for efficiency variables")); + return weight; + } weight = hWeights->GetBinContent(bin); } @@ -112,7 +212,7 @@ class EfficiencyCorrection return fmt::format("[EFFICIENCY CORRECTION] {}", msg); } - static auto isHistEmpty(HistType* hist) -> bool + static auto isHistEmpty(TH3* hist) -> bool { if (!hist) { return true; @@ -125,9 +225,9 @@ class EfficiencyCorrection return true; } - auto loadHistFromCCDB(const int64_t timestamp) const -> HistType* + auto loadHistFromCCDB(const int64_t timestamp) const -> TH3* { - auto hWeights = ccdb.getForTimeStamp(config->confEffCorCCDBPath, timestamp); + auto hWeights = ccdb.getForTimeStamp(config->confEffCorCCDBPath, timestamp); if (!hWeights || hWeights->IsZombie()) { LOGF(error, notify("Could not load histogram \"%s/%ld\""), config->confEffCorCCDBPath.value, timestamp); return nullptr; @@ -143,10 +243,15 @@ class EfficiencyCorrection EffCorConfigurableGroup* config{}; - bool shouldApplyCorrection = false; + bool shouldApplyCorrection{false}; + bool shouldFillHistograms{false}; o2::ccdb::BasicCCDBManager& ccdb{o2::ccdb::BasicCCDBManager::instance()}; - std::array hLoaded{}; + std::array hLoaded{}; + + framework::HistogramRegistry* histRegistry{}; + static constexpr std::string_view histDirectory{"EfficiencyCorrection"}; + static constexpr std::string_view histSuffix[2]{"one", "two"}; }; } // namespace o2::analysis::femto_universe::efficiency_correction diff --git a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx index 27e865ce552..c6aed0637f2 100644 --- a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx +++ b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx @@ -105,7 +105,6 @@ struct FemtoUniversePairTaskTrackTrackExtended { /// Histogramming for particle 1 FemtoUniverseParticleHisto trackHistoPartOne; - FemtoUniverseParticleHisto hMCTruth1; /// Particle 2 Configurable confIsSame{"confIsSame", false, "Pairs of the same particle"}; @@ -129,7 +128,6 @@ struct FemtoUniversePairTaskTrackTrackExtended { /// Histogramming for particle 2 FemtoUniverseParticleHisto trackHistoPartTwo; - FemtoUniverseParticleHisto hMCTruth2; /// Histogramming for Event FemtoUniverseEventHisto eventHisto; @@ -176,9 +174,10 @@ struct FemtoUniversePairTaskTrackTrackExtended { HistogramRegistry qaRegistry{"TrackQA", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry resultRegistry{"Correlations", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry mixQaRegistry{"mixQaRegistry", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry effCorrRegistry{"EfficiencyCorrection", {}, OutputObjHandlingPolicy::AnalysisObject}; EffCorConfigurableGroup effCorConfGroup; - EfficiencyCorrection effCorrection{&effCorConfGroup}; + EfficiencyCorrection effCorrection{&effCorConfGroup}; /// @brief Counter for particle swapping int fNeventsProcessed = 0; @@ -322,18 +321,13 @@ struct FemtoUniversePairTaskTrackTrackExtended { void init(InitContext&) { - if (twotracksconfigs.confIsMC) { - hMCTruth1.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarPDGBins, false, trackonefilter.confPDGCodePartOne, false); - if (!confIsSame) { - hMCTruth2.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarPDGBins, false, tracktwofilter.confPDGCodePartTwo, false); - } - } - effCorrection.init(); + auto axis = std::vector{{240, 0, 6}, {29, -2, 2}, {2000, 0, 20000}}; + effCorrection.init(&effCorrRegistry, axis); eventHisto.init(&qaRegistry); - trackHistoPartOne.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, trackonefilter.confPDGCodePartOne, true); // last true = isDebug + trackHistoPartOne.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, trackonefilter.confPDGCodePartOne, true, std::nullopt); // last true = isDebug if (!confIsSame) { - trackHistoPartTwo.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, tracktwofilter.confPDGCodePartTwo, true); // last true = isDebug + trackHistoPartTwo.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, tracktwofilter.confPDGCodePartTwo, true, std::nullopt); // last true = isDebug } mixQaRegistry.add("MixingQA/hSECollisionBins", ";bin;Entries", kTH1F, {{120, -0.5, 119.5}}); @@ -361,23 +355,12 @@ struct FemtoUniversePairTaskTrackTrackExtended { } template - requires isOneOrTwo - auto doMCTruth(FemtoUniverseParticleHisto hist, auto parts) -> void + requires IsOneOrTwo + auto doMCTruth(auto parts, int partPDG, int partCharge) -> void { - auto expectedPDG = 0; - auto expectedCharge = 0.0l; - - if constexpr (N == ParticleNo::ONE) { - expectedPDG = trackonefilter.confPDGCodePartOne; - expectedCharge = trackonefilter.confChargePart1; - } else if constexpr (N == ParticleNo::TWO) { - expectedPDG = tracktwofilter.confPDGCodePartTwo; - expectedCharge = tracktwofilter.confChargePart2; - } - for (const auto& particle : parts) { auto pdgCode = static_cast(particle.pidCut()); - if (pdgCode != expectedPDG) { + if (pdgCode != partPDG) { continue; } @@ -386,9 +369,11 @@ struct FemtoUniversePairTaskTrackTrackExtended { continue; } - if (sign(pdgParticle->Charge()) == sign(expectedCharge)) { - hist.template fillQA(particle); + if (sign(pdgParticle->Charge()) != sign(partCharge)) { + continue; } + + effCorrection.fillTruthHist(particle); } } @@ -432,6 +417,9 @@ struct FemtoUniversePairTaskTrackTrackExtended { } trackHistoPartOne.fillQA(part); + if constexpr (isMC) { + effCorrection.fillRecoHist(part, trackonefilter.confPDGCodePartOne); + } } if (!confIsSame) { @@ -456,6 +444,9 @@ struct FemtoUniversePairTaskTrackTrackExtended { } } trackHistoPartTwo.fillQA(part); + if constexpr (isMC) { + effCorrection.fillRecoHist(part, tracktwofilter.confPDGCodePartTwo); + } } /// Now build the combinations for non-identical particle pairs @@ -505,9 +496,9 @@ struct FemtoUniversePairTaskTrackTrackExtended { continue; } - float weight = effCorrection.getWeight(ParticleNo::ONE, p1.pt()); + float weight = effCorrection.getWeight(ParticleNo::ONE, p1); if (!confIsSame) { - weight *= effCorrection.getWeight(ParticleNo::TWO, p2.pt()); + weight *= effCorrection.getWeight(ParticleNo::TWO, p2); } if (swpart) @@ -566,9 +557,9 @@ struct FemtoUniversePairTaskTrackTrackExtended { continue; } - float weight = effCorrection.getWeight(ParticleNo::ONE, p1.pt()); + float weight = effCorrection.getWeight(ParticleNo::ONE, p1); if (!confIsSame) { - weight *= effCorrection.getWeight(ParticleNo::TWO, p2.pt()); + weight *= effCorrection.getWeight(ParticleNo::TWO, p2); } sameEventCont.setPair(p1, p2, multCol, twotracksconfigs.confUse3D, weight); @@ -602,11 +593,11 @@ struct FemtoUniversePairTaskTrackTrackExtended { fillCollision(col); auto groupMCTruth1 = partsOneMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); - doMCTruth<1>(hMCTruth1, groupMCTruth1); + doMCTruth(groupMCTruth1, trackonefilter.confPDGCodePartOne, trackonefilter.confChargePart1); if (!confIsSame) { auto groupMCTruth2 = partsTwoMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); - doMCTruth<2>(hMCTruth2, groupMCTruth2); + doMCTruth(groupMCTruth2, tracktwofilter.confPDGCodePartTwo, tracktwofilter.confChargePart2); } auto groupMCReco1 = partsOneMCReco->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); @@ -674,9 +665,9 @@ struct FemtoUniversePairTaskTrackTrackExtended { } } - float weight = effCorrection.getWeight(ParticleNo::ONE, p1.pt()); + float weight = effCorrection.getWeight(ParticleNo::ONE, p1); if (!confIsSame) { - weight *= effCorrection.getWeight(ParticleNo::TWO, p2.pt()); + weight *= effCorrection.getWeight(ParticleNo::TWO, p2); } if (swpart) From 7a90125d908c8c22c063e079a64cd12ed40a1e00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dawid=20Karpi=C5=84ski?= Date: Sun, 18 May 2025 21:47:13 +0200 Subject: [PATCH 3/9] feat: macro with projection and working efficiency correction calculation --- .../Core/FemtoUniverseEfficiencyCorrection.h | 49 ++++--- .../Macros/calculateEfficiency.cxx | 138 ++++++++++++++++++ .../femto_universe_efficiency_calculator.py | 0 ...emto_universe_efficiency_phi_calculator.py | 0 .../femtoUniverseProducerTask.cxx | 4 +- 5 files changed, 167 insertions(+), 24 deletions(-) create mode 100644 PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx rename PWGCF/FemtoUniverse/{Scripts => Macros}/femto_universe_efficiency_calculator.py (100%) rename PWGCF/FemtoUniverse/{Scripts => Macros}/femto_universe_efficiency_phi_calculator.py (100%) diff --git a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h index da751a24a52..59cbb5b81c2 100644 --- a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h +++ b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h @@ -25,7 +25,6 @@ #include #include #include -#include #include #include @@ -97,7 +96,17 @@ class EfficiencyCorrection continue; } - hLoaded[idx] = timestamp > 0 ? loadHistFromCCDB(timestamp) : nullptr; + if (timestamp > 0) { + if (config->confEffCorVariables->size() == 1) { + hLoaded[idx] = loadHistFromCCDB(timestamp); + } else if (config->confEffCorVariables->size() == 2) { + hLoaded[idx] = loadHistFromCCDB(timestamp); + } else { + LOGF(fatal, notify("unknown configuration for efficiency variables")); + } + } else { + hLoaded[idx] = nullptr; + } } } } @@ -131,12 +140,11 @@ class EfficiencyCorrection auto mcParticle = particle.fdMCParticle(); if (mcParticle.pdgMCTruth() == particlePDG) { - // TODO question: fill with particle vs mcParticle, pt, eta, multV0M? switch (mcParticle.partOriginMCTruth()) { case (o2::aod::femtouniverse_mc_particle::kPrimary): histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hPrimary"), mcParticle.pt(), - particle.eta(), + mcParticle.eta(), particle.fdCollision().multV0M()); break; @@ -145,28 +153,28 @@ class EfficiencyCorrection case (o2::aod::femtouniverse_mc_particle::kDaughterSigmaplus): histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hSecondary"), mcParticle.pt(), - particle.eta(), + mcParticle.eta(), particle.fdCollision().multV0M()); break; case (o2::aod::femtouniverse_mc_particle::kMaterial): histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hMaterial"), mcParticle.pt(), - particle.eta(), + mcParticle.eta(), particle.fdCollision().multV0M()); break; case (o2::aod::femtouniverse_mc_particle::kFake): histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hFake"), mcParticle.pt(), - particle.eta(), + mcParticle.eta(), particle.fdCollision().multV0M()); break; default: histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hOther"), mcParticle.pt(), - particle.eta(), + mcParticle.eta(), particle.fdCollision().multV0M()); break; } @@ -179,23 +187,19 @@ class EfficiencyCorrection auto hWeights = hLoaded[partNo - 1]; if (shouldApplyCorrection && hWeights) { - auto dim = hWeights->GetDimension(); - - if (dim != 3) { - LOGF(fatal, notify("Histogram \"%s\" has wrong dimension %d != 3"), config->confEffCorCCDBPath.value, dim); + auto dim = static_cast(hWeights->GetDimension()); + if (dim != config->confEffCorVariables.value.size()) { + LOGF(fatal, notify("Histogram \"%s\" has wrong dimension %d != %d"), config->confEffCorCCDBPath.value, dim, config->confEffCorVariables.value.size()); return weight; } auto bin = -1; if (config->confEffCorVariables.value == std::vector{"pt"}) { - auto projected = hWeights->Project3D("x"); - bin = projected->FindBin(particle.pt()); + bin = hLoaded[partNo - 1]->FindBin(particle.pt()); } else if (config->confEffCorVariables.value == std::vector{"pt", "eta"}) { - auto projected = hWeights->Project3D("xy"); - bin = projected->FindBin(particle.pt(), particle.eta()); + bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.eta()); } else if (config->confEffCorVariables.value == std::vector{"pt", "cent-mult"}) { - auto projected = hWeights->Project3D("xz"); - bin = projected->FindBin(particle.pt(), particle.fdCollision().multV0M()); + bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.fdCollision().multV0M()); } else { LOGF(fatal, notify("unknown configuration for efficiency variables")); return weight; @@ -212,7 +216,7 @@ class EfficiencyCorrection return fmt::format("[EFFICIENCY CORRECTION] {}", msg); } - static auto isHistEmpty(TH3* hist) -> bool + static auto isHistEmpty(TH1* hist) -> bool { if (!hist) { return true; @@ -225,9 +229,10 @@ class EfficiencyCorrection return true; } - auto loadHistFromCCDB(const int64_t timestamp) const -> TH3* + template + auto loadHistFromCCDB(const int64_t timestamp) const -> H* { - auto hWeights = ccdb.getForTimeStamp(config->confEffCorCCDBPath, timestamp); + auto hWeights = ccdb.getForTimeStamp(config->confEffCorCCDBPath, timestamp); if (!hWeights || hWeights->IsZombie()) { LOGF(error, notify("Could not load histogram \"%s/%ld\""), config->confEffCorCCDBPath.value, timestamp); return nullptr; @@ -247,7 +252,7 @@ class EfficiencyCorrection bool shouldFillHistograms{false}; o2::ccdb::BasicCCDBManager& ccdb{o2::ccdb::BasicCCDBManager::instance()}; - std::array hLoaded{}; + std::array hLoaded{}; framework::HistogramRegistry* histRegistry{}; static constexpr std::string_view histDirectory{"EfficiencyCorrection"}; diff --git a/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx b/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx new file mode 100644 index 00000000000..2f19bb2de81 --- /dev/null +++ b/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx @@ -0,0 +1,138 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fs = std::filesystem; + +auto* getHistogram(TFile* file, const std::string& name, const std::string& projection) +{ + // TODO: I think axes after projection are transposed + return dynamic_cast(file->Get(name.c_str()))->Project3D(projection.c_str()); +} + +template +auto* cloneHistogram(H* hist, const std::string& name) +{ + return dynamic_cast(hist->Clone(name.c_str())); +} + +void forEachBin(TH1* hist, auto func) +{ + if (hist->GetDimension() == 1) { + for (auto x{1}; x <= hist->GetNbinsX(); ++x) { + func(x, 0); + } + } else if (hist->GetDimension() == 2) { + for (auto x{1}; x <= hist->GetNbinsX(); ++x) { + for (auto y{1}; y <= hist->GetNbinsY(); ++y) { + func(x, y); + } + } + } else { + assert(false && "should not happen"); + } +} + +void calculateEfficiency(const fs::path& resultsPath, const fs::path& histPath, const std::string& projection) +{ + assert(projection == "x" || projection == "xy" || projection == "xz"); + + auto isAlien{false}; + if (resultsPath.string().starts_with("alien://")) { + TGrid::Connect("alien://"); + isAlien = true; + } + + auto* resultFile{TFile::Open(resultsPath.c_str())}; + assert(resultFile != nullptr && !resultFile->IsZombie()); + + using namespace std::chrono; + auto now{duration_cast(system_clock::now().time_since_epoch()).count()}; + auto outputPath{isAlien ? std::filesystem::current_path() : resultsPath.parent_path()}; + outputPath /= std::format("{}-effcor-{}.root", resultsPath.stem().string(), now); + + auto* outputFile{TFile::Open(outputPath.c_str(), "RECREATE")}; + assert(outputFile != nullptr && !outputFile->IsZombie()); + + auto* histTruth{getHistogram(resultFile, histPath / "hMCTruth", projection)}; + assert(histTruth); + + auto* histPrimary{getHistogram(resultFile, histPath / "hPrimary", projection)}; + assert(histPrimary); + + auto* histSecondary{getHistogram(resultFile, histPath / "hSecondary", projection)}; + assert(histSecondary); + + auto* histTotal{cloneHistogram(histPrimary, "hTotal")}; + histTotal->Add(histSecondary); + + auto* histEfficiency{cloneHistogram(histPrimary, "hEfficiency")}; + histEfficiency->Reset(); + + auto* histWeights{cloneHistogram(histPrimary, "hWeights")}; + histWeights->Reset(); + + forEachBin(histPrimary, [&](int x, int y) { + auto primVal{histPrimary->GetBinContent(x, y)}; + auto primErr{histPrimary->GetBinError(x, y)}; + + auto secVal{histSecondary->GetBinContent(x, y)}; + auto secErr{histSecondary->GetBinError(x, y)}; + + auto truthVal{histTruth->GetBinContent(x, y)}; + auto truthErr{histTruth->GetBinError(x, y)}; + + auto effVal{0.}; + auto effErr{0.}; + if (truthVal > 0) { + effVal = primVal / truthVal; + effErr = std::sqrt(std::pow(primErr / truthVal, 2) + std::pow((primVal * truthErr / std::pow(truthVal, 2)), 2)); + } + + histEfficiency->SetBinContent(x, y, effVal); + histEfficiency->SetBinError(x, y, effErr); + + auto totalVal{primVal + secVal}; + auto totalErr{std::hypot(primErr, secErr)}; + + auto contVal{0.}; + auto contErr{0.}; + if (totalVal > 0) { + contVal = secVal / totalVal; + contErr = std::sqrt(std::pow(secErr / totalVal, 2) + std::pow((secVal * totalErr / std::pow(totalVal, 2)), 2)); + } + + auto weightVal{0.}; + auto weightErr{0.}; + if (effVal > 0) { + weightVal = (1 - contVal) / effVal; + weightErr = std::sqrt(std::pow(contErr / effVal, 2) + std::pow((1 - contVal) * effErr / std::pow(effVal, 2), 2)); + } + + histWeights->SetBinContent(x, y, weightVal); + histWeights->SetBinError(x, y, weightErr); + }); + + outputFile->WriteTObject(histEfficiency); + outputFile->WriteTObject(histWeights); + + outputFile->Close(); + resultFile->Close(); +} + +int main(int argc, char** argv) +{ + assert(argc == 4); + calculateEfficiency(argv[1], argv[2], argv[3]); + return 0; +} diff --git a/PWGCF/FemtoUniverse/Scripts/femto_universe_efficiency_calculator.py b/PWGCF/FemtoUniverse/Macros/femto_universe_efficiency_calculator.py similarity index 100% rename from PWGCF/FemtoUniverse/Scripts/femto_universe_efficiency_calculator.py rename to PWGCF/FemtoUniverse/Macros/femto_universe_efficiency_calculator.py diff --git a/PWGCF/FemtoUniverse/Scripts/femto_universe_efficiency_phi_calculator.py b/PWGCF/FemtoUniverse/Macros/femto_universe_efficiency_phi_calculator.py similarity index 100% rename from PWGCF/FemtoUniverse/Scripts/femto_universe_efficiency_phi_calculator.py rename to PWGCF/FemtoUniverse/Macros/femto_universe_efficiency_phi_calculator.py diff --git a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx index 1f532014b93..55a19bf897c 100644 --- a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx +++ b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx @@ -903,8 +903,8 @@ struct FemtoUniverseProducerTask { { for (const auto& c : col) { const auto vtxZ = c.posZ(); - float mult = 0; - int multNtr = 0; + float mult = confIsRun3 ? c.multFV0M() : 0.5 * (c.multFV0M()); + int multNtr = confIsRun3 ? c.multNTracksPV() : c.multTracklets(); if (std::abs(vtxZ) > confEvtZvtx) { continue; From 0c9edb877e5b1e843c05931024ef565dd897f05c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dawid=20Karpi=C5=84ski?= Date: Mon, 19 May 2025 13:04:42 +0200 Subject: [PATCH 4/9] fix: project the other way around --- PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx b/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx index 2f19bb2de81..9883314c012 100644 --- a/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx +++ b/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx @@ -21,7 +21,7 @@ auto* getHistogram(TFile* file, const std::string& name, const std::string& proj } template -auto* cloneHistogram(H* hist, const std::string& name) +H* cloneHistogram(H* hist, const std::string& name) { return dynamic_cast(hist->Clone(name.c_str())); } @@ -45,7 +45,7 @@ void forEachBin(TH1* hist, auto func) void calculateEfficiency(const fs::path& resultsPath, const fs::path& histPath, const std::string& projection) { - assert(projection == "x" || projection == "xy" || projection == "xz"); + assert(projection == "x" || projection == "yx" || projection == "zx"); auto isAlien{false}; if (resultsPath.string().starts_with("alien://")) { From 8b0d8244f5d5eab3c404ce4666ff03e8ea059361 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dawid=20Karpi=C5=84ski?= Date: Mon, 19 May 2025 13:50:06 +0200 Subject: [PATCH 5/9] feat: add 3d hist to loading from ccdb --- .../Core/FemtoUniverseEfficiencyCorrection.h | 4 ++ .../Macros/calculateEfficiency.cxx | 43 ++++++++++++------- 2 files changed, 31 insertions(+), 16 deletions(-) diff --git a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h index 59cbb5b81c2..fbb7113abc2 100644 --- a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h +++ b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h @@ -101,6 +101,8 @@ class EfficiencyCorrection hLoaded[idx] = loadHistFromCCDB(timestamp); } else if (config->confEffCorVariables->size() == 2) { hLoaded[idx] = loadHistFromCCDB(timestamp); + } else if (config->confEffCorVariables->size() == 3) { + hLoaded[idx] = loadHistFromCCDB(timestamp); } else { LOGF(fatal, notify("unknown configuration for efficiency variables")); } @@ -200,6 +202,8 @@ class EfficiencyCorrection bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.eta()); } else if (config->confEffCorVariables.value == std::vector{"pt", "cent-mult"}) { bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.fdCollision().multV0M()); + } else if (config->confEffCorVariables.value == std::vector{"pt", "eta", "cent-mult"}) { + bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.eta(), particle.fdCollision().multV0M()); } else { LOGF(fatal, notify("unknown configuration for efficiency variables")); return weight; diff --git a/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx b/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx index 9883314c012..80da583e4c0 100644 --- a/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx +++ b/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx @@ -16,8 +16,11 @@ namespace fs = std::filesystem; auto* getHistogram(TFile* file, const std::string& name, const std::string& projection) { - // TODO: I think axes after projection are transposed - return dynamic_cast(file->Get(name.c_str()))->Project3D(projection.c_str()); + auto* hist{dynamic_cast(file->Get(name.c_str()))}; + if (projection == "none") { + return hist; + } + return dynamic_cast(hist->Project3D(projection.c_str())); } template @@ -30,12 +33,20 @@ void forEachBin(TH1* hist, auto func) { if (hist->GetDimension() == 1) { for (auto x{1}; x <= hist->GetNbinsX(); ++x) { - func(x, 0); + func(x, 0, 0); } } else if (hist->GetDimension() == 2) { for (auto x{1}; x <= hist->GetNbinsX(); ++x) { for (auto y{1}; y <= hist->GetNbinsY(); ++y) { - func(x, y); + func(x, y, 0); + } + } + } else if (hist->GetDimension() == 3) { + for (auto x{1}; x <= hist->GetNbinsX(); ++x) { + for (auto y{1}; y <= hist->GetNbinsY(); ++y) { + for (auto z{1}; z <= hist->GetNbinsZ(); ++z) { + func(x, y, z); + } } } } else { @@ -45,7 +56,7 @@ void forEachBin(TH1* hist, auto func) void calculateEfficiency(const fs::path& resultsPath, const fs::path& histPath, const std::string& projection) { - assert(projection == "x" || projection == "yx" || projection == "zx"); + assert(projection == "x" || projection == "yx" || projection == "zx" || projection == "none"); auto isAlien{false}; if (resultsPath.string().starts_with("alien://")) { @@ -82,15 +93,15 @@ void calculateEfficiency(const fs::path& resultsPath, const fs::path& histPath, auto* histWeights{cloneHistogram(histPrimary, "hWeights")}; histWeights->Reset(); - forEachBin(histPrimary, [&](int x, int y) { - auto primVal{histPrimary->GetBinContent(x, y)}; - auto primErr{histPrimary->GetBinError(x, y)}; + forEachBin(histPrimary, [&](int x, int y, int z) { + auto primVal{histPrimary->GetBinContent(x, y, z)}; + auto primErr{histPrimary->GetBinError(x, y, z)}; - auto secVal{histSecondary->GetBinContent(x, y)}; - auto secErr{histSecondary->GetBinError(x, y)}; + auto secVal{histSecondary->GetBinContent(x, y, z)}; + auto secErr{histSecondary->GetBinError(x, y, z)}; - auto truthVal{histTruth->GetBinContent(x, y)}; - auto truthErr{histTruth->GetBinError(x, y)}; + auto truthVal{histTruth->GetBinContent(x, y, z)}; + auto truthErr{histTruth->GetBinError(x, y, z)}; auto effVal{0.}; auto effErr{0.}; @@ -99,8 +110,8 @@ void calculateEfficiency(const fs::path& resultsPath, const fs::path& histPath, effErr = std::sqrt(std::pow(primErr / truthVal, 2) + std::pow((primVal * truthErr / std::pow(truthVal, 2)), 2)); } - histEfficiency->SetBinContent(x, y, effVal); - histEfficiency->SetBinError(x, y, effErr); + histEfficiency->SetBinContent(x, y, z, effVal); + histEfficiency->SetBinError(x, y, z, effErr); auto totalVal{primVal + secVal}; auto totalErr{std::hypot(primErr, secErr)}; @@ -119,8 +130,8 @@ void calculateEfficiency(const fs::path& resultsPath, const fs::path& histPath, weightErr = std::sqrt(std::pow(contErr / effVal, 2) + std::pow((1 - contVal) * effErr / std::pow(effVal, 2), 2)); } - histWeights->SetBinContent(x, y, weightVal); - histWeights->SetBinError(x, y, weightErr); + histWeights->SetBinContent(x, y, z, weightVal); + histWeights->SetBinError(x, y, z, weightErr); }); outputFile->WriteTObject(histEfficiency); From 70677726f58b2408dc64ae86699d591b54300bd9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dawid=20Karpi=C5=84ski?= Date: Mon, 19 May 2025 22:03:36 +0200 Subject: [PATCH 6/9] feat: use configurables for axis specifications --- .../femtoUniversePairTaskTrackTrackExtended.cxx | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx index c6aed0637f2..bb17fc05d2c 100644 --- a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx +++ b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx @@ -321,9 +321,6 @@ struct FemtoUniversePairTaskTrackTrackExtended { void init(InitContext&) { - auto axis = std::vector{{240, 0, 6}, {29, -2, 2}, {2000, 0, 20000}}; - effCorrection.init(&effCorrRegistry, axis); - eventHisto.init(&qaRegistry); trackHistoPartOne.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, trackonefilter.confPDGCodePartOne, true, std::nullopt); // last true = isDebug if (!confIsSame) { @@ -345,6 +342,15 @@ struct FemtoUniversePairTaskTrackTrackExtended { vPIDPartOne = trackonefilter.confPIDPartOne.value; vPIDPartTwo = tracktwofilter.confPIDPartTwo.value; kNsigma = twotracksconfigs.confTrkPIDnSigmaMax.value; + + effCorrection.init( + &effCorrRegistry, + { + static_cast(confTempFitVarpTBins), + {confEtaBins, -2, 2}, + confMultBins, + } // + ); } template From 576087799dd193eefcd0a93382668ebc97d1adfd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dawid=20Karpi=C5=84ski?= Date: Tue, 20 May 2025 19:22:29 +0200 Subject: [PATCH 7/9] feat: add flags to macro, add axes titles --- .../Core/FemtoUniverseEfficiencyCorrection.h | 12 +- .../FemtoUniverse/Macros/calculateEffCor.cxx | 239 ++++++++++++++++++ .../Macros/calculateEfficiency.cxx | 149 ----------- 3 files changed, 245 insertions(+), 155 deletions(-) create mode 100644 PWGCF/FemtoUniverse/Macros/calculateEffCor.cxx delete mode 100644 PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx diff --git a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h index fbb7113abc2..b8d8cf531bf 100644 --- a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h +++ b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h @@ -67,13 +67,13 @@ class EfficiencyCorrection for (const auto& suffix : histSuffix) { auto path = std::format("{}/{}", histDirectory, suffix); - registry->add((path + "/hMCTruth").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); + registry->add((path + "/hMCTruth").c_str(), "MCTruth; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); - registry->add((path + "/hPrimary").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); - registry->add((path + "/hSecondary").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); - registry->add((path + "/hMaterial").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); - registry->add((path + "/hFake").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); - registry->add((path + "/hOther").c_str(), "; ; Entries", framework::kTH3F, axisSpecs); + registry->add((path + "/hPrimary").c_str(), "Primary; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hSecondary").c_str(), "Secondary; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hMaterial").c_str(), "Material; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hFake").c_str(), "Fake; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hOther").c_str(), "Other; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); } } diff --git a/PWGCF/FemtoUniverse/Macros/calculateEffCor.cxx b/PWGCF/FemtoUniverse/Macros/calculateEffCor.cxx new file mode 100644 index 00000000000..a5c7625881e --- /dev/null +++ b/PWGCF/FemtoUniverse/Macros/calculateEffCor.cxx @@ -0,0 +1,239 @@ +// Copyright 2019-2020 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 calculateEffCor.cxx +/// \brief Macro for calculating efficiency corrections based on 3D histograms +/// \author Dawid Karpiński, WUT Warsaw, dawid.karpinski@cern.ch + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fs = std::filesystem; + +auto* getHistogram(TFile* file, const std::string& name) +{ + return dynamic_cast(file->Get(name.c_str())); +} + +auto* projectHistogram(TH3* hist, const std::string& projection) +{ + return hist->Project3D(projection.c_str()); +} + +template +H* cloneHistogram(H* hist, const std::string& name) +{ + return dynamic_cast(hist->Clone(name.c_str())); +} + +auto forEachBin(TH1* hist, auto func) -> void +{ + if (hist->GetDimension() == 1) { + for (auto x{1}; x <= hist->GetNbinsX(); ++x) { + func(x, 0, 0); + } + } else if (hist->GetDimension() == 2) { + for (auto x{1}; x <= hist->GetNbinsX(); ++x) { + for (auto y{1}; y <= hist->GetNbinsY(); ++y) { + func(x, y, 0); + } + } + } else if (hist->GetDimension() == 3) { + for (auto x{1}; x <= hist->GetNbinsX(); ++x) { + for (auto y{1}; y <= hist->GetNbinsY(); ++y) { + for (auto z{1}; z <= hist->GetNbinsZ(); ++z) { + func(x, y, z); + } + } + } + } else { + assert(false && "should not happen"); + } +} + +auto setAxisTitles(TH1* hist, const std::string& projection) -> void +{ + auto* xAxis{hist->GetXaxis()}; + auto* yAxis{hist->GetYaxis()}; + auto* zAxis{hist->GetZaxis()}; + + xAxis->SetTitle("#it{p}_{T} (GeV/#it{c})"); + + if (hist->GetDimension() == 2) { + if (projection == "yx") { + yAxis->SetTitle("#it{#eta} (rad)"); + } else if (projection == "zx") { + yAxis->SetTitle("Cent/Mult"); + } + } else if (hist->GetDimension() == 3) { + yAxis->SetTitle("#it{#eta} (rad)"); + zAxis->SetTitle("Cent/Mult"); + } +} + +auto calculateEfficiency(const fs::path& resultsPath, const fs::path& histPath, const std::string& projection) -> void +{ + assert(!resultsPath.empty() && !histPath.empty()); + if (projection != "" && projection != "x" && projection != "yx" && projection != "zx") { + std::cerr << "Error: projection must be one of: x, yx, zx\n"; + std::exit(1); + return; + } + + auto isAlien{false}; + if (resultsPath.string().starts_with("alien://")) { + TGrid::Connect("alien://"); + isAlien = true; + } + + auto* resultFile{TFile::Open(resultsPath.c_str())}; + assert(resultFile != nullptr && !resultFile->IsZombie()); + + using namespace std::chrono; + auto now{duration_cast(system_clock::now().time_since_epoch()).count()}; + auto outputPath{isAlien ? std::filesystem::current_path() : resultsPath.parent_path()}; + outputPath /= std::format("{}-effcor-{}.root", resultsPath.stem().string(), now); + + auto* outputFile{TFile::Open(outputPath.c_str(), "RECREATE")}; + assert(outputFile != nullptr && !outputFile->IsZombie()); + + auto* histTruthBase{getHistogram(resultFile, histPath / "hMCTruth")}; + auto* histPrimaryBase{getHistogram(resultFile, histPath / "hPrimary")}; + auto* histSecondaryBase{getHistogram(resultFile, histPath / "hSecondary")}; + + assert(histTruthBase); + assert(histPrimaryBase); + assert(histSecondaryBase); + + TH1* histTruth{histTruthBase}; + TH1* histPrimary{histPrimaryBase}; + TH1* histSecondary{histSecondaryBase}; + + if (projection != "") { + histPrimary = projectHistogram(histPrimaryBase, projection); + histSecondary = projectHistogram(histSecondaryBase, projection); + histTruth = projectHistogram(histTruthBase, projection); + } + + auto* histTotal{cloneHistogram(histPrimary, "hTotal")}; + histTotal->Add(histSecondary); + + auto* histEfficiency{cloneHistogram(histPrimary, "hEfficiency")}; + histEfficiency->Reset(); + setAxisTitles(histEfficiency, projection); + + auto* histWeights{cloneHistogram(histPrimary, "hWeights")}; + histWeights->Reset(); + setAxisTitles(histWeights, projection); + + forEachBin(histPrimary, [&](int x, int y, int z) { + auto primVal{histPrimary->GetBinContent(x, y, z)}; + auto primErr{histPrimary->GetBinError(x, y, z)}; + + auto secVal{histSecondary->GetBinContent(x, y, z)}; + auto secErr{histSecondary->GetBinError(x, y, z)}; + + auto truthVal{histTruth->GetBinContent(x, y, z)}; + auto truthErr{histTruth->GetBinError(x, y, z)}; + + auto effVal{0.}; + auto effErr{0.}; + if (truthVal > 0) { + effVal = primVal / truthVal; + effErr = std::sqrt(std::pow(primErr / truthVal, 2) + std::pow((primVal * truthErr / std::pow(truthVal, 2)), 2)); + } + + histEfficiency->SetBinContent(x, y, z, effVal); + histEfficiency->SetBinError(x, y, z, effErr); + + auto totalVal{primVal + secVal}; + auto totalErr{std::hypot(primErr, secErr)}; + + auto contVal{0.}; + auto contErr{0.}; + if (totalVal > 0) { + contVal = secVal / totalVal; + contErr = std::sqrt(std::pow(secErr / totalVal, 2) + std::pow((secVal * totalErr / std::pow(totalVal, 2)), 2)); + } + + auto weightVal{0.}; + auto weightErr{0.}; + if (effVal > 0) { + weightVal = (1 - contVal) / effVal; + weightErr = std::sqrt(std::pow(contErr / effVal, 2) + std::pow((1 - contVal) * effErr / std::pow(effVal, 2), 2)); + } + + histWeights->SetBinContent(x, y, z, weightVal); + histWeights->SetBinError(x, y, z, weightErr); + }); + + outputFile->WriteTObject(histEfficiency); + outputFile->WriteTObject(histWeights); + + outputFile->Close(); + resultFile->Close(); +} + +auto printUsage(const char* name) -> void +{ + std::cerr << "Usage: " << name << "\n" + << " -f \n" + << " -d \n" + << " -p [optional, default: no projection, 3D histogram]\n" + << " Available projections:\n" + << " x - projection onto pT axis (1D histogram)\n" + << " yx - projection onto pT, eta (2D histogram)\n" + << " zx - projection onto pT, centrality/multiplicity (2D histogram)\n"; +} + +int main(int argc, char** argv) +{ + std::string results{""}, hist{""}, proj{""}; + auto flag{0}; + + while ((flag = getopt(argc, argv, "f:d:p:")) != -1) { + switch (flag) { + case 'f': + results = optarg; + break; + case 'd': + hist = optarg; + break; + case 'p': + proj = optarg; + break; + default: + printUsage(argv[0]); + return 1; + } + } + + if (results.empty() || hist.empty()) { + printUsage(argv[0]); + return 1; + } + + calculateEfficiency(results, hist, proj); + return 0; +} diff --git a/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx b/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx deleted file mode 100644 index 80da583e4c0..00000000000 --- a/PWGCF/FemtoUniverse/Macros/calculateEfficiency.cxx +++ /dev/null @@ -1,149 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace fs = std::filesystem; - -auto* getHistogram(TFile* file, const std::string& name, const std::string& projection) -{ - auto* hist{dynamic_cast(file->Get(name.c_str()))}; - if (projection == "none") { - return hist; - } - return dynamic_cast(hist->Project3D(projection.c_str())); -} - -template -H* cloneHistogram(H* hist, const std::string& name) -{ - return dynamic_cast(hist->Clone(name.c_str())); -} - -void forEachBin(TH1* hist, auto func) -{ - if (hist->GetDimension() == 1) { - for (auto x{1}; x <= hist->GetNbinsX(); ++x) { - func(x, 0, 0); - } - } else if (hist->GetDimension() == 2) { - for (auto x{1}; x <= hist->GetNbinsX(); ++x) { - for (auto y{1}; y <= hist->GetNbinsY(); ++y) { - func(x, y, 0); - } - } - } else if (hist->GetDimension() == 3) { - for (auto x{1}; x <= hist->GetNbinsX(); ++x) { - for (auto y{1}; y <= hist->GetNbinsY(); ++y) { - for (auto z{1}; z <= hist->GetNbinsZ(); ++z) { - func(x, y, z); - } - } - } - } else { - assert(false && "should not happen"); - } -} - -void calculateEfficiency(const fs::path& resultsPath, const fs::path& histPath, const std::string& projection) -{ - assert(projection == "x" || projection == "yx" || projection == "zx" || projection == "none"); - - auto isAlien{false}; - if (resultsPath.string().starts_with("alien://")) { - TGrid::Connect("alien://"); - isAlien = true; - } - - auto* resultFile{TFile::Open(resultsPath.c_str())}; - assert(resultFile != nullptr && !resultFile->IsZombie()); - - using namespace std::chrono; - auto now{duration_cast(system_clock::now().time_since_epoch()).count()}; - auto outputPath{isAlien ? std::filesystem::current_path() : resultsPath.parent_path()}; - outputPath /= std::format("{}-effcor-{}.root", resultsPath.stem().string(), now); - - auto* outputFile{TFile::Open(outputPath.c_str(), "RECREATE")}; - assert(outputFile != nullptr && !outputFile->IsZombie()); - - auto* histTruth{getHistogram(resultFile, histPath / "hMCTruth", projection)}; - assert(histTruth); - - auto* histPrimary{getHistogram(resultFile, histPath / "hPrimary", projection)}; - assert(histPrimary); - - auto* histSecondary{getHistogram(resultFile, histPath / "hSecondary", projection)}; - assert(histSecondary); - - auto* histTotal{cloneHistogram(histPrimary, "hTotal")}; - histTotal->Add(histSecondary); - - auto* histEfficiency{cloneHistogram(histPrimary, "hEfficiency")}; - histEfficiency->Reset(); - - auto* histWeights{cloneHistogram(histPrimary, "hWeights")}; - histWeights->Reset(); - - forEachBin(histPrimary, [&](int x, int y, int z) { - auto primVal{histPrimary->GetBinContent(x, y, z)}; - auto primErr{histPrimary->GetBinError(x, y, z)}; - - auto secVal{histSecondary->GetBinContent(x, y, z)}; - auto secErr{histSecondary->GetBinError(x, y, z)}; - - auto truthVal{histTruth->GetBinContent(x, y, z)}; - auto truthErr{histTruth->GetBinError(x, y, z)}; - - auto effVal{0.}; - auto effErr{0.}; - if (truthVal > 0) { - effVal = primVal / truthVal; - effErr = std::sqrt(std::pow(primErr / truthVal, 2) + std::pow((primVal * truthErr / std::pow(truthVal, 2)), 2)); - } - - histEfficiency->SetBinContent(x, y, z, effVal); - histEfficiency->SetBinError(x, y, z, effErr); - - auto totalVal{primVal + secVal}; - auto totalErr{std::hypot(primErr, secErr)}; - - auto contVal{0.}; - auto contErr{0.}; - if (totalVal > 0) { - contVal = secVal / totalVal; - contErr = std::sqrt(std::pow(secErr / totalVal, 2) + std::pow((secVal * totalErr / std::pow(totalVal, 2)), 2)); - } - - auto weightVal{0.}; - auto weightErr{0.}; - if (effVal > 0) { - weightVal = (1 - contVal) / effVal; - weightErr = std::sqrt(std::pow(contErr / effVal, 2) + std::pow((1 - contVal) * effErr / std::pow(effVal, 2), 2)); - } - - histWeights->SetBinContent(x, y, z, weightVal); - histWeights->SetBinError(x, y, z, weightErr); - }); - - outputFile->WriteTObject(histEfficiency); - outputFile->WriteTObject(histWeights); - - outputFile->Close(); - resultFile->Close(); -} - -int main(int argc, char** argv) -{ - assert(argc == 4); - calculateEfficiency(argv[1], argv[2], argv[3]); - return 0; -} From 7489af0f65c90d50dfd3f138cf983dbad9fca924 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dawid=20Karpi=C5=84ski?= Date: Tue, 20 May 2025 20:41:15 +0200 Subject: [PATCH 8/9] feat: change variables configurable to string --- .../Core/FemtoUniverseEfficiencyCorrection.h | 58 +++++++++++-------- ....cxx => calculateEfficiencyCorrection.cxx} | 19 +++--- 2 files changed, 44 insertions(+), 33 deletions(-) rename PWGCF/FemtoUniverse/Macros/{calculateEffCor.cxx => calculateEfficiencyCorrection.cxx} (94%) diff --git a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h index b8d8cf531bf..b4a1ff97539 100644 --- a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h +++ b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h @@ -29,6 +29,7 @@ #include #include #include +#include #include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h" @@ -48,7 +49,7 @@ struct EffCorConfigurableGroup : framework::ConfigurableGroup { framework::Configurable confEffCorCCDBUrl{"confEffCorCCDBUrl", "http://alice-ccdb.cern.ch", "[Efficiency Correction] CCDB URL to use"}; framework::Configurable confEffCorCCDBPath{"confEffCorCCDBPath", "", "[Efficiency Correction] CCDB path to histograms"}; framework::Configurable> confEffCorCCDBTimestamps{"confEffCorCCDBTimestamps", {}, "[Efficiency Correction] Timestamps of histograms in CCDB (0 can be used as a placeholder, e.g. when running subwagons)"}; - framework::Configurable> confEffCorVariables{"confEffCorVariables", {"pt"}, "[Efficiency Correction] Variables for efficiency correction histogram dimensions (available: pt, eta, cent-mult)"}; + framework::Configurable confEffCorVariables{"confEffCorVariables", "pt", "[Efficiency Correction] Variables for efficiency correction histogram dimensions (available: 'pt'; 'pt,eta'; 'pt,mult'; 'pt,eta,mult')"}; }; class EfficiencyCorrection @@ -66,14 +67,12 @@ class EfficiencyCorrection if (shouldFillHistograms) { for (const auto& suffix : histSuffix) { auto path = std::format("{}/{}", histDirectory, suffix); - - registry->add((path + "/hMCTruth").c_str(), "MCTruth; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); - - registry->add((path + "/hPrimary").c_str(), "Primary; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); - registry->add((path + "/hSecondary").c_str(), "Secondary; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); - registry->add((path + "/hMaterial").c_str(), "Material; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); - registry->add((path + "/hFake").c_str(), "Fake; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); - registry->add((path + "/hOther").c_str(), "Other; pT; Eta; Cent/Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hMCTruth").c_str(), "MCTruth; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hPrimary").c_str(), "Primary; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hSecondary").c_str(), "Secondary; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hMaterial").c_str(), "Material; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hFake").c_str(), "Fake; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hOther").c_str(), "Other; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); } } @@ -97,14 +96,19 @@ class EfficiencyCorrection } if (timestamp > 0) { - if (config->confEffCorVariables->size() == 1) { - hLoaded[idx] = loadHistFromCCDB(timestamp); - } else if (config->confEffCorVariables->size() == 2) { - hLoaded[idx] = loadHistFromCCDB(timestamp); - } else if (config->confEffCorVariables->size() == 3) { - hLoaded[idx] = loadHistFromCCDB(timestamp); - } else { - LOGF(fatal, notify("unknown configuration for efficiency variables")); + switch (getDimensionFromVariables()) { + case 1: + hLoaded[idx] = loadHistFromCCDB(timestamp); + break; + case 2: + hLoaded[idx] = loadHistFromCCDB(timestamp); + break; + case 3: + hLoaded[idx] = loadHistFromCCDB(timestamp); + break; + default: + LOGF(fatal, notify("Unknown configuration for efficiency variables")); + break; } } else { hLoaded[idx] = nullptr; @@ -183,29 +187,29 @@ class EfficiencyCorrection } } - auto getWeight(ParticleNo partNo, auto particle) const -> float + auto getWeight(ParticleNo partNo, auto particle) -> float { auto weight = 1.0f; auto hWeights = hLoaded[partNo - 1]; if (shouldApplyCorrection && hWeights) { auto dim = static_cast(hWeights->GetDimension()); - if (dim != config->confEffCorVariables.value.size()) { + if (dim != getDimensionFromVariables()) { LOGF(fatal, notify("Histogram \"%s\" has wrong dimension %d != %d"), config->confEffCorCCDBPath.value, dim, config->confEffCorVariables.value.size()); return weight; } auto bin = -1; - if (config->confEffCorVariables.value == std::vector{"pt"}) { + if (config->confEffCorVariables.value == "pt") { bin = hLoaded[partNo - 1]->FindBin(particle.pt()); - } else if (config->confEffCorVariables.value == std::vector{"pt", "eta"}) { + } else if (config->confEffCorVariables.value == "pt,eta") { bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.eta()); - } else if (config->confEffCorVariables.value == std::vector{"pt", "cent-mult"}) { + } else if (config->confEffCorVariables.value == "pt,mult") { bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.fdCollision().multV0M()); - } else if (config->confEffCorVariables.value == std::vector{"pt", "eta", "cent-mult"}) { + } else if (config->confEffCorVariables.value == "pt,eta,mult") { bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.eta(), particle.fdCollision().multV0M()); } else { - LOGF(fatal, notify("unknown configuration for efficiency variables")); + LOGF(fatal, notify("Unknown configuration for efficiency variables")); return weight; } weight = hWeights->GetBinContent(bin); @@ -250,6 +254,12 @@ class EfficiencyCorrection return hWeights; } + auto getDimensionFromVariables() -> size_t + { + auto parts = std::views::split(config->confEffCorVariables.value, ','); + return std::ranges::distance(parts); + } + EffCorConfigurableGroup* config{}; bool shouldApplyCorrection{false}; diff --git a/PWGCF/FemtoUniverse/Macros/calculateEffCor.cxx b/PWGCF/FemtoUniverse/Macros/calculateEfficiencyCorrection.cxx similarity index 94% rename from PWGCF/FemtoUniverse/Macros/calculateEffCor.cxx rename to PWGCF/FemtoUniverse/Macros/calculateEfficiencyCorrection.cxx index a5c7625881e..a2f23fa2697 100644 --- a/PWGCF/FemtoUniverse/Macros/calculateEffCor.cxx +++ b/PWGCF/FemtoUniverse/Macros/calculateEfficiencyCorrection.cxx @@ -9,7 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file calculateEffCor.cxx +/// \file calculateEfficiencyCorrection.cxx /// \brief Macro for calculating efficiency corrections based on 3D histograms /// \author Dawid Karpiński, WUT Warsaw, dawid.karpinski@cern.ch @@ -21,13 +21,14 @@ #include #include #include +#include + #include #include #include #include #include #include -#include namespace fs = std::filesystem; @@ -82,17 +83,17 @@ auto setAxisTitles(TH1* hist, const std::string& projection) -> void if (hist->GetDimension() == 2) { if (projection == "yx") { - yAxis->SetTitle("#it{#eta} (rad)"); + yAxis->SetTitle("#it{#eta}"); } else if (projection == "zx") { - yAxis->SetTitle("Cent/Mult"); + yAxis->SetTitle("mult"); } } else if (hist->GetDimension() == 3) { - yAxis->SetTitle("#it{#eta} (rad)"); - zAxis->SetTitle("Cent/Mult"); + yAxis->SetTitle("#it{#eta}"); + zAxis->SetTitle("mult"); } } -auto calculateEfficiency(const fs::path& resultsPath, const fs::path& histPath, const std::string& projection) -> void +auto calculateEfficiencyCorrection(const fs::path& resultsPath, const fs::path& histPath, const std::string& projection) -> void { assert(!resultsPath.empty() && !histPath.empty()); if (projection != "" && projection != "x" && projection != "yx" && projection != "zx") { @@ -204,7 +205,7 @@ auto printUsage(const char* name) -> void << " Available projections:\n" << " x - projection onto pT axis (1D histogram)\n" << " yx - projection onto pT, eta (2D histogram)\n" - << " zx - projection onto pT, centrality/multiplicity (2D histogram)\n"; + << " zx - projection onto pT, mult (2D histogram)\n"; } int main(int argc, char** argv) @@ -234,6 +235,6 @@ int main(int argc, char** argv) return 1; } - calculateEfficiency(results, hist, proj); + calculateEfficiencyCorrection(results, hist, proj); return 0; } From 3ba9c08dec070cbd49d0e6ee1c7dd4a19f85d311 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dawid=20Karpi=C5=84ski?= Date: Tue, 20 May 2025 22:45:11 +0200 Subject: [PATCH 9/9] fix: megalinter --- PWGCF/FemtoUniverse/Macros/calculateEfficiencyCorrection.cxx | 2 +- .../Tasks/femtoUniversePairTaskTrackTrackExtended.cxx | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/PWGCF/FemtoUniverse/Macros/calculateEfficiencyCorrection.cxx b/PWGCF/FemtoUniverse/Macros/calculateEfficiencyCorrection.cxx index a2f23fa2697..e99331231e2 100644 --- a/PWGCF/FemtoUniverse/Macros/calculateEfficiencyCorrection.cxx +++ b/PWGCF/FemtoUniverse/Macros/calculateEfficiencyCorrection.cxx @@ -25,7 +25,7 @@ #include #include -#include +#include // NOLINT #include #include #include diff --git a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx index bb17fc05d2c..b85e4832089 100644 --- a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx +++ b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx @@ -349,8 +349,7 @@ struct FemtoUniversePairTaskTrackTrackExtended { static_cast(confTempFitVarpTBins), {confEtaBins, -2, 2}, confMultBins, - } // - ); + }); } template