Skip to content

Commit b2ce618

Browse files
Mingze129Mingze129
committed
Optimise and update the logic for MC analysis
1 parent 9ab1b1f commit b2ce618

1 file changed

Lines changed: 57 additions & 59 deletions

File tree

‎PWGHF/D2H/Tasks/taskCharmPolarisation.cxx‎

Lines changed: 57 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@
2222
#include "PWGHF/Core/HfHelper.h"
2323
#include "PWGHF/DataModel/CandidateReconstructionTables.h"
2424
#include "PWGHF/DataModel/CandidateSelectionTables.h"
25+
#include "PWGLF/DataModel/mcCentrality.h"
2526

2627
#include "Common/Core/EventPlaneHelper.h"
2728
#include "Common/Core/RecoDecay.h"
@@ -218,8 +219,6 @@ struct HfTaskCharmPolarisation {
218219

219220
using CollisionsWithMcLabels = soa::SmallGroups<soa::Join<aod::Collisions, aod::McCollisionLabels>>;
220221
using CollsWithQVecs = soa::Join<aod::Collisions, aod::EvSels, aod::QvectorFT0Cs, aod::QvectorFT0Ms, aod::QvectorFV0As, aod::CentFT0Ms, aod::CentFT0Cs>;
221-
using CollsWithQVecsWithMcLabels = soa::SmallGroups<soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::QvectorFT0Cs, aod::QvectorFT0Ms, aod::QvectorFV0As, aod::CentFT0Ms, aod::CentFT0Cs>>;
222-
using GenCollisWithQvecs = soa::Join<aod::McCollisions, aod::QvectorFT0Cs, aod::QvectorFT0Ms, aod::QvectorFV0As, aod::CentFT0Ms, aod::CentFT0Cs>;
223222
using TracksWithMcLabels = soa::Join<aod::Tracks, aod::TracksExtra, aod::McTrackLabels>;
224223
using TracksWithExtra = soa::Join<aod::Tracks, aod::TracksExtra>;
225224

@@ -247,6 +246,8 @@ struct HfTaskCharmPolarisation {
247246
Preslice<FilteredCandDstarWSelFlagAndMc> dstarWithMcPerCollision = aod::hf_cand::collisionId;
248247
Preslice<FilteredCandDstarWSelFlagAndMcAndMl> dstarWithMcAndMlPerCollision = aod::hf_cand::collisionId;
249248

249+
PresliceUnsorted<aod::McCollisionLabels> colPerMcCollision = aod::mcparticle::mcCollisionId;
250+
250251
Preslice<FilteredCandLcToPKPiWSelFlag> lcToPKPiPerCollision = aod::hf_cand::collisionId;
251252
Preslice<FilteredCandLcToPKPiWSelFlagAndMl> lcToPKPiWithMlPerCollision = aod::hf_cand::collisionId;
252253
Preslice<FilteredCandLcToPKPiWSelFlagAndMc> lcToPKPiWithMcPerCollision = aod::hf_cand::collisionId;
@@ -273,6 +274,7 @@ struct HfTaskCharmPolarisation {
273274
ConfigurableAxis configThnAxisNumItsClsMin{"configThnAxisNumItsClsMin", {4, 3.5f, 7.5f}, "min #it{N}_{cls ITS}"};
274275
ConfigurableAxis configThnAxisNumTpcClsMin{"configThnAxisNumTpcClsMin", {3, 79.5f, 140.5f}, "min #it{N}_{cls TPC}"};
275276
ConfigurableAxis configThnAxisCharge{"configThnAxisCharge", {2, -2.f, 2.f}, "electric charge"};
277+
ConfigurableAxis configThnAxisCentrality{"configThnAxisCentrality", {100, 0.f, 100.f}, "centrality (%)"};
276278

277279
HistogramRegistry registry{"registry", {}};
278280

@@ -309,12 +311,12 @@ struct HfTaskCharmPolarisation {
309311
}
310312
}
311313

312-
if (activatePartRecoDstar && !(doprocessDstarMc || doprocessDstarMcWithMl)) {
314+
if (activatePartRecoDstar && !(doprocessDstarMc || doprocessDstarMcWithMl || doprocessDstarMcInPbPb || doprocessDstarMcWithMlInPbPb)) {
313315
LOGP(fatal, "Check on partly reconstructed D* mesons only possible for processDstarMc and processDstarMcWithMl");
314316
}
315317

316318
// check bkg rotation for MC (not supported currently)
317-
if (nBkgRotations > 0 && (doprocessDstarMc || doprocessDstarMcWithMl || doprocessLcToPKPiMc || doprocessLcToPKPiMcWithMl || doprocessLcToPKPiBackgroundMcWithMl)) {
319+
if (nBkgRotations > 0 && (doprocessDstarMc || doprocessDstarMcWithMl || doprocessLcToPKPiMc || doprocessLcToPKPiMcWithMl || doprocessLcToPKPiBackgroundMcWithMl || doprocessDstarMcInPbPb || doprocessDstarMcWithMlInPbPb)) {
318320
LOGP(fatal, "No background rotation supported for MC.");
319321
}
320322

@@ -351,6 +353,7 @@ struct HfTaskCharmPolarisation {
351353
const AxisSpec thnAxisInvMass2PKLcMonitoring{lcPKPiChannels.configThnAxisInvMass2PKLcMonitoring, "#it{M}^{2}(pK) from #Lambda_{c}^{+} (GeV/#it{c}^{2})"};
352354
const AxisSpec thnAxisInvMass2PPiLcMonitoring{lcPKPiChannels.configThnAxisInvMass2PPiLcMonitoring, "#it{M}^{2}(p#pi) from #Lambda_{c}^{+} (GeV/#it{c}^{2})"};
353355
const AxisSpec thnAxisTHnAxisEulerPhi{configTHnAxisEulerPhi, "Euler polar angle #phi"};
356+
const AxisSpec thnAxisCentrality{configThnAxisCentrality, "centrality (%)"};
354357

355358
auto invMassBins = thnAxisInvMass.binEdges;
356359
minInvMass = invMassBins.front();
@@ -359,6 +362,9 @@ struct HfTaskCharmPolarisation {
359362
registry.add("hNumPvContributorsAll", "Number of PV contributors for all events ;num. PV contributors; counts", HistType::kTH1F, {thnAxisNumPvContributors});
360363
registry.add("hNumPvContributorsCand", "Number of PV contributors for events with candidates;num. PV contributors; counts", HistType::kTH1F, {thnAxisNumPvContributors});
361364
registry.add("hNumPvContributorsCandInMass", "Number of PV contributors for events with candidates in the signal region;num. PV contributors; counts", HistType::kTH1F, {thnAxisNumPvContributors});
365+
if (doprocessDstarInPbPb || doprocessDstarMcInPbPb || doprocessDstarWithMlInPbPb || doprocessDstarMcWithMlInPbPb) {
366+
registry.add("hCentrality", "Centrality distribution for D*+ candidates;centrality (%); counts", HistType::kTH1F, {thnAxisCentrality});
367+
}
362368

363369
if (activateTHnSparseCosThStarHelicity) {
364370
std::vector<AxisSpec> hHelicityaxes = {thnAxisInvMass, thnAxisPt, thnAxisNumPvContributors, thnAxisY};
@@ -639,46 +645,22 @@ struct HfTaskCharmPolarisation {
639645
}
640646
}
641647

642-
if (activateTHnSparseCosThStarEP && (!doprocessDstarInPbPb || !doprocessDstarMcInPbPb || !doprocessDstarWithMlInPbPb || !doprocessDstarMcWithMlInPbPb)) {
648+
if (activateTHnSparseCosThStarEP && !(doprocessDstarInPbPb || doprocessDstarWithMlInPbPb)) {
643649
LOGP(fatal, "THnSparse with cosThStar w.r.t. event plane axis is not supported for pp analysis, please check the configuration!");
644650
} else if (activateTHnSparseCosThStarEP) {
645651
std::vector<AxisSpec> hEPaxes = {thnAxisInvMass, thnAxisPt, thnAxisNumPvContributors, thnAxisY};
646-
if (doprocessDstar || doprocessDstarMc || doprocessDstarWithMl || doprocessDstarMcWithMl || doprocessDstarInPbPb || doprocessDstarMcInPbPb || doprocessDstarWithMlInPbPb || doprocessDstarMcWithMlInPbPb) {
647-
652+
if (doprocessDstarInPbPb || doprocessDstarWithMlInPbPb) {
648653
hEPaxes.insert(hEPaxes.end(), {thnAxisInvMassD0, thnAxisCosThetaStarEP});
649-
if (doprocessDstarWithMl || doprocessDstarMcWithMl || doprocessDstarWithMlInPbPb || doprocessDstarMcWithMlInPbPb) {
654+
if (doprocessDstarWithMlInPbPb) {
650655
hEPaxes.insert(hEPaxes.end(), {thnAxisMlBkg, thnAxisMlNonPrompt});
651656
}
652657
if (activateTrackingSys) {
653658
hEPaxes.insert(hEPaxes.end(), {thnAxisAbsEtaTrackMin, thnAxisNumItsClsMin, thnAxisNumTpcClsMin});
654659
}
655-
if (doprocessDstarMc || doprocessDstarMcWithMl || doprocessDstarMcInPbPb || doprocessDstarMcWithMlInPbPb) {
656-
std::vector<AxisSpec> hRecoPromptEPAxes(hEPaxes);
657-
hRecoPromptEPAxes.insert(hRecoPromptEPAxes.end(), {thnAxisDauToMuons});
658-
std::vector<AxisSpec> hRecoNonPromptEPAxes(hEPaxes);
659-
hRecoNonPromptEPAxes.insert(hRecoNonPromptEPAxes.end(), {thnAxisDauToMuons, thnAxisPtB});
660-
registry.add("hRecoPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for reconstructed prompt D*+ candidates", HistType::kTHnSparseF, hRecoPromptEPAxes);
661-
registry.add("hRecoNonPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for reconstructed non-prompt D*+ candidates", HistType::kTHnSparseF, hRecoNonPromptEPAxes);
662-
if (activatePartRecoDstar) {
663-
registry.add("hPartRecoPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for partially reconstructed prompt D*+ candidates", HistType::kTHnSparseF, hRecoPromptEPAxes);
664-
registry.add("hPartRecoNonPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for partially reconstructed non-prompt D*+ candidates", HistType::kTHnSparseF, hRecoNonPromptEPAxes);
665-
}
666-
} else {
667-
if (nBkgRotations > 0) {
668-
hEPaxes.push_back(thnAxisIsRotatedCandidate);
669-
}
670-
registry.add("hEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores", HistType::kTHnSparseF, hEPaxes);
671-
}
672-
}
673-
if (doprocessDstarMc || doprocessDstarMcWithMl || doprocessDstarMcInPbPb || doprocessDstarMcWithMlInPbPb || doprocessLcToPKPiMc || doprocessLcToPKPiMcWithMl || doprocessLcToPKPiBackgroundMcWithMl) {
674-
std::vector<AxisSpec> hgenPromptAxes = {thnAxisPt, thnAxisNumPvContributors, thnAxisY, thnAxisCosThetaStarEP, thnAxisDausAcc, thnAxisResoChannelLc, thnAxisCharge};
675-
std::vector<AxisSpec> hgenNonPromptAxes = {thnAxisPt, thnAxisNumPvContributors, thnAxisY, thnAxisCosThetaStarEP, thnAxisPtB, thnAxisDausAcc, thnAxisResoChannelLc, thnAxisCharge};
676-
registry.add("hGenPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes);
677-
registry.add("hGenNonPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes);
678-
if (activatePartRecoDstar) {
679-
registry.add("hPartRecoGenPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for partially reconstructed generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes);
680-
registry.add("hPartRecoGenNonPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for partially reconstructed generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes);
660+
if (nBkgRotations > 0) {
661+
hEPaxes.push_back(thnAxisIsRotatedCandidate);
681662
}
663+
registry.add("hEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores", HistType::kTHnSparseF, hEPaxes);
682664
}
683665
}
684666

@@ -1494,8 +1476,7 @@ struct HfTaskCharmPolarisation {
14941476

14951477
/// Get the Q vector
14961478
/// \param collision is the collision with the Q vector information
1497-
template <typename CollsWithQVecs>
1498-
std::vector<float> getQVec(CollsWithQVecs const& collision)
1479+
std::vector<float> getQVec(CollsWithQVecs::iterator const& collision)
14991480
{
15001481
float xQVec = -999.;
15011482
float yQVec = -999.;
@@ -2070,9 +2051,18 @@ struct HfTaskCharmPolarisation {
20702051
/// \param mcParticle is the Mc particle
20712052
/// \param mcParticles is the table of Mc particles
20722053
/// \param numPvContributors is the number of PV contributors in the associated reco collision
2073-
template <charm_polarisation::DecayChannel channel, typename Part, typename Particles>
2074-
void runMcGenPolarisationAnalysis(Part const& mcParticle, Particles const& mcParticles, int numPvContributors)
2054+
template <charm_polarisation::DecayChannel channel, bool withCent = false, typename Part, typename Particles, typename Cent = void>
2055+
void runMcGenPolarisationAnalysis(Part const& mcParticle, Particles const& mcParticles, int numPvContributors, Cent const* centrality = nullptr)
20752056
{
2057+
if constexpr (withCent) {
2058+
assert(qVecs && "Centrality analysis requested but Cent == nullptr");
2059+
}
2060+
if constexpr (withCent) {
2061+
if (*centrality < centralityMin || *centrality > centralityMax) {
2062+
return; // skip this collision if outside of the centrality range
2063+
}
2064+
}
2065+
20762066
int8_t origin{RecoDecay::OriginType::None};
20772067
std::vector<int> listDaughters{};
20782068
float massDau{0.f}, massCharmHad{0.f};
@@ -2316,14 +2306,16 @@ struct HfTaskCharmPolarisation {
23162306
PROCESS_SWITCH(HfTaskCharmPolarisation, processDstarMcWithMl, "Process Dstar candidates in MC with ML", false);
23172307

23182308
void processDstarInPbPb(CollsWithQVecs::iterator const& collision,
2319-
FilteredCandDstarWSelFlagAndMl const& dstarCandidates,
2309+
FilteredCandDstarWSelFlag const& dstarCandidates,
23202310
TracksWithExtra const& tracks)
23212311
{
23222312
float centrality = {-1.f};
23232313
centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator);
23242314
if (centrality < centralityMin || centrality > centralityMax) {
23252315
return; // skip this collision if outside of the centrality range
23262316
}
2317+
registry.fill(HIST("hCentrality"), centrality);
2318+
23272319
auto thisCollId = collision.globalIndex();
23282320
int numPvContributors = collision.numContrib();
23292321
auto groupedDstarCandidates = dstarCandidates.sliceBy(dstarPerCollision, thisCollId);
@@ -2350,6 +2342,8 @@ struct HfTaskCharmPolarisation {
23502342
if (centrality < centralityMin || centrality > centralityMax) {
23512343
return; // skip this collision if outside of the centrality range
23522344
}
2345+
registry.fill(HIST("hCentrality"), centrality);
2346+
23532347
auto thisCollId = collision.globalIndex();
23542348
int numPvContributors = collision.numContrib();
23552349
auto groupedDstarCandidates = dstarCandidates.sliceBy(dstarWithMlPerCollision, thisCollId);
@@ -2367,22 +2361,22 @@ struct HfTaskCharmPolarisation {
23672361
}
23682362
PROCESS_SWITCH(HfTaskCharmPolarisation, processDstarWithMlInPbPb, "Process Dstar candidates with ML in PbPb collisions", false);
23692363

2370-
void processDstarMcInPbPb(GenCollisWithQvecs::iterator const& collision,
2364+
void processDstarMcInPbPb(aod::McCollisions::iterator const&,
23712365
McParticlesDstarMatched const& mcParticles,
2372-
CollsWithQVecsWithMcLabels const& collisions, // this is grouped with SmallGroupsCollisionsWithMcLabels const& collisions,
2366+
CollisionsWithMcLabels const& collisions, // this is grouped with SmallGroupsCollisionsWithMcLabels const& collisions,
23732367
FilteredCandDstarWSelFlagAndMc const& dstarCandidates,
23742368
TracksWithExtra const& tracks)
23752369
{
23762370
float centrality = {-1.f};
2377-
centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator);
2378-
if (centrality < centralityMin || centrality > centralityMax) {
2379-
return; // skip this collision if outside of the centrality range
2380-
}
23812371
int numPvContributorsGen{0};
23822372

2383-
std::vector<float> qVecs = getQVec(collision);
2384-
23852373
for (const auto& collision : collisions) { // loop over reco collisions associated to this gen collision
2374+
centrality = o2::hf_centrality::getCentralityColl(collision);
2375+
if (centrality < centralityMin || centrality > centralityMax) {
2376+
return; // skip this collision if outside of the centrality range
2377+
}
2378+
registry.fill(HIST("hCentrality"), centrality);
2379+
23862380
auto thisCollId = collision.globalIndex();
23872381
int numPvContributors = collision.numContrib();
23882382
auto groupedDstarCandidates = dstarCandidates.sliceBy(dstarWithMcPerCollision, thisCollId);
@@ -2394,34 +2388,36 @@ struct HfTaskCharmPolarisation {
23942388

23952389
for (const auto& dstarCandidate : groupedDstarCandidates) {
23962390
nCands++;
2397-
if (runPolarisationAnalysis<charm_polarisation::DecayChannel::DstarToDzeroPi, false, true, false, true>(dstarCandidate, 0, numPvContributors, -1 /*MC particles*/, tracks, &qVecs)) {
2391+
if (runPolarisationAnalysis<charm_polarisation::DecayChannel::DstarToDzeroPi, false, true>(dstarCandidate, 0, numPvContributors, -1 /*MC particles*/, tracks)) {
23982392
nCandsInSignalRegion++;
23992393
}
24002394
}
24012395
fillMultHistos(numPvContributors, nCands, nCandsInSignalRegion);
24022396
}
24032397
for (const auto& mcParticle : mcParticles) {
2404-
runMcGenPolarisationAnalysis<charm_polarisation::DecayChannel::DstarToDzeroPi>(mcParticle, mcParticles, numPvContributorsGen);
2398+
const auto& recoCollsPerMcColl = collisions.sliceBy(colPerMcCollision, mcParticle.mcCollision().globalIndex());
2399+
float cent = o2::hf_centrality::getCentralityGenColl(recoCollsPerMcColl);
2400+
runMcGenPolarisationAnalysis<charm_polarisation::DecayChannel::DstarToDzeroPi, true>(mcParticle, mcParticles, numPvContributorsGen, &cent);
24052401
}
24062402
}
24072403
PROCESS_SWITCH(HfTaskCharmPolarisation, processDstarMcInPbPb, "Process Dstar candidates in PbPb MC without ML", false);
24082404

2409-
void processDstarMcWithMlInPbPb(GenCollisWithQvecs::iterator const& collision,
2405+
void processDstarMcWithMlInPbPb(aod::McCollisions::iterator const&,
24102406
McParticlesDstarMatched const& mcParticles,
2411-
CollsWithQVecsWithMcLabels const& collisions, // this is grouped with SmallGroupsCollisionsWithMcLabels const& collisions,
2407+
CollisionsWithMcLabels const& collisions, // this is grouped with SmallGroupsCollisionsWithMcLabels const& collisions,
24122408
FilteredCandDstarWSelFlagAndMcAndMl const& dstarCandidates,
24132409
TracksWithExtra const& tracks)
24142410
{
24152411
float centrality = {-1.f};
2416-
centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator);
2417-
if (centrality < centralityMin || centrality > centralityMax) {
2418-
return; // skip this collision if outside of the centrality range
2419-
}
24202412
int numPvContributorsGen{0};
24212413

2422-
std::vector<float> qVecs = getQVec(collision);
2423-
24242414
for (const auto& collision : collisions) { // loop over reco collisions associated to this gen collision
2415+
centrality = o2::hf_centrality::getCentralityColl(collision);
2416+
if (centrality < centralityMin || centrality > centralityMax) {
2417+
return; // skip this collision if outside of the centrality range
2418+
}
2419+
registry.fill(HIST("hCentrality"), centrality);
2420+
24252421
auto thisCollId = collision.globalIndex();
24262422
int numPvContributors = collision.numContrib();
24272423
auto groupedDstarCandidates = dstarCandidates.sliceBy(dstarWithMcAndMlPerCollision, thisCollId);
@@ -2433,14 +2429,16 @@ struct HfTaskCharmPolarisation {
24332429

24342430
for (const auto& dstarCandidate : groupedDstarCandidates) {
24352431
nCands++;
2436-
if (runPolarisationAnalysis<charm_polarisation::DecayChannel::DstarToDzeroPi, true, true, false, true>(dstarCandidate, 0, numPvContributors, -1 /*MC particles*/, tracks, &qVecs)) {
2432+
if (runPolarisationAnalysis<charm_polarisation::DecayChannel::DstarToDzeroPi, true, true>(dstarCandidate, 0, numPvContributors, -1 /*MC particles*/, tracks)) {
24372433
nCandsInSignalRegion++;
24382434
}
24392435
}
24402436
fillMultHistos(numPvContributors, nCands, nCandsInSignalRegion);
24412437
}
24422438
for (const auto& mcParticle : mcParticles) {
2443-
runMcGenPolarisationAnalysis<charm_polarisation::DecayChannel::DstarToDzeroPi>(mcParticle, mcParticles, numPvContributorsGen);
2439+
const auto& recoCollsPerMcColl = collisions.sliceBy(colPerMcCollision, mcParticle.mcCollision().globalIndex());
2440+
float cent = o2::hf_centrality::getCentralityGenColl(recoCollsPerMcColl);
2441+
runMcGenPolarisationAnalysis<charm_polarisation::DecayChannel::DstarToDzeroPi, true>(mcParticle, mcParticles, numPvContributorsGen, &cent);
24442442
}
24452443
}
24462444
PROCESS_SWITCH(HfTaskCharmPolarisation, processDstarMcWithMlInPbPb, "Process Dstar candidates in PbPb MC with ML", false);

0 commit comments

Comments
 (0)