Skip to content

Commit 07588b1

Browse files
davkkdavkk
authored
[PWGCF/FemtoUniverse] Efficiency correction in 3D (#11270)
1 parent 8c73854 commit 07588b1

File tree

6 files changed

+435
-74
lines changed

6 files changed

+435
-74
lines changed

‎PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h‎

Lines changed: 156 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -16,15 +16,22 @@
1616
#ifndef PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSEEFFICIENCYCORRECTION_H_
1717
#define PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSEEFFICIENCYCORRECTION_H_
1818

19+
#include <CCDB/BasicCCDBManager.h>
20+
#include <Framework/Configurable.h>
21+
#include <Framework/ConfigurableKinds.h>
22+
#include <Framework/HistogramRegistry.h>
23+
#include <Framework/HistogramSpec.h>
24+
25+
#include <Framework/O2DatabasePDGPlugin.h>
26+
#include <TH1.h>
27+
#include <TH2.h>
28+
1929
#include <vector>
2030
#include <string>
2131
#include <algorithm>
32+
#include <ranges>
2233

23-
#include "Framework/Configurable.h"
24-
#include "CCDB/BasicCCDBManager.h"
25-
#include "TH1.h"
26-
#include "TH2.h"
27-
#include "TH3.h"
34+
#include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h"
2835

2936
namespace o2::analysis::femto_universe::efficiency_correction
3037
{
@@ -34,39 +41,41 @@ enum ParticleNo : size_t {
3441
};
3542

3643
template <size_t T>
37-
concept isOneOrTwo = T == ParticleNo::ONE || T == ParticleNo::TWO;
38-
39-
template <typename T>
40-
consteval auto getHistDim() -> int
41-
{
42-
if (std::is_same_v<T, TH1>)
43-
return 1;
44-
else if (std::is_same_v<T, TH2>)
45-
return 2;
46-
else if (std::is_same_v<T, TH3>)
47-
return 3;
48-
else
49-
return -1;
50-
}
44+
concept IsOneOrTwo = T == ParticleNo::ONE || T == ParticleNo::TWO;
5145

5246
struct EffCorConfigurableGroup : framework::ConfigurableGroup {
5347
framework::Configurable<bool> confEffCorApply{"confEffCorApply", false, "[Efficiency Correction] Should apply efficiency corrections"};
48+
framework::Configurable<bool> confEffCorFillHist{"confEffCorFillHist", false, "[Efficiency Correction] Should fill histograms for efficiency corrections"};
5449
framework::Configurable<std::string> confEffCorCCDBUrl{"confEffCorCCDBUrl", "http://alice-ccdb.cern.ch", "[Efficiency Correction] CCDB URL to use"};
5550
framework::Configurable<std::string> confEffCorCCDBPath{"confEffCorCCDBPath", "", "[Efficiency Correction] CCDB path to histograms"};
5651
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)"};
52+
framework::Configurable<std::string> confEffCorVariables{"confEffCorVariables", "pt", "[Efficiency Correction] Variables for efficiency correction histogram dimensions (available: 'pt'; 'pt,eta'; 'pt,mult'; 'pt,eta,mult')"};
5753
};
5854

59-
template <typename HistType>
60-
requires std::is_base_of_v<TH1, HistType>
6155
class EfficiencyCorrection
6256
{
6357
public:
6458
explicit EfficiencyCorrection(EffCorConfigurableGroup* config) : config(config) // o2-linter: disable=name/function-variable
6559
{
6660
}
6761

68-
auto init() -> void
62+
auto init(framework::HistogramRegistry* registry, std::vector<framework::AxisSpec> axisSpecs) -> void
6963
{
64+
shouldFillHistograms = config->confEffCorFillHist;
65+
66+
histRegistry = registry;
67+
if (shouldFillHistograms) {
68+
for (const auto& suffix : histSuffix) {
69+
auto path = std::format("{}/{}", histDirectory, suffix);
70+
registry->add((path + "/hMCTruth").c_str(), "MCTruth; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs);
71+
registry->add((path + "/hPrimary").c_str(), "Primary; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs);
72+
registry->add((path + "/hSecondary").c_str(), "Secondary; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs);
73+
registry->add((path + "/hMaterial").c_str(), "Material; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs);
74+
registry->add((path + "/hFake").c_str(), "Fake; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs);
75+
registry->add((path + "/hOther").c_str(), "Other; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs);
76+
}
77+
}
78+
7079
ccdb.setURL(config->confEffCorCCDBUrl);
7180
ccdb.setLocalObjectValidityChecking();
7281
ccdb.setFatalWhenNull(false);
@@ -86,20 +95,123 @@ class EfficiencyCorrection
8695
continue;
8796
}
8897

89-
hLoaded[idx] = timestamp > 0 ? loadHistFromCCDB(timestamp) : nullptr;
98+
if (timestamp > 0) {
99+
switch (getDimensionFromVariables()) {
100+
case 1:
101+
hLoaded[idx] = loadHistFromCCDB<TH1>(timestamp);
102+
break;
103+
case 2:
104+
hLoaded[idx] = loadHistFromCCDB<TH2>(timestamp);
105+
break;
106+
case 3:
107+
hLoaded[idx] = loadHistFromCCDB<TH3>(timestamp);
108+
break;
109+
default:
110+
LOGF(fatal, notify("Unknown configuration for efficiency variables"));
111+
break;
112+
}
113+
} else {
114+
hLoaded[idx] = nullptr;
115+
}
90116
}
91117
}
92118
}
93119

94-
template <typename... BinVars>
95-
requires(sizeof...(BinVars) == getHistDim<HistType>())
96-
auto getWeight(ParticleNo partNo, const BinVars&... binVars) const -> float
120+
template <uint8_t N>
121+
requires IsOneOrTwo<N>
122+
void fillTruthHist(auto particle)
123+
{
124+
if (!shouldFillHistograms) {
125+
return;
126+
}
127+
128+
histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hMCTruth"),
129+
particle.pt(),
130+
particle.eta(),
131+
particle.fdCollision().multV0M());
132+
}
133+
134+
template <uint8_t N>
135+
requires IsOneOrTwo<N>
136+
void fillRecoHist(auto particle, int particlePDG)
137+
{
138+
if (!shouldFillHistograms) {
139+
return;
140+
}
141+
142+
if (!particle.has_fdMCParticle()) {
143+
return;
144+
}
145+
146+
auto mcParticle = particle.fdMCParticle();
147+
148+
if (mcParticle.pdgMCTruth() == particlePDG) {
149+
switch (mcParticle.partOriginMCTruth()) {
150+
case (o2::aod::femtouniverse_mc_particle::kPrimary):
151+
histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hPrimary"),
152+
mcParticle.pt(),
153+
mcParticle.eta(),
154+
particle.fdCollision().multV0M());
155+
break;
156+
157+
case (o2::aod::femtouniverse_mc_particle::kDaughter):
158+
case (o2::aod::femtouniverse_mc_particle::kDaughterLambda):
159+
case (o2::aod::femtouniverse_mc_particle::kDaughterSigmaplus):
160+
histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hSecondary"),
161+
mcParticle.pt(),
162+
mcParticle.eta(),
163+
particle.fdCollision().multV0M());
164+
break;
165+
166+
case (o2::aod::femtouniverse_mc_particle::kMaterial):
167+
histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hMaterial"),
168+
mcParticle.pt(),
169+
mcParticle.eta(),
170+
particle.fdCollision().multV0M());
171+
break;
172+
173+
case (o2::aod::femtouniverse_mc_particle::kFake):
174+
histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hFake"),
175+
mcParticle.pt(),
176+
mcParticle.eta(),
177+
particle.fdCollision().multV0M());
178+
break;
179+
180+
default:
181+
histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hOther"),
182+
mcParticle.pt(),
183+
mcParticle.eta(),
184+
particle.fdCollision().multV0M());
185+
break;
186+
}
187+
}
188+
}
189+
190+
auto getWeight(ParticleNo partNo, auto particle) -> float
97191
{
98192
auto weight = 1.0f;
99193
auto hWeights = hLoaded[partNo - 1];
100194

101195
if (shouldApplyCorrection && hWeights) {
102-
auto bin = hWeights->FindBin(binVars...);
196+
auto dim = static_cast<size_t>(hWeights->GetDimension());
197+
if (dim != getDimensionFromVariables()) {
198+
LOGF(fatal, notify("Histogram \"%s\" has wrong dimension %d != %d"), config->confEffCorCCDBPath.value, dim, config->confEffCorVariables.value.size());
199+
return weight;
200+
}
201+
202+
auto bin = -1;
203+
if (config->confEffCorVariables.value == "pt") {
204+
bin = hLoaded[partNo - 1]->FindBin(particle.pt());
205+
} else if (config->confEffCorVariables.value == "pt,eta") {
206+
bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.eta());
207+
} else if (config->confEffCorVariables.value == "pt,mult") {
208+
bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.fdCollision().multV0M());
209+
} else if (config->confEffCorVariables.value == "pt,eta,mult") {
210+
bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.eta(), particle.fdCollision().multV0M());
211+
} else {
212+
LOGF(fatal, notify("Unknown configuration for efficiency variables"));
213+
return weight;
214+
}
103215
weight = hWeights->GetBinContent(bin);
104216
}
105217

@@ -112,7 +224,7 @@ class EfficiencyCorrection
112224
return fmt::format("[EFFICIENCY CORRECTION] {}", msg);
113225
}
114226

115-
static auto isHistEmpty(HistType* hist) -> bool
227+
static auto isHistEmpty(TH1* hist) -> bool
116228
{
117229
if (!hist) {
118230
return true;
@@ -125,9 +237,10 @@ class EfficiencyCorrection
125237
return true;
126238
}
127239

128-
auto loadHistFromCCDB(const int64_t timestamp) const -> HistType*
240+
template <typename H>
241+
auto loadHistFromCCDB(const int64_t timestamp) const -> H*
129242
{
130-
auto hWeights = ccdb.getForTimeStamp<HistType>(config->confEffCorCCDBPath, timestamp);
243+
auto hWeights = ccdb.getForTimeStamp<H>(config->confEffCorCCDBPath, timestamp);
131244
if (!hWeights || hWeights->IsZombie()) {
132245
LOGF(error, notify("Could not load histogram \"%s/%ld\""), config->confEffCorCCDBPath.value, timestamp);
133246
return nullptr;
@@ -141,12 +254,23 @@ class EfficiencyCorrection
141254
return hWeights;
142255
}
143256

257+
auto getDimensionFromVariables() -> size_t
258+
{
259+
auto parts = std::views::split(config->confEffCorVariables.value, ',');
260+
return std::ranges::distance(parts);
261+
}
262+
144263
EffCorConfigurableGroup* config{};
145264

146-
bool shouldApplyCorrection = false;
265+
bool shouldApplyCorrection{false};
266+
bool shouldFillHistograms{false};
147267

148268
o2::ccdb::BasicCCDBManager& ccdb{o2::ccdb::BasicCCDBManager::instance()};
149-
std::array<HistType*, 2> hLoaded{};
269+
std::array<TH1*, 2> hLoaded{};
270+
271+
framework::HistogramRegistry* histRegistry{};
272+
static constexpr std::string_view histDirectory{"EfficiencyCorrection"};
273+
static constexpr std::string_view histSuffix[2]{"one", "two"};
150274
};
151275

152276
} // namespace o2::analysis::femto_universe::efficiency_correction

0 commit comments

Comments
 (0)