Skip to content

Commit 7018e5b

Browse files
fgrosafgrosa
authored
PWGHF: Adapt process function for generated particles for cent cuts (#6996)
1 parent eb6cbbe commit 7018e5b

1 file changed

Lines changed: 133 additions & 40 deletions

File tree

‎PWGHF/Tasks/taskMcValidation.cxx‎

Lines changed: 133 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -22,11 +22,13 @@
2222
#include "Framework/AnalysisTask.h"
2323
#include "Framework/HistogramRegistry.h"
2424
#include "Framework/runDataProcessing.h"
25+
#include "Framework/RunningWorkflowInfo.h"
2526
#include "Framework/StaticFor.h"
2627

2728
#include "CCDB/BasicCCDBManager.h"
2829
#include "Common/DataModel/CollisionAssociationTables.h"
2930

31+
#include "PWGHF/Core/CentralityEstimation.h"
3032
#include "PWGHF/DataModel/CandidateReconstructionTables.h"
3133
#include "PWGHF/DataModel/CandidateSelectionTables.h"
3234
#include "PWGHF/Utils/utilsEvSelHf.h"
@@ -37,6 +39,7 @@ using namespace o2::aod;
3739
using namespace o2::framework;
3840
using namespace o2::framework::expressions;
3941
using namespace o2::hf_evsel;
42+
using namespace o2::hf_centrality;
4043

4144
namespace
4245
{
@@ -82,13 +85,26 @@ static constexpr std::string_view originNames[nOriginTypes] = {"Prompt", "NonPro
8285
/// - Number of candidates per collision
8386
/// - Momentum Conservation for these particles
8487
struct HfTaskMcValidationGen {
88+
89+
using BCsInfo = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels>;
90+
using CollisionsNoCents = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels>;
91+
using CollisionsFT0Cs = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::CentFT0Cs>;
92+
using CollisionsFT0Ms = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::CentFT0Ms>;
93+
PresliceUnsorted<CollisionsNoCents> colPerMcCollision = aod::mccollisionlabel::mcCollisionId;
94+
PresliceUnsorted<CollisionsFT0Cs> colPerMcCollisionFT0C = aod::mccollisionlabel::mcCollisionId;
95+
PresliceUnsorted<CollisionsFT0Ms> colPerMcCollisionFT0M = aod::mccollisionlabel::mcCollisionId;
96+
Preslice<aod::McParticles> mcParticlesPerMcCollision = aod::mcparticle::mcCollisionId;
97+
8598
Configurable<double> xVertexMin{"xVertexMin", -100., "min. x of generated primary vertex [cm]"};
8699
Configurable<double> xVertexMax{"xVertexMax", 100., "max. x of generated primary vertex [cm]"};
87100
Configurable<double> yVertexMin{"yVertexMin", -100., "min. y of generated primary vertex [cm]"};
88101
Configurable<double> yVertexMax{"yVertexMax", 100., "max. y of generated primary vertex [cm]"};
89102
Configurable<double> zVertexMin{"zVertexMin", -100., "min. z of generated primary vertex [cm]"};
90103
Configurable<double> zVertexMax{"zVertexMax", 100., "max. z of generated primary vertex [cm]"};
91104
Configurable<int> eventGeneratorType{"eventGeneratorType", -1, "If positive, enable event selection using subGeneratorId information. The value indicates which events to keep (0 = MB, 4 = charm triggered, 5 = beauty triggered)"};
105+
Configurable<bool> rejectParticlesFromBkgEvent{"rejectParticlesFromBkgEvent", true, "Reject particles"};
106+
107+
HfEventSelectionMc hfEvSelMc; // mc event selection and monitoring
92108

93109
AxisSpec axisNhadrons{10, -0.5, 9.5};
94110
AxisSpec axisNquarks{20, -0.5, 19.5};
@@ -126,8 +142,14 @@ struct HfTaskMcValidationGen {
126142
{"NonPromptCharmBaryons/hNonPromptBaryonsYDistr", "Y distribution vs non-prompt charm baryon; ; #it{y}^{gen}", {HistType::kTH2F, {axisBaryonSpecies, axisY}}},
127143
{"NonPromptCharmBaryons/hNonPromptBaryonsDecLenDistr", "Decay length distribution vs non-prompt charm baryon; ; decay length (#mum)", {HistType::kTH2F, {axisBaryonSpecies, axisDecLen}}}}};
128144

129-
void init(InitContext&)
145+
void init(InitContext& initContext)
130146
{
147+
148+
std::array<bool, 3> processes = {doprocessNoCentSel, doprocessCentFT0C, doprocessCentFT0M};
149+
if (std::accumulate(processes.begin(), processes.end(), 0) > 1) {
150+
LOGP(fatal, "At most one process function for generated particles can be enabled at a time.");
151+
}
152+
131153
// add per species histograms
132154
for (size_t iOrigin = 0; iOrigin < nOriginTypes; iOrigin++) {
133155
for (int iChannel = 0; iChannel < nMesonChannels; iChannel++) {
@@ -156,52 +178,59 @@ struct HfTaskMcValidationGen {
156178
registry.get<TH2>(HIST("NonPromptCharmBaryons/hNonPromptBaryonsYDistr"))->GetXaxis()->SetBinLabel(iBin, labels[iBin + nMesonChannels - 1].data());
157179
registry.get<TH2>(HIST("NonPromptCharmBaryons/hNonPromptBaryonsDecLenDistr"))->GetXaxis()->SetBinLabel(iBin, labels[iBin + nMesonChannels - 1].data());
158180
}
159-
}
160181

161-
/// Primary-vertex selection
162-
/// \param collision mcCollision table row
163-
template <typename Col>
164-
bool selectVertex(const Col& collision)
165-
{
166-
// x position
167-
if (collision.posX() < xVertexMin || collision.posX() > xVertexMax) {
168-
return false;
169-
}
170-
// y position
171-
if (collision.posY() < yVertexMin || collision.posY() > yVertexMax) {
172-
return false;
173-
}
174-
// z position
175-
if (collision.posZ() < zVertexMin || collision.posZ() > zVertexMax) {
176-
return false;
182+
// inspect for which particle species the candidates were created and which zPvPosMax cut was set for reconstructed
183+
const auto& workflows = initContext.services().get<RunningWorkflowInfo const>();
184+
for (const DeviceSpec& device : workflows.devices) {
185+
if (device.name.compare("hf-task-mc-validation-rec") == 0) {
186+
hfEvSelMc.configureFromDevice(device);
187+
break;
188+
}
177189
}
178-
return true;
190+
191+
hfEvSelMc.addHistograms(registry); // particles monitoring
179192
}
180193

181-
void process(aod::McCollision const& mcCollision,
182-
aod::McParticles const& mcParticles)
194+
template <o2::hf_centrality::CentralityEstimator centEstimator, typename GenColl, typename Particles, typename RecoColls>
195+
void runCheckGenParticles(GenColl const& mcCollision, Particles const& mcParticles, RecoColls const& recoCollisions, BCsInfo const&, std::array<int, nChannels>& counterPrompt, std::array<int, nChannels>& counterNonPrompt)
183196
{
184197
if (eventGeneratorType >= 0 && mcCollision.getSubGeneratorId() != eventGeneratorType) {
185198
return;
186199
}
187-
if (!selectVertex(mcCollision)) {
200+
201+
// Slice the collisions table to get the collision info for the current MC collision
202+
float centrality{-1.f};
203+
uint16_t rejectionMask{0};
204+
if constexpr (centEstimator == CentralityEstimator::FT0C) {
205+
rejectionMask = hfEvSelMc.getHfMcCollisionRejectionMask<BCsInfo, centEstimator>(mcCollision, recoCollisions, centrality);
206+
} else if constexpr (centEstimator == CentralityEstimator::FT0M) {
207+
rejectionMask = hfEvSelMc.getHfMcCollisionRejectionMask<BCsInfo, centEstimator>(mcCollision, recoCollisions, centrality);
208+
} else if constexpr (centEstimator == CentralityEstimator::None) {
209+
rejectionMask = hfEvSelMc.getHfMcCollisionRejectionMask<BCsInfo, centEstimator>(mcCollision, recoCollisions, centrality);
210+
}
211+
hfEvSelMc.fillHistograms(rejectionMask);
212+
if (rejectionMask != 0) {
188213
return;
189214
}
190215

191216
int cPerCollision = 0;
192217
int cBarPerCollision = 0;
193218
int bPerCollision = 0;
194219
int bBarPerCollision = 0;
195-
std::array<int, nChannels> counterPrompt{0}, counterNonPrompt{0};
196220

197221
for (const auto& particle : mcParticles) {
222+
223+
if (rejectParticlesFromBkgEvent && particle.fromBackgroundEvent()) {
224+
continue;
225+
}
226+
198227
if (!particle.has_mothers()) {
199228
continue;
200229
}
201230

202231
int particlePdgCode = particle.pdgCode();
203232
bool isDiffFromMothers = true;
204-
for (const auto& mother : particle.mothers_as<aod::McParticles>()) {
233+
for (const auto& mother : particle.template mothers_as<Particles>()) {
205234
if (particlePdgCode == mother.pdgCode()) {
206235
isDiffFromMothers = false;
207236
break;
@@ -324,7 +353,7 @@ struct HfTaskMcValidationGen {
324353
counterNonPrompt[iD]++;
325354
}
326355

327-
auto daughter0 = particle.daughters_as<aod::McParticles>().begin();
356+
auto daughter0 = particle.template daughters_as<Particles>().begin();
328357
double vertexDau[3] = {daughter0.vx(), daughter0.vy(), daughter0.vz()};
329358
double vertexPrimary[3] = {mcCollision.posX(), mcCollision.posY(), mcCollision.posZ()};
330359
auto decayLength = RecoDecay::distance(vertexPrimary, vertexDau);
@@ -364,17 +393,79 @@ struct HfTaskMcValidationGen {
364393
registry.fill(HIST("Quarks/hCountB"), bPerCollision);
365394
registry.fill(HIST("Quarks/hCountCbar"), cBarPerCollision);
366395
registry.fill(HIST("Quarks/hCountBbar"), bBarPerCollision);
367-
static_for<0, nMesonChannels - 1>([&](auto i) {
368-
constexpr int index = i.value;
369-
registry.fill(HIST("PromptCharmMesons/hCountPrompt") + HIST(particleNames[index]), counterPrompt[index]);
370-
registry.fill(HIST("NonPromptCharmMesons/hCountNonPrompt") + HIST(particleNames[index]), counterNonPrompt[index]);
371-
});
372-
static_for<nMesonChannels, nChannels - 1>([&](auto i) {
373-
constexpr int index = i.value;
374-
registry.fill(HIST("PromptCharmBaryons/hCountPrompt") + HIST(particleNames[index]), counterPrompt[index]);
375-
registry.fill(HIST("NonPromptCharmBaryons/hCountNonPrompt") + HIST(particleNames[index]), counterNonPrompt[index]);
376-
});
377-
};
396+
}
397+
398+
void processNoCentSel(aod::McCollisions const& mcCollisions,
399+
aod::McParticles const& mcParticles,
400+
CollisionsNoCents const& recoCollisions,
401+
BCsInfo const& bcInfo)
402+
{
403+
for (const auto& mcCollision : mcCollisions) {
404+
const auto recoCollsPerMcColl = recoCollisions.sliceBy(colPerMcCollision, mcCollision.globalIndex());
405+
const auto mcParticlesPerMcColl = mcParticles.sliceBy(mcParticlesPerMcCollision, mcCollision.globalIndex());
406+
std::array<int, nChannels> counterPrompt{0}, counterNonPrompt{0};
407+
runCheckGenParticles<o2::hf_centrality::CentralityEstimator::None>(mcCollision, mcParticlesPerMcColl, recoCollsPerMcColl, bcInfo, counterPrompt, counterNonPrompt);
408+
static_for<0, nMesonChannels - 1>([&](auto i) {
409+
constexpr int index = i.value;
410+
registry.fill(HIST("PromptCharmMesons/hCountPrompt") + HIST(particleNames[index]), counterPrompt[index]);
411+
registry.fill(HIST("NonPromptCharmMesons/hCountNonPrompt") + HIST(particleNames[index]), counterNonPrompt[index]);
412+
});
413+
static_for<nMesonChannels, nChannels - 1>([&](auto i) {
414+
constexpr int index = i.value;
415+
registry.fill(HIST("PromptCharmBaryons/hCountPrompt") + HIST(particleNames[index]), counterPrompt[index]);
416+
registry.fill(HIST("NonPromptCharmBaryons/hCountNonPrompt") + HIST(particleNames[index]), counterNonPrompt[index]);
417+
});
418+
}
419+
} // end processNoCentSel
420+
PROCESS_SWITCH(HfTaskMcValidationGen, processNoCentSel, "Process generated collisions information without centrality selection", true);
421+
422+
void processCentFT0C(aod::McCollisions const& mcCollisions,
423+
aod::McParticles const& mcParticles,
424+
CollisionsFT0Cs const& recoCollisions,
425+
BCsInfo const& bcInfo)
426+
{
427+
for (const auto& mcCollision : mcCollisions) {
428+
const auto recoCollsPerMcColl = recoCollisions.sliceBy(colPerMcCollisionFT0C, mcCollision.globalIndex());
429+
const auto mcParticlesPerMcColl = mcParticles.sliceBy(mcParticlesPerMcCollision, mcCollision.globalIndex());
430+
std::array<int, nChannels> counterPrompt{0}, counterNonPrompt{0};
431+
runCheckGenParticles<o2::hf_centrality::CentralityEstimator::FT0C>(mcCollision, mcParticlesPerMcColl, recoCollsPerMcColl, bcInfo, counterPrompt, counterNonPrompt);
432+
static_for<0, nMesonChannels - 1>([&](auto i) {
433+
constexpr int index = i.value;
434+
registry.fill(HIST("PromptCharmMesons/hCountPrompt") + HIST(particleNames[index]), counterPrompt[index]);
435+
registry.fill(HIST("NonPromptCharmMesons/hCountNonPrompt") + HIST(particleNames[index]), counterNonPrompt[index]);
436+
});
437+
static_for<nMesonChannels, nChannels - 1>([&](auto i) {
438+
constexpr int index = i.value;
439+
registry.fill(HIST("PromptCharmBaryons/hCountPrompt") + HIST(particleNames[index]), counterPrompt[index]);
440+
registry.fill(HIST("NonPromptCharmBaryons/hCountNonPrompt") + HIST(particleNames[index]), counterNonPrompt[index]);
441+
});
442+
}
443+
} // end processCentFT0C
444+
PROCESS_SWITCH(HfTaskMcValidationGen, processCentFT0C, "Process generated collisions information with centrality selection using FT0C", false);
445+
446+
void processCentFT0M(aod::McCollisions const& mcCollisions,
447+
aod::McParticles const& mcParticles,
448+
CollisionsFT0Ms const& recoCollisions,
449+
BCsInfo const& bcInfo)
450+
{
451+
for (const auto& mcCollision : mcCollisions) {
452+
const auto recoCollsPerMcColl = recoCollisions.sliceBy(colPerMcCollisionFT0M, mcCollision.globalIndex());
453+
const auto mcParticlesPerMcColl = mcParticles.sliceBy(mcParticlesPerMcCollision, mcCollision.globalIndex());
454+
std::array<int, nChannels> counterPrompt{0}, counterNonPrompt{0};
455+
runCheckGenParticles<o2::hf_centrality::CentralityEstimator::FT0M>(mcCollision, mcParticlesPerMcColl, recoCollsPerMcColl, bcInfo, counterPrompt, counterNonPrompt);
456+
static_for<0, nMesonChannels - 1>([&](auto i) {
457+
constexpr int index = i.value;
458+
registry.fill(HIST("PromptCharmMesons/hCountPrompt") + HIST(particleNames[index]), counterPrompt[index]);
459+
registry.fill(HIST("NonPromptCharmMesons/hCountNonPrompt") + HIST(particleNames[index]), counterNonPrompt[index]);
460+
});
461+
static_for<nMesonChannels, nChannels - 1>([&](auto i) {
462+
constexpr int index = i.value;
463+
registry.fill(HIST("PromptCharmBaryons/hCountPrompt") + HIST(particleNames[index]), counterPrompt[index]);
464+
registry.fill(HIST("NonPromptCharmBaryons/hCountNonPrompt") + HIST(particleNames[index]), counterNonPrompt[index]);
465+
});
466+
}
467+
} // end processCentFT0M
468+
PROCESS_SWITCH(HfTaskMcValidationGen, processCentFT0M, "Process generated collisions information with centrality selection using FT0M", false);
378469
};
379470

380471
/// Reconstruction Level Validation
@@ -397,6 +488,8 @@ struct HfTaskMcValidationRec {
397488
using HfCand2ProngWithMCRec = soa::Join<aod::HfCand2Prong, aod::HfCand2ProngMcRec>;
398489
using HfCand3ProngWithMCRec = soa::Join<aod::HfCand3Prong, aod::HfCand3ProngMcRec>;
399490
using CollisionsWithMCLabels = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels>;
491+
using CollisionsWithMCLabelsAndCentFT0C = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, CentFT0Cs>;
492+
using CollisionsWithMCLabelsAndCentFT0M = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, CentFT0Ms>;
400493
using TracksWithSel = soa::Join<aod::TracksWMc, aod::TracksExtra, aod::TrackSelection, aod::TrackCompColls>;
401494

402495
Partition<TracksWithSel> tracksFilteredGlobalTrackWoDCA = requireGlobalTrackWoDCAInFilter();
@@ -721,15 +814,15 @@ struct HfTaskMcValidationRec {
721814
} // end process
722815
PROCESS_SWITCH(HfTaskMcValidationRec, processColl, "Process collision information without centrality selection", true);
723816

724-
void processCollWithCentFTOC(soa::Join<CollisionsWithMCLabels, aod::CentFT0Cs>::iterator const& collision,
817+
void processCollWithCentFTOC(CollisionsWithMCLabelsAndCentFT0C::iterator const& collision,
725818
aod::McCollisions const& mcCollisions,
726819
aod::BCsWithTimestamps const& bcs)
727820
{
728821
checkCollisions<o2::hf_centrality::CentralityEstimator::FT0C>(collision, mcCollisions, bcs);
729822
} // end process
730823
PROCESS_SWITCH(HfTaskMcValidationRec, processCollWithCentFTOC, "Process collision information with centrality selection with FT0C", false);
731824

732-
void processCollWithCentFTOM(soa::Join<CollisionsWithMCLabels, aod::CentFT0Ms>::iterator const& collision,
825+
void processCollWithCentFTOM(CollisionsWithMCLabelsAndCentFT0M::iterator const& collision,
733826
aod::McCollisions const& mcCollisions,
734827
aod::BCsWithTimestamps const& bcs)
735828
{
@@ -747,7 +840,7 @@ struct HfTaskMcValidationRec {
747840
}
748841
PROCESS_SWITCH(HfTaskMcValidationRec, processCollAssoc, "Process collision-association information, requires extra table from TrackToCollisionAssociation task (fillTableOfCollIdsPerTrack=true)", false);
749842

750-
void processCollAssocWithCentFTOC(soa::Join<CollisionsWithMCLabels, aod::CentFT0Cs> const& collisions,
843+
void processCollAssocWithCentFTOC(CollisionsWithMCLabelsAndCentFT0C const& collisions,
751844
TracksWithSel const& tracks,
752845
aod::McParticles const& mcParticles,
753846
aod::McCollisions const& mcCollisions,
@@ -757,7 +850,7 @@ struct HfTaskMcValidationRec {
757850
}
758851
PROCESS_SWITCH(HfTaskMcValidationRec, processCollAssocWithCentFTOC, "Process collision-association information with centrality selection with FT0C, requires extra table from TrackToCollisionAssociation task (fillTableOfCollIdsPerTrack=true)", false);
759852

760-
void processCollAssocWithCentFTOM(soa::Join<CollisionsWithMCLabels, aod::CentFT0Ms> const& collisions,
853+
void processCollAssocWithCentFTOM(CollisionsWithMCLabelsAndCentFT0M const& collisions,
761854
TracksWithSel const& tracks,
762855
aod::McParticles const& mcParticles,
763856
aod::McCollisions const& mcCollisions,

0 commit comments

Comments
 (0)