Skip to content

Commit 1cc634e

Browse files
BongHwiBongHwi
authored
[PWGCF] Add jFlucEfficiencyTask for efficiency calculation in JCorran (#10066)
1 parent 9a34a1d commit 1cc634e

File tree

2 files changed

+328
-0
lines changed

2 files changed

+328
-0
lines changed

‎PWGCF/JCorran/Tasks/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,3 +38,8 @@ o2physics_add_dpl_workflow(epdzeroflow-analysis
3838
SOURCES jEPDzeroFlowAnalysis.cxx
3939
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore
4040
COMPONENT_NAME Analysis)
41+
42+
o2physics_add_dpl_workflow(j-fluc-efficiency-task
43+
SOURCES jFlucEfficiencyTask.cxx
44+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::JCorran
45+
COMPONENT_NAME Analysis)

‎PWGCF/JCorran/Tasks/jFlucEfficiencyTask.cxx‎

Lines changed: 323 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,323 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
///
12+
/// \file jFlucEfficiencyTask.cxx
13+
/// \brief Task to calculate the efficiency of the cf-derived tracks/particles
14+
/// \author DongJo Kim, Jasper Parkkila, Bong-Hwi Lim (djkim@cern.ch, jparkkil@cern.ch, bong-hwi.lim@cern.ch)
15+
/// \since March 2024
16+
17+
#include <vector>
18+
#include "Framework/AnalysisTask.h"
19+
#include "Framework/HistogramRegistry.h"
20+
#include "Framework/runDataProcessing.h"
21+
#include "Common/Core/TrackSelection.h"
22+
#include "Common/DataModel/Centrality.h"
23+
#include "Common/DataModel/EventSelection.h"
24+
#include "Common/DataModel/Multiplicity.h"
25+
#include "Common/DataModel/PIDResponse.h"
26+
#include "Common/DataModel/TrackSelectionTables.h"
27+
#include "PWGCF/DataModel/CorrelationsDerived.h"
28+
29+
using namespace o2;
30+
using namespace o2::framework;
31+
using namespace o2::framework::expressions;
32+
33+
struct JFlucEfficiencyTask {
34+
// Add the pT binning array as a static member
35+
static constexpr std::array<double, 94> PttJacek = {
36+
0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
37+
0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
38+
1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
39+
2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8,
40+
4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 8.0, 9.0, 10.0,
41+
11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
42+
26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0, 60.0,
43+
70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0,
44+
170.0, 180.0, 190.0, 200.0, 210.0, 220.0, 230.0, 240.0, 250.0, 260.0,
45+
270.0, 280.0, 290.0, 300.0};
46+
47+
// Update the axisPt configuration with proper vector initialization
48+
ConfigurableAxis axisPt{"axisPt", std::vector<double>(PttJacek.begin(), PttJacek.end()), "pT axis"};
49+
50+
// Configurable for track selection
51+
Configurable<float> cfgPtMin{"cfgPtMin", 0.2f, "Minimum transverse momentum"};
52+
Configurable<float> cfgPtMax{"cfgPtMax", 300.0f, "Maximum transverse momentum"};
53+
Configurable<float> cfgEtaMin{"cfgEtaMin", -1.0f, "Minimum pseudorapidity"};
54+
Configurable<float> cfgEtaMax{"cfgEtaMax", 1.0f, "Maximum pseudorapidity"};
55+
Configurable<float> cfgCutVertex{"cfgCutVertex", 10.0f, "Vertex cut"};
56+
Configurable<float> cfgCentMin{"cfgCentMin", 0.0f, "Min centrality"};
57+
Configurable<float> cfgCentMax{"cfgCentMax", 100.0f, "Max centrality"};
58+
Configurable<uint8_t> cfgTrackBitMask{"cfgTrackBitMask", 0, "BitMask for track selection systematics"};
59+
60+
// Configurable axes
61+
ConfigurableAxis axisMultiplicity{"axisMultiplicity", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100}, "multiplicity / centrality axis"};
62+
63+
// Filter declarations
64+
Filter cfCollisionFilter = (nabs(aod::collision::posZ) < cfgCutVertex) &&
65+
(aod::cfcollision::multiplicity > cfgCentMin) &&
66+
(aod::cfcollision::multiplicity < cfgCentMax);
67+
Filter cfMCCollisionFilter = (nabs(aod::mccollision::posZ) < cfgCutVertex) &&
68+
(aod::cfmccollision::multiplicity > cfgCentMin) &&
69+
(aod::cfmccollision::multiplicity < cfgCentMax);
70+
Filter cfMCParticleFilter = (aod::cfmcparticle::pt >= cfgPtMin) &&
71+
(aod::cfmcparticle::pt <= cfgPtMax) &&
72+
(aod::cfmcparticle::eta >= cfgEtaMin) &&
73+
(aod::cfmcparticle::eta <= cfgEtaMax);
74+
Filter cfTrackFilter = (aod::cftrack::pt >= cfgPtMin) &&
75+
(aod::cftrack::pt <= cfgPtMax) &&
76+
(aod::cftrack::eta >= cfgEtaMin) &&
77+
(aod::cftrack::eta <= cfgEtaMax) &&
78+
((aod::track::trackType & (uint8_t)cfgTrackBitMask) == (uint8_t)cfgTrackBitMask);
79+
80+
Filter trackFilter = (nabs(aod::track::eta) < cfgEtaMax) &&
81+
(aod::track::pt > cfgPtMin) &&
82+
((requireGlobalTrackInFilter()) ||
83+
(aod::track::isGlobalTrackSDD == (uint8_t) true));
84+
85+
Configurable<bool> cfgEfficiencyFromData{"cfgEfficiencyFromData", false, "Calculate efficiency using data events as reference"};
86+
Configurable<int> cfgVerbosity{"cfgVerbosity", 0, "Verbosity level"};
87+
88+
// Histogram Registry
89+
HistogramRegistry registry{
90+
"registry",
91+
{{"hEventCounterMC", "Event counter MC;Counter;Counts", {HistType::kTH1F, {{3, -0.5, 2.5}}}},
92+
{"hEventCounterReco", "Event counter Reco;Counter;Counts", {HistType::kTH1F, {{3, -0.5, 2.5}}}},
93+
{"hZVertexMC", "MC Z vertex distribution;Z vertex (cm);Centrality (%)", {HistType::kTH2F, {{200, -20, 20}, {axisMultiplicity}}}},
94+
{"hZVertexReco", "Reconstructed Z vertex distribution;Z vertex (cm);Centrality (%)", {HistType::kTH2F, {{200, -20, 20}, {axisMultiplicity}}}},
95+
{"hZVertexCorrelation", "Z vertex correlation;MC Z vertex (cm);Reco Z vertex (cm)", {HistType::kTH2F, {{200, -20, 20}, {200, -20, 20}}}}}};
96+
97+
// Configurable for debugging
98+
Configurable<bool> debugMode{"debugMode", false, "Debug mode"};
99+
100+
void init(InitContext const&)
101+
{
102+
if (debugMode) {
103+
LOGF(info, "Initializing JFlucEfficiencyTask");
104+
}
105+
106+
if (doprocessMC) {
107+
registry.add("hPtGen", "Generated p_{T} (all);p_{T} (GeV/c);Centrality (%);Counts",
108+
o2::framework::HistType::kTH2F, {AxisSpec(axisPt), AxisSpec(axisMultiplicity)});
109+
registry.add("hEtaGen", "Generated #eta (all);#eta;Centrality (%);Counts",
110+
o2::framework::HistType::kTH2F, {AxisSpec(100, -1, 1), AxisSpec(axisMultiplicity)});
111+
registry.add("hPtGenPos", "Generated p_{T} (positive);p_{T} (GeV/c);Centrality (%);Counts",
112+
o2::framework::HistType::kTH2F, {AxisSpec(axisPt), AxisSpec(axisMultiplicity)});
113+
114+
registry.add("hPtGenNeg", "Generated p_{T} (negative);p_{T} (GeV/c);Centrality (%);Counts",
115+
o2::framework::HistType::kTH2F, {AxisSpec(axisPt), AxisSpec(axisMultiplicity)});
116+
}
117+
118+
if (doprocessData) {
119+
registry.add("hPtRec", "Reconstructed p_{T} (all);p_{T} (GeV/c);Centrality (%);Counts",
120+
o2::framework::HistType::kTH2F, {AxisSpec(axisPt), AxisSpec(axisMultiplicity)});
121+
122+
registry.add("hEtaRec", "Reconstructed #eta (all);#eta;Centrality (%);Counts",
123+
o2::framework::HistType::kTH2F, {AxisSpec(100, -1, 1), AxisSpec(axisMultiplicity)});
124+
125+
registry.add("hPtRecPos", "Reconstructed p_{T} (positive);p_{T} (GeV/c);Centrality (%);Counts",
126+
o2::framework::HistType::kTH2F, {AxisSpec(axisPt), AxisSpec(axisMultiplicity)});
127+
128+
registry.add("hPtRecNeg", "Reconstructed p_{T} (negative);p_{T} (GeV/c);Centrality (%);Counts",
129+
o2::framework::HistType::kTH2F, {AxisSpec(axisPt), AxisSpec(axisMultiplicity)});
130+
}
131+
132+
if (cfgEfficiencyFromData) {
133+
registry.add("hPtGenData", "Generated p_{T} from data events (all);p_{T} (GeV/c);Centrality (%);Counts",
134+
{HistType::kTH2F, {axisPt, axisMultiplicity}});
135+
registry.add("hEtaGenData", "Generated #eta from data events (all);#eta;Centrality (%);Counts",
136+
{HistType::kTH2F, {AxisSpec(100, -1, 1), axisMultiplicity}});
137+
registry.add("hPtGenDataPos", "Generated p_{T} from data events (positive);p_{T} (GeV/c);Centrality (%);Counts",
138+
{HistType::kTH2F, {axisPt, axisMultiplicity}});
139+
registry.add("hPtGenDataNeg", "Generated p_{T} from data events (negative);p_{T} (GeV/c);Centrality (%);Counts",
140+
{HistType::kTH2F, {axisPt, axisMultiplicity}});
141+
registry.add("hPtRecData", "Reconstructed p_{T} (all);p_{T} (GeV/c);Centrality (%);Counts",
142+
o2::framework::HistType::kTH2F, {AxisSpec(axisPt), AxisSpec(axisMultiplicity)});
143+
registry.add("hEtaRecData", "Reconstructed #eta (all);#eta;Centrality (%);Counts",
144+
o2::framework::HistType::kTH2F, {AxisSpec(100, -1, 1), AxisSpec(axisMultiplicity)});
145+
registry.add("hPtRecDataPos", "Reconstructed p_{T} (positive);p_{T} (GeV/c);Centrality (%);Counts",
146+
o2::framework::HistType::kTH2F, {AxisSpec(axisPt), AxisSpec(axisMultiplicity)});
147+
registry.add("hPtRecDataNeg", "Reconstructed p_{T} (negative);p_{T} (GeV/c);Centrality (%);Counts",
148+
o2::framework::HistType::kTH2F, {AxisSpec(axisPt), AxisSpec(axisMultiplicity)});
149+
}
150+
151+
// Initialize histogram labels
152+
auto h1 = registry.get<TH1>(HIST("hEventCounterMC"));
153+
auto h2 = registry.get<TH1>(HIST("hEventCounterReco"));
154+
155+
if (h1 && h2) {
156+
h1->GetXaxis()->SetBinLabel(1, "All MC Events");
157+
h1->GetXaxis()->SetBinLabel(2, "Selected MC Events");
158+
h1->GetXaxis()->SetBinLabel(3, "Analyzed MC Events");
159+
160+
h2->GetXaxis()->SetBinLabel(1, "All Reco Events");
161+
h2->GetXaxis()->SetBinLabel(2, "Selected Reco Events");
162+
h2->GetXaxis()->SetBinLabel(3, "Analyzed Reco Events");
163+
} else {
164+
LOGF(error, "Failed to get histograms from registry");
165+
}
166+
}
167+
168+
void processMC(soa::Filtered<aod::CFMcCollisions>::iterator const& mcCollision, soa::Filtered<aod::CFMcParticles> const& mcParticles)
169+
{
170+
float centrality = mcCollision.multiplicity();
171+
172+
for (const auto& particle : mcParticles) {
173+
if (!particle.isPhysicalPrimary()) {
174+
continue;
175+
}
176+
177+
registry.fill(HIST("hPtGen"), particle.pt(), centrality);
178+
registry.fill(HIST("hEtaGen"), particle.eta(), centrality);
179+
180+
if (particle.sign() > 0) { // Positive particles
181+
registry.fill(HIST("hPtGenPos"), particle.pt(), centrality);
182+
} else if (particle.sign() < 0) { // Negative particles
183+
registry.fill(HIST("hPtGenNeg"), particle.pt(), centrality);
184+
}
185+
}
186+
}
187+
188+
void processData(soa::Filtered<aod::CFCollisions>::iterator const& cfCollision, soa::Filtered<aod::CFTracks> const& cfTracks)
189+
{
190+
float centrality = cfCollision.multiplicity();
191+
192+
if (centrality < cfgCentMin || centrality > cfgCentMax) {
193+
return;
194+
}
195+
196+
for (const auto& track : cfTracks) {
197+
registry.fill(HIST("hPtRec"), track.pt(), centrality);
198+
registry.fill(HIST("hEtaRec"), track.eta(), centrality);
199+
200+
if (track.sign() > 0) { // Positive tracks
201+
registry.fill(HIST("hPtRecPos"), track.pt(), centrality);
202+
} else if (track.sign() < 0) { // Negative tracks
203+
registry.fill(HIST("hPtRecNeg"), track.pt(), centrality);
204+
}
205+
}
206+
}
207+
208+
template <typename TCollision, typename TTracks>
209+
void fillQA(const TCollision& /*collision*/, float multiplicity, const TTracks& tracks)
210+
{
211+
registry.fill(HIST("multiplicity"), multiplicity);
212+
for (const auto& track : tracks) {
213+
registry.fill(HIST("yields"), multiplicity, track.pt(), track.eta());
214+
registry.fill(HIST("etaphi"), track.eta(), track.phi());
215+
}
216+
}
217+
218+
// NOTE SmallGroups includes soa::Filtered always
219+
Preslice<aod::CFTracksWithLabel> perCollision = aod::cftrack::cfCollisionId;
220+
void processEfficiency(soa::Filtered<aod::CFMcCollisions>::iterator const& mcCollision,
221+
aod::CFMcParticles const& mcParticles,
222+
soa::SmallGroups<aod::CFCollisionsWithLabel> const& collisions,
223+
aod::CFTracksWithLabel const& tracks)
224+
{
225+
try {
226+
// Count MC events and fill MC z-vertex with centrality
227+
registry.fill(HIST("hEventCounterMC"), 0);
228+
registry.fill(HIST("hZVertexMC"), mcCollision.posZ(), mcCollision.multiplicity());
229+
230+
if (debugMode) {
231+
LOGF(info, "Processing MC collision %d at z = %.3f", mcCollision.globalIndex(), mcCollision.posZ());
232+
}
233+
234+
// Fill MC particle histograms
235+
for (const auto& mcParticle : mcParticles) {
236+
if (!mcParticle.isPhysicalPrimary())
237+
continue;
238+
239+
// Fill generated particle histograms
240+
registry.fill(HIST("hPtGen"), mcParticle.pt(), mcCollision.multiplicity());
241+
registry.fill(HIST("hEtaGen"), mcParticle.eta(), mcCollision.multiplicity());
242+
243+
if (mcParticle.sign() > 0) {
244+
registry.fill(HIST("hPtGenPos"), mcParticle.pt(), mcCollision.multiplicity());
245+
} else if (mcParticle.sign() < 0) {
246+
registry.fill(HIST("hPtGenNeg"), mcParticle.pt(), mcCollision.multiplicity());
247+
}
248+
249+
if (cfgEfficiencyFromData) {
250+
registry.fill(HIST("hPtGenData"), mcParticle.pt(), mcCollision.multiplicity());
251+
registry.fill(HIST("hEtaGenData"), mcParticle.eta(), mcCollision.multiplicity());
252+
if (mcParticle.sign() > 0) {
253+
registry.fill(HIST("hPtGenDataPos"), mcParticle.pt(), mcCollision.multiplicity());
254+
} else if (mcParticle.sign() < 0) {
255+
registry.fill(HIST("hPtGenDataNeg"), mcParticle.pt(), mcCollision.multiplicity());
256+
}
257+
}
258+
}
259+
260+
registry.fill(HIST("hEventCounterMC"), 1);
261+
262+
// Check reconstructed collisions
263+
if (collisions.size() == 0) {
264+
if (debugMode) {
265+
LOGF(info, "No reconstructed collisions found for MC collision %d", mcCollision.globalIndex());
266+
}
267+
return;
268+
}
269+
270+
// Process reconstructed events
271+
for (const auto& collision : collisions) {
272+
registry.fill(HIST("hEventCounterReco"), 0);
273+
registry.fill(HIST("hZVertexReco"), collision.posZ(), collision.multiplicity());
274+
registry.fill(HIST("hZVertexCorrelation"), mcCollision.posZ(), collision.posZ());
275+
276+
if (debugMode) {
277+
LOGF(info, "Processing reconstructed collision %d at z = %.3f",
278+
collision.globalIndex(), collision.posZ());
279+
}
280+
281+
// Fill track histograms
282+
for (const auto& track : tracks) {
283+
if (!track.has_cfMCParticle()) {
284+
if (debugMode) {
285+
LOGF(debug, "Track without MC particle found");
286+
}
287+
continue;
288+
}
289+
290+
registry.fill(HIST("hPtRecData"), track.pt(), collision.multiplicity());
291+
registry.fill(HIST("hEtaRecData"), track.eta(), collision.multiplicity());
292+
293+
if (track.sign() > 0) {
294+
registry.fill(HIST("hPtRecDataPos"), track.pt(), collision.multiplicity());
295+
} else if (track.sign() < 0) {
296+
registry.fill(HIST("hPtRecDatfaNeg"), track.pt(), collision.multiplicity());
297+
}
298+
}
299+
300+
// Count selected and analyzed events
301+
registry.fill(HIST("hEventCounterReco"), 1);
302+
registry.fill(HIST("hEventCounterReco"), 2);
303+
}
304+
305+
registry.fill(HIST("hEventCounterMC"), 2);
306+
307+
} catch (const std::exception& e) {
308+
LOGF(error, "Exception caught in processEfficiency: %s", e.what());
309+
} catch (...) {
310+
LOGF(error, "Unknown exception caught in processEfficiency");
311+
}
312+
}
313+
314+
PROCESS_SWITCH(JFlucEfficiencyTask, processMC, "Process MC only", false);
315+
PROCESS_SWITCH(JFlucEfficiencyTask, processData, "Process data only", false);
316+
PROCESS_SWITCH(JFlucEfficiencyTask, processEfficiency, "Process efficiency task", true);
317+
};
318+
319+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
320+
{
321+
return WorkflowSpec{
322+
adaptAnalysisTask<JFlucEfficiencyTask>(cfgc)};
323+
}

0 commit comments

Comments
 (0)