Skip to content

Commit 1229bbe

Browse files
zhangbiao-phyzhangbiao-phy
authored
Add optional event-status and generated-centrality filling to taskD0
Added event status handling and centrality tracking for reconstructed collisions.
1 parent f55863b commit 1229bbe

1 file changed

Lines changed: 83 additions & 18 deletions

File tree

‎PWGHF/D2H/Tasks/taskD0.cxx‎

Lines changed: 83 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@
2626
#include "PWGHF/DataModel/TrackIndexSkimmingTables.h"
2727
#include "PWGHF/Utils/utilsEvSelHf.h"
2828
#include "PWGHF/Utils/utilsUpcHf.h"
29+
#include "PWGLF/DataModel/mcCentrality.h"
2930
#include "PWGUD/Core/UPCHelpers.h"
3031

3132
#include "Common/CCDB/ctpRateFetcher.h"
@@ -58,6 +59,7 @@
5859
#include <cstdint>
5960
#include <numeric>
6061
#include <string>
62+
#include <unordered_map>
6163
#include <vector>
6264

6365
using namespace o2;
@@ -82,6 +84,12 @@ enum CandTypeSel {
8284
};
8385
} // namespace
8486
struct HfTaskD0 {
87+
enum RecoEventStatus {
88+
NoRecoCollision = 0,
89+
NoSelectedRecoCollision,
90+
SelectedRecoCollision
91+
};
92+
8593
Configurable<int> selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"};
8694
Configurable<int> selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"};
8795
Configurable<double> yCandGenMax{"yCandGenMax", 0.5, "max. gen particle rapidity"};
@@ -92,6 +100,8 @@ struct HfTaskD0 {
92100
Configurable<int> selectionPid{"selectionPid", 1, "Selection Flag for reco PID candidates"};
93101
Configurable<std::vector<double>> binsPt{"binsPt", std::vector<double>{hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits"};
94102
Configurable<int> centEstimator{"centEstimator", 0, "Centrality estimation (None: 0, FT0C: 2, FT0M: 3)"};
103+
Configurable<bool> fillEventStatus{"fillEventStatus", false, "Use generated centrality in hSparseAcc and append reco-event status (requires storeCentrality)"};
104+
Configurable<int> centEstimatorGen{"centEstimatorGen", 2, "Generated centrality for signal loss (FT0C: 2, FT0M: 3); use the calibrated OO estimator"};
95105
Configurable<int> occEstimator{"occEstimator", 0, "Occupancy estimation (None: 0, ITS: 1, FT0C: 2)"};
96106
Configurable<bool> storeCentrality{"storeCentrality", false, "Flag to store centrality information"};
97107
Configurable<bool> storeOccupancyAndIR{"storeOccupancyAndIR", false, "Flag to store occupancy information and interaction rate"};
@@ -125,6 +135,7 @@ struct HfTaskD0 {
125135
using CollisionsCent = soa::Join<aod::Collisions, aod::EvSels, aod::PVMults, aod::CentFT0Ms, aod::CentFT0Cs>;
126136
using CollisionsWithMcLabels = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::PVMults>;
127137
using CollisionsWithMcLabelsCent = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::PVMults, aod::CentFT0Ms, aod::CentFT0Cs>;
138+
using McCollisionsWithCentrality = soa::Join<aod::McCollisions, aod::McCentFT0Cs, aod::McCentFT0Ms>;
128139
using TracksSelQuality = soa::Join<aod::TracksExtra, aod::TracksWMc>;
129140
using TracksWPid = soa::Join<o2::aod::FullTracks, aod::TracksDCA, o2::aod::TrackSelection, aod::TracksPidPi, aod::PidTpcTofFullPi, aod::TracksPidKa, aod::PidTpcTofFullKa, aod::TracksPidPr, aod::PidTpcTofFullPr>;
130141
// using TracksWithExtra = o2::soa::Join<o2::aod::FullTracks, o2::aod::TrackSelection>;
@@ -331,13 +342,22 @@ struct HfTaskD0 {
331342
std::vector<AxisSpec> axesAcc = {thnAxisGenPtD, thnAxisGenPtB, thnAxisY, thnAxisOrigin, thnAxisNumPvContr};
332343

333344
if (storeCentrality) {
334-
axesAcc.push_back(thnAxisCent);
345+
axesAcc.emplace_back(thnConfigAxisCent, fillEventStatus ? "Generated centrality (%)" : "Centrality");
335346
}
336347
// interaction rate only store in Data and MC Reco. Level
337348
if (storeOccupancyAndIR) {
338349
axesAcc.push_back(thnAxisOccupancy);
339350
}
340351

352+
if (fillEventStatus) {
353+
if (!storeCentrality) {
354+
LOGP(fatal, "fillEventStatus requires storeCentrality=true.");
355+
}
356+
if (centEstimatorGen != CentralityEstimator::FT0C && centEstimatorGen != CentralityEstimator::FT0M) {
357+
LOGP(fatal, "centEstimatorGen must be FT0C (2) or FT0M (3).");
358+
}
359+
axesAcc.emplace_back(3, -0.5, 2.5, "Reco event status (0: no reco, 1: none selected, 2: selected)");
360+
}
341361
registry.add("hSparseAcc", "Thn for generated D0 from charm and beauty", HistType::kTHnSparseD, axesAcc);
342362
registry.get<THnSparse>(HIST("hSparseAcc"))->Sumw2();
343363
}
@@ -484,7 +504,11 @@ struct HfTaskD0 {
484504
registry.add("QAtracks/hDCAxy_GapC", "Gap C; DCA xy", {HistType::kTH1F, {{400, -2, 2.}}});
485505
registry.add("QAtracks/hDCAz_GapC", "Gap C; DCA z", {HistType::kTH1F, {{400, -4, 4.}}});
486506

487-
hfEvSel.addHistograms(registry);
507+
if (fillEventStatus && (doprocessMcWithDCAFitterN || doprocessMcWithDCAFitterNCent || doprocessMcWithKFParticle || doprocessMcWithDCAFitterNMl || doprocessMcWithDCAFitterNMlCent || doprocessMcWithKFParticleMl)) {
508+
hfEvSel.init(registry);
509+
} else {
510+
hfEvSel.addHistograms(registry);
511+
}
488512

489513
ccdb->setURL(ccdbUrl);
490514
ccdb->setCaching(true);
@@ -970,7 +994,7 @@ struct HfTaskD0 {
970994
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
971995
TracksSelQuality const&,
972996
CollType const& collisions,
973-
aod::McCollisions const&,
997+
McCollisionsWithCentrality const&,
974998
BCsType const&)
975999
{
9761000
// MC rec.
@@ -1255,7 +1279,31 @@ struct HfTaskD0 {
12551279
}
12561280
}
12571281
}
1258-
// MC gen.
1282+
// Classify reco associations once per collision, not once per generated particle.
1283+
// The task's hfEvSel.* configuration must match the candidate creator's cuts.
1284+
std::unordered_map<int64_t, int> recoEventStatus;
1285+
if (fillEventStatus) {
1286+
for (const auto& collision : collisions) {
1287+
if (!collision.has_mcCollision()) {
1288+
continue;
1289+
}
1290+
float unusedCentrality{-1.f};
1291+
const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask<true, CentralityEstimator::None, BCsType>(collision, unusedCentrality, ccdb, registry);
1292+
bool selected = rejectionMask == 0;
1293+
if constexpr (std::is_same_v<CollType, CollisionsWithMcLabelsCent>) {
1294+
if (centEstimator != CentralityEstimator::None) {
1295+
const auto recoCentrality = getCentralityColl(collision, centEstimator);
1296+
selected = selected && recoCentrality >= hfEvSel.centralityMin && recoCentrality <= hfEvSel.centralityMax;
1297+
}
1298+
}
1299+
auto& status = recoEventStatus[collision.mcCollisionId()];
1300+
const int currentStatus = selected ? SelectedRecoCollision : NoSelectedRecoCollision;
1301+
status = std::max(status, currentStatus);
1302+
}
1303+
}
1304+
1305+
// MC gen. Each particle is filled once, even for split or unreconstructed events.
1306+
// Upstream event rejection can clear flagMcMatchGen: this is not an unfiltered truth sample.
12591307
for (const auto& particle : mcParticles) {
12601308
if (std::abs(particle.flagMcMatchGen()) == o2::hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK) {
12611309
if (yCandGenMax >= 0. && std::abs(RecoDecay::y(particle.pVector(), o2::constants::physics::MassD0)) > yCandGenMax) {
@@ -1267,6 +1315,8 @@ struct HfTaskD0 {
12671315
registry.fill(HIST("hPtGen"), ptGen);
12681316
registry.fill(HIST("hPtVsYGen"), ptGen, yGen);
12691317

1318+
int eventStatus = NoRecoCollision;
1319+
12701320
unsigned maxNumContrib = 0;
12711321
float cent{-1.f};
12721322
float occ{-1.f};
@@ -1288,32 +1338,47 @@ struct HfTaskD0 {
12881338
}
12891339
}
12901340

1341+
if (fillEventStatus) {
1342+
const auto mcCollision = particle.template mcCollision_as<McCollisionsWithCentrality>();
1343+
const auto eventEntry = recoEventStatus.find(mcCollision.globalIndex());
1344+
eventStatus = eventEntry == recoEventStatus.end() ? NoRecoCollision : eventEntry->second;
1345+
// Replace the existing centrality coordinate, including zero-reco events.
1346+
cent = centEstimatorGen == CentralityEstimator::FT0C ? mcCollision.centFT0C() : mcCollision.centFT0M();
1347+
}
1348+
const auto fillGeneratedSparse = [&](auto... coordinates) {
1349+
if (fillEventStatus) {
1350+
registry.fill(HIST("hSparseAcc"), coordinates..., eventStatus);
1351+
} else {
1352+
registry.fill(HIST("hSparseAcc"), coordinates...);
1353+
}
1354+
};
1355+
12911356
if (particle.originMcGen() == RecoDecay::OriginType::Prompt) {
12921357
registry.fill(HIST("hPtGenPrompt"), ptGen);
12931358
registry.fill(HIST("hYGenPrompt"), yGen);
12941359
registry.fill(HIST("hPtVsYGenPrompt"), ptGen, yGen);
12951360
if (storeCentrality && storeOccupancyAndIR) {
1296-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 1, maxNumContrib, cent, occ);
1361+
fillGeneratedSparse(ptGen, ptGenB, yGen, 1, maxNumContrib, cent, occ);
12971362
} else if (storeCentrality && !storeOccupancyAndIR) {
1298-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 1, maxNumContrib, cent);
1363+
fillGeneratedSparse(ptGen, ptGenB, yGen, 1, maxNumContrib, cent);
12991364
} else if (!storeCentrality && storeOccupancyAndIR) {
1300-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 1, maxNumContrib, occ);
1365+
fillGeneratedSparse(ptGen, ptGenB, yGen, 1, maxNumContrib, occ);
13011366
} else {
1302-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 1, maxNumContrib);
1367+
fillGeneratedSparse(ptGen, ptGenB, yGen, 1, maxNumContrib);
13031368
}
13041369
} else {
13051370
ptGenB = mcParticles.rawIteratorAt(particle.idxBhadMotherPart()).pt();
13061371
registry.fill(HIST("hPtGenNonPrompt"), ptGen);
13071372
registry.fill(HIST("hYGenNonPrompt"), yGen);
13081373
registry.fill(HIST("hPtVsYGenNonPrompt"), ptGen, yGen);
13091374
if (storeCentrality && storeOccupancyAndIR) {
1310-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 2, maxNumContrib, cent, occ);
1375+
fillGeneratedSparse(ptGen, ptGenB, yGen, 2, maxNumContrib, cent, occ);
13111376
} else if (storeCentrality && !storeOccupancyAndIR) {
1312-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 2, maxNumContrib, cent);
1377+
fillGeneratedSparse(ptGen, ptGenB, yGen, 2, maxNumContrib, cent);
13131378
} else if (!storeCentrality && storeOccupancyAndIR) {
1314-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 2, maxNumContrib, occ);
1379+
fillGeneratedSparse(ptGen, ptGenB, yGen, 2, maxNumContrib, occ);
13151380
} else {
1316-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 2, maxNumContrib);
1381+
fillGeneratedSparse(ptGen, ptGenB, yGen, 2, maxNumContrib);
13171382
}
13181383
}
13191384
registry.fill(HIST("hEtaGen"), particle.eta());
@@ -1325,7 +1390,7 @@ struct HfTaskD0 {
13251390
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13261391
TracksSelQuality const& tracks,
13271392
CollisionsWithMcLabels const& collisions,
1328-
aod::McCollisions const& mcCollisions,
1393+
McCollisionsWithCentrality const& mcCollisions,
13291394
aod::BcFullInfos const& bcs)
13301395
{
13311396
processMc<aod::hf_cand::VertexerType::DCAFitter, false>(selectedD0CandidatesMc, mcParticles, tracks, collisions, mcCollisions, bcs);
@@ -1336,7 +1401,7 @@ struct HfTaskD0 {
13361401
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13371402
TracksSelQuality const& tracks,
13381403
CollisionsWithMcLabelsCent const& collisions,
1339-
aod::McCollisions const& mcCollisions,
1404+
McCollisionsWithCentrality const& mcCollisions,
13401405
aod::BcFullInfos const& bcs)
13411406
{
13421407
processMc<aod::hf_cand::VertexerType::DCAFitter, false>(selectedD0CandidatesMc, mcParticles, tracks, collisions, mcCollisions, bcs);
@@ -1347,7 +1412,7 @@ struct HfTaskD0 {
13471412
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13481413
TracksSelQuality const& tracks,
13491414
CollisionsWithMcLabels const& collisions,
1350-
aod::McCollisions const& mcCollisions,
1415+
McCollisionsWithCentrality const& mcCollisions,
13511416
aod::BcFullInfos const& bcs)
13521417
{
13531418
processMc<aod::hf_cand::VertexerType::KfParticle, false>(selectedD0CandidatesMcKF, mcParticles, tracks, collisions, mcCollisions, bcs);
@@ -1359,7 +1424,7 @@ struct HfTaskD0 {
13591424
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13601425
TracksSelQuality const& tracks,
13611426
CollisionsWithMcLabels const& collisions,
1362-
aod::McCollisions const& mcCollisions,
1427+
McCollisionsWithCentrality const& mcCollisions,
13631428
aod::BcFullInfos const& bcs)
13641429
{
13651430
processMc<aod::hf_cand::VertexerType::DCAFitter, true>(selectedD0CandidatesMlMc, mcParticles, tracks, collisions, mcCollisions, bcs);
@@ -1370,7 +1435,7 @@ struct HfTaskD0 {
13701435
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13711436
TracksSelQuality const& tracks,
13721437
CollisionsWithMcLabelsCent const& collisions,
1373-
aod::McCollisions const& mcCollisions,
1438+
McCollisionsWithCentrality const& mcCollisions,
13741439
aod::BcFullInfos const& bcs)
13751440
{
13761441
processMc<aod::hf_cand::VertexerType::DCAFitter, true>(selectedD0CandidatesMlMc, mcParticles, tracks, collisions, mcCollisions, bcs);
@@ -1381,7 +1446,7 @@ struct HfTaskD0 {
13811446
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13821447
TracksSelQuality const& tracks,
13831448
CollisionsWithMcLabels const& collisions,
1384-
aod::McCollisions const& mcCollisions,
1449+
McCollisionsWithCentrality const& mcCollisions,
13851450
aod::BcFullInfos const& bcs)
13861451
{
13871452
processMc<aod::hf_cand::VertexerType::KfParticle, true>(selectedD0CandidatesMlMcKF, mcParticles, tracks, collisions, mcCollisions, bcs);

0 commit comments

Comments
 (0)