Skip to content

Commit 1b8c02e

Browse files
zhangbiao-phyzhangbiao-phy
authored
[PWGHF] Add optional event-status filling to taskD0 (#17820)
1 parent e46fd70 commit 1b8c02e

1 file changed

Lines changed: 82 additions & 18 deletions

File tree

‎PWGHF/D2H/Tasks/taskD0.cxx‎

Lines changed: 82 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,7 @@ 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)"};
95104
Configurable<int> occEstimator{"occEstimator", 0, "Occupancy estimation (None: 0, ITS: 1, FT0C: 2)"};
96105
Configurable<bool> storeCentrality{"storeCentrality", false, "Flag to store centrality information"};
97106
Configurable<bool> storeOccupancyAndIR{"storeOccupancyAndIR", false, "Flag to store occupancy information and interaction rate"};
@@ -125,6 +134,7 @@ struct HfTaskD0 {
125134
using CollisionsCent = soa::Join<aod::Collisions, aod::EvSels, aod::PVMults, aod::CentFT0Ms, aod::CentFT0Cs>;
126135
using CollisionsWithMcLabels = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::PVMults>;
127136
using CollisionsWithMcLabelsCent = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::PVMults, aod::CentFT0Ms, aod::CentFT0Cs>;
137+
using McCollisionsWithCentrality = soa::Join<aod::McCollisions, aod::McCentFT0Cs, aod::McCentFT0Ms>;
128138
using TracksSelQuality = soa::Join<aod::TracksExtra, aod::TracksWMc>;
129139
using TracksWPid = soa::Join<o2::aod::FullTracks, aod::TracksDCA, o2::aod::TrackSelection, aod::TracksPidPi, aod::PidTpcTofFullPi, aod::TracksPidKa, aod::PidTpcTofFullKa, aod::TracksPidPr, aod::PidTpcTofFullPr>;
130140
// using TracksWithExtra = o2::soa::Join<o2::aod::FullTracks, o2::aod::TrackSelection>;
@@ -331,13 +341,22 @@ struct HfTaskD0 {
331341
std::vector<AxisSpec> axesAcc = {thnAxisGenPtD, thnAxisGenPtB, thnAxisY, thnAxisOrigin, thnAxisNumPvContr};
332342

333343
if (storeCentrality) {
334-
axesAcc.push_back(thnAxisCent);
344+
axesAcc.emplace_back(thnConfigAxisCent, fillEventStatus ? "Generated centrality (%)" : "Centrality");
335345
}
336346
// interaction rate only store in Data and MC Reco. Level
337347
if (storeOccupancyAndIR) {
338348
axesAcc.push_back(thnAxisOccupancy);
339349
}
340350

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

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

489512
ccdb->setURL(ccdbUrl);
490513
ccdb->setCaching(true);
@@ -970,7 +993,7 @@ struct HfTaskD0 {
970993
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
971994
TracksSelQuality const&,
972995
CollType const& collisions,
973-
aod::McCollisions const&,
996+
McCollisionsWithCentrality const&,
974997
BCsType const&)
975998
{
976999
// MC rec.
@@ -1255,7 +1278,31 @@ struct HfTaskD0 {
12551278
}
12561279
}
12571280
}
1258-
// MC gen.
1281+
// Classify reco associations once per collision, not once per generated particle.
1282+
// The task's hfEvSel.* configuration must match the candidate creator's cuts.
1283+
std::unordered_map<int64_t, int> recoEventStatus;
1284+
if (fillEventStatus) {
1285+
for (const auto& collision : collisions) {
1286+
if (!collision.has_mcCollision()) {
1287+
continue;
1288+
}
1289+
float unusedCentrality{-1.f};
1290+
const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask<true, CentralityEstimator::None, BCsType>(collision, unusedCentrality, ccdb, registry);
1291+
bool selected = rejectionMask == 0;
1292+
if constexpr (std::is_same_v<CollType, CollisionsWithMcLabelsCent>) {
1293+
if (centEstimator != CentralityEstimator::None) {
1294+
const auto recoCentrality = getCentralityColl(collision, centEstimator);
1295+
selected = selected && recoCentrality >= hfEvSel.centralityMin && recoCentrality <= hfEvSel.centralityMax;
1296+
}
1297+
}
1298+
auto& status = recoEventStatus[collision.mcCollisionId()];
1299+
const int currentStatus = selected ? SelectedRecoCollision : NoSelectedRecoCollision;
1300+
status = std::max(status, currentStatus);
1301+
}
1302+
}
1303+
1304+
// MC gen. Each particle is filled once, even for split or unreconstructed events.
1305+
// Upstream event rejection can clear flagMcMatchGen: this is not an unfiltered truth sample.
12591306
for (const auto& particle : mcParticles) {
12601307
if (std::abs(particle.flagMcMatchGen()) == o2::hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK) {
12611308
if (yCandGenMax >= 0. && std::abs(RecoDecay::y(particle.pVector(), o2::constants::physics::MassD0)) > yCandGenMax) {
@@ -1267,6 +1314,8 @@ struct HfTaskD0 {
12671314
registry.fill(HIST("hPtGen"), ptGen);
12681315
registry.fill(HIST("hPtVsYGen"), ptGen, yGen);
12691316

1317+
int eventStatus = NoRecoCollision;
1318+
12701319
unsigned maxNumContrib = 0;
12711320
float cent{-1.f};
12721321
float occ{-1.f};
@@ -1288,32 +1337,47 @@ struct HfTaskD0 {
12881337
}
12891338
}
12901339

1340+
if (fillEventStatus) {
1341+
const auto mcCollision = particle.template mcCollision_as<McCollisionsWithCentrality>();
1342+
const auto eventEntry = recoEventStatus.find(mcCollision.globalIndex());
1343+
eventStatus = eventEntry == recoEventStatus.end() ? NoRecoCollision : eventEntry->second;
1344+
// Replace the existing centrality coordinate, including zero-reco events.
1345+
cent = centEstimator == CentralityEstimator::FT0C ? mcCollision.centFT0C() : mcCollision.centFT0M();
1346+
}
1347+
const auto fillGeneratedSparse = [&](auto... coordinates) {
1348+
if (fillEventStatus) {
1349+
registry.fill(HIST("hSparseAcc"), coordinates..., eventStatus);
1350+
} else {
1351+
registry.fill(HIST("hSparseAcc"), coordinates...);
1352+
}
1353+
};
1354+
12911355
if (particle.originMcGen() == RecoDecay::OriginType::Prompt) {
12921356
registry.fill(HIST("hPtGenPrompt"), ptGen);
12931357
registry.fill(HIST("hYGenPrompt"), yGen);
12941358
registry.fill(HIST("hPtVsYGenPrompt"), ptGen, yGen);
12951359
if (storeCentrality && storeOccupancyAndIR) {
1296-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 1, maxNumContrib, cent, occ);
1360+
fillGeneratedSparse(ptGen, ptGenB, yGen, 1, maxNumContrib, cent, occ);
12971361
} else if (storeCentrality && !storeOccupancyAndIR) {
1298-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 1, maxNumContrib, cent);
1362+
fillGeneratedSparse(ptGen, ptGenB, yGen, 1, maxNumContrib, cent);
12991363
} else if (!storeCentrality && storeOccupancyAndIR) {
1300-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 1, maxNumContrib, occ);
1364+
fillGeneratedSparse(ptGen, ptGenB, yGen, 1, maxNumContrib, occ);
13011365
} else {
1302-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 1, maxNumContrib);
1366+
fillGeneratedSparse(ptGen, ptGenB, yGen, 1, maxNumContrib);
13031367
}
13041368
} else {
13051369
ptGenB = mcParticles.rawIteratorAt(particle.idxBhadMotherPart()).pt();
13061370
registry.fill(HIST("hPtGenNonPrompt"), ptGen);
13071371
registry.fill(HIST("hYGenNonPrompt"), yGen);
13081372
registry.fill(HIST("hPtVsYGenNonPrompt"), ptGen, yGen);
13091373
if (storeCentrality && storeOccupancyAndIR) {
1310-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 2, maxNumContrib, cent, occ);
1374+
fillGeneratedSparse(ptGen, ptGenB, yGen, 2, maxNumContrib, cent, occ);
13111375
} else if (storeCentrality && !storeOccupancyAndIR) {
1312-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 2, maxNumContrib, cent);
1376+
fillGeneratedSparse(ptGen, ptGenB, yGen, 2, maxNumContrib, cent);
13131377
} else if (!storeCentrality && storeOccupancyAndIR) {
1314-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 2, maxNumContrib, occ);
1378+
fillGeneratedSparse(ptGen, ptGenB, yGen, 2, maxNumContrib, occ);
13151379
} else {
1316-
registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 2, maxNumContrib);
1380+
fillGeneratedSparse(ptGen, ptGenB, yGen, 2, maxNumContrib);
13171381
}
13181382
}
13191383
registry.fill(HIST("hEtaGen"), particle.eta());
@@ -1325,7 +1389,7 @@ struct HfTaskD0 {
13251389
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13261390
TracksSelQuality const& tracks,
13271391
CollisionsWithMcLabels const& collisions,
1328-
aod::McCollisions const& mcCollisions,
1392+
McCollisionsWithCentrality const& mcCollisions,
13291393
aod::BcFullInfos const& bcs)
13301394
{
13311395
processMc<aod::hf_cand::VertexerType::DCAFitter, false>(selectedD0CandidatesMc, mcParticles, tracks, collisions, mcCollisions, bcs);
@@ -1336,7 +1400,7 @@ struct HfTaskD0 {
13361400
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13371401
TracksSelQuality const& tracks,
13381402
CollisionsWithMcLabelsCent const& collisions,
1339-
aod::McCollisions const& mcCollisions,
1403+
McCollisionsWithCentrality const& mcCollisions,
13401404
aod::BcFullInfos const& bcs)
13411405
{
13421406
processMc<aod::hf_cand::VertexerType::DCAFitter, false>(selectedD0CandidatesMc, mcParticles, tracks, collisions, mcCollisions, bcs);
@@ -1347,7 +1411,7 @@ struct HfTaskD0 {
13471411
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13481412
TracksSelQuality const& tracks,
13491413
CollisionsWithMcLabels const& collisions,
1350-
aod::McCollisions const& mcCollisions,
1414+
McCollisionsWithCentrality const& mcCollisions,
13511415
aod::BcFullInfos const& bcs)
13521416
{
13531417
processMc<aod::hf_cand::VertexerType::KfParticle, false>(selectedD0CandidatesMcKF, mcParticles, tracks, collisions, mcCollisions, bcs);
@@ -1359,7 +1423,7 @@ struct HfTaskD0 {
13591423
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13601424
TracksSelQuality const& tracks,
13611425
CollisionsWithMcLabels const& collisions,
1362-
aod::McCollisions const& mcCollisions,
1426+
McCollisionsWithCentrality const& mcCollisions,
13631427
aod::BcFullInfos const& bcs)
13641428
{
13651429
processMc<aod::hf_cand::VertexerType::DCAFitter, true>(selectedD0CandidatesMlMc, mcParticles, tracks, collisions, mcCollisions, bcs);
@@ -1370,7 +1434,7 @@ struct HfTaskD0 {
13701434
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13711435
TracksSelQuality const& tracks,
13721436
CollisionsWithMcLabelsCent const& collisions,
1373-
aod::McCollisions const& mcCollisions,
1437+
McCollisionsWithCentrality const& mcCollisions,
13741438
aod::BcFullInfos const& bcs)
13751439
{
13761440
processMc<aod::hf_cand::VertexerType::DCAFitter, true>(selectedD0CandidatesMlMc, mcParticles, tracks, collisions, mcCollisions, bcs);
@@ -1381,7 +1445,7 @@ struct HfTaskD0 {
13811445
soa::Join<aod::McParticles, aod::HfCand2ProngMcGen> const& mcParticles,
13821446
TracksSelQuality const& tracks,
13831447
CollisionsWithMcLabels const& collisions,
1384-
aod::McCollisions const& mcCollisions,
1448+
McCollisionsWithCentrality const& mcCollisions,
13851449
aod::BcFullInfos const& bcs)
13861450
{
13871451
processMc<aod::hf_cand::VertexerType::KfParticle, true>(selectedD0CandidatesMlMcKF, mcParticles, tracks, collisions, mcCollisions, bcs);

0 commit comments

Comments
 (0)