Skip to content

Commit 908f3de

Browse files
MarcellocostiMarcellocosti
committed
[WIP] adding cent and occ for mc gen
1 parent 2e080cf commit 908f3de

File tree

2 files changed

+151
-72
lines changed

2 files changed

+151
-72
lines changed

‎PWGHF/D2H/Tasks/taskDplus.cxx‎

Lines changed: 149 additions & 70 deletions
Original file line numberDiff line numberDiff line change
@@ -72,9 +72,13 @@ struct HfTaskDplus {
7272
using CandDplusDataWithMl = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelDplusToPiKPi, aod::HfMlDplusToPiKPi>>;
7373
using CandDplusMcReco = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelDplusToPiKPi, aod::HfCand3ProngMcRec>>;
7474
using CandDplusMcRecoWithMl = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelDplusToPiKPi, aod::HfCand3ProngMcRec, aod::HfMlDplusToPiKPi>>;
75-
using McParticles = soa::Join<aod::McParticles, aod::HfCand3ProngMcGen>;
75+
using CandDplusMcGen = soa::Join<aod::McParticles, aod::HfCand3ProngMcGen>;
76+
7677
using CollisionsCent = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Ms, aod::CentFT0Cs>;
77-
using McCollisionsCent = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::CentFT0Ms, aod::CentFT0Cs>;
78+
using McRecoCollisionsCent = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::CentFT0Ms, aod::CentFT0Cs>;
79+
80+
PresliceUnsorted<aod::McCollisionLabels> recoColPerMcCollision = aod::mccollisionlabel::mcCollisionId;
81+
Preslice<aod::McParticles> mcParticlesPerMcCollision = aod::mcparticle::mcCollisionId;
7882

7983
Filter filterDplusFlag = (o2::aod::hf_track_index::hfflag & static_cast<uint8_t>(BIT(aod::hf_cand_3prong::DecayType::DplusToPiKPi))) != static_cast<uint8_t>(0);
8084

@@ -255,7 +259,7 @@ struct HfTaskDplus {
255259
float& centrality,
256260
float& occupancy)
257261
{
258-
LOG(info) << "info " << centrality << occupancy << ptbhad;
262+
LOG(info) << "info " << centrality << " " << occupancy << " " << ptbhad;
259263
std::vector<float> outputMl = {-999., -999., -999.};
260264
for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) {
261265
outputMl[iclass] = candidate.mlProbDplusToPiKPi()[classMl->at(iclass)];
@@ -397,13 +401,14 @@ struct HfTaskDplus {
397401
/// \param flagGenB transverse momentum of beauty mother for nonprompt candidates
398402
/// \param occupancy collision occupancy
399403
/// \param centrality collision centrality
400-
template <typename T2>
401-
void fillSparseMCGen(const T2& particle,
404+
template <typename T1>
405+
void fillSparseMCGen(const T1& particle,
402406
float& ptGenB,
403407
int& flagGenB,
404408
float& centrality,
405409
float& occupancy)
406410
{
411+
LOG(info) << "Filling sparse MC gen: " << " " << centrality << " " << occupancy;
407412
auto yGen = RecoDecay::y(particle.pVector(), o2::constants::physics::MassDPlus);
408413
if (particle.originMcGen() == RecoDecay::OriginType::Prompt) {
409414
if (storeCentrality && storeOccupancy) {
@@ -468,9 +473,8 @@ struct HfTaskDplus {
468473

469474
// Run analysis for the reconstructed Dplus candidates with MC matching
470475
/// \param candidates are reconstructed candidates
471-
/// \param mcParticles are particles with MC information
472-
template <bool fillMl, typename T1, typename T2>
473-
void runMCAnalysis(const T1& /*recoCandidates*/, const T2& mcParticles, McCollisionsCent const&)
476+
template <bool fillMl, typename T1>
477+
void runMCRecoAnalysis(const T1& /*recoCandidates*/, McRecoCollisionsCent const&)
474478
{
475479
float cent{-1};
476480
float occ{-1};
@@ -502,7 +506,7 @@ struct HfTaskDplus {
502506
flagBhad = getBHadMotherFlag(candidate.pdgBhadMotherPart());
503507

504508
if (storeCentrality || storeOccupancy) {
505-
auto collision = candidate.template collision_as<McCollisionsCent>();
509+
auto collision = candidate.template collision_as<McRecoCollisionsCent>();
506510
if (storeCentrality && centEstimator != CentralityEstimator::None) {
507511
cent = getCentrality(collision);
508512
}
@@ -522,44 +526,69 @@ struct HfTaskDplus {
522526
if ((yCandRecoMax >= 0. && std::abs(hfHelper.yDplus(candidate)) > yCandRecoMax)) {
523527
continue;
524528
}
525-
auto collision = candidate.template collision_as<McCollisionsCent>();
529+
auto collision = candidate.template collision_as<McRecoCollisionsCent>();
530+
if (storeCentrality && centEstimator != CentralityEstimator::None) {
531+
cent = getCentrality(collision);
532+
}
533+
if (storeOccupancy && occEstimator != OccupancyEstimator::None) {
534+
occ = getOccupancy(collision);
535+
}
526536
fillHistoMCRec<false>(candidate);
527537
fillSparseML<true, false>(candidate, ptBhad, flagBhad, cent, occ);
528538
}
529539
}
540+
}
541+
542+
// Run analysis for the reconstructed Dplus candidates with MC matching
543+
/// \param candidates are reconstructed candidates
544+
/// \param mcParticles are particles with MC information
545+
template <typename Cand>
546+
void runMCGenAnalysis(aod::McCollisions const& mcGenCollisions,
547+
McRecoCollisionsCent const& mcRecoCollisions,
548+
const Cand& mcGenParticles)
549+
{
530550
// MC gen.
531-
cent = -1.;
532-
occ = -1.;
551+
float cent{-1.};
552+
float occ{-1.};
533553
float ptGenB{-1.};
534554
int flagGenB{-1};
535-
for (const auto& particle : mcParticles) {
536-
auto yGen = RecoDecay::y(particle.pVector(), o2::constants::physics::MassDPlus);
537-
if ((yCandGenMax >= 0. && std::abs(yGen) > yCandGenMax) || (std::abs(particle.flagMcMatchGen()) != 1 << aod::hf_cand_3prong::DecayType::DplusToPiKPi)) {
538-
continue;
539-
}
540-
if (particle.originMcGen() == RecoDecay::OriginType::NonPrompt) {
541-
auto bHadMother = mcParticles.rawIteratorAt(particle.idxBhadMotherPart() - mcParticles.offset());
542-
flagGenB = getBHadMotherFlag(bHadMother.pdgCode());
543-
ptGenB = bHadMother.pt();
555+
556+
for (const auto& mcGenCollision : mcGenCollisions) {
557+
const auto recoCollsPerGenMcColl = mcRecoCollisions.sliceBy(recoColPerMcCollision, mcGenCollision.globalIndex());
558+
const auto mcParticlesPerGenMcColl = mcGenParticles.sliceBy(mcParticlesPerMcCollision, mcGenCollision.globalIndex());
559+
if (storeCentrality || storeOccupancy) {
560+
cent = -1.;
561+
occ = -1.;
562+
if (storeCentrality && centEstimator != CentralityEstimator::None) {
563+
cent = getMcGenCollCentrality(recoCollsPerGenMcColl);
564+
}
565+
if (storeOccupancy && occEstimator != OccupancyEstimator::None) {
566+
occ = getMcGenCollOccupancy(recoCollsPerGenMcColl);
567+
}
568+
569+
for (const auto& particle : mcParticlesPerGenMcColl) {
570+
ptGenB = -1;
571+
flagGenB = -1;
572+
auto yGen = RecoDecay::y(particle.pVector(), o2::constants::physics::MassDPlus);
573+
if ((yCandGenMax >= 0. && std::abs(yGen) > yCandGenMax) || (std::abs(particle.flagMcMatchGen()) != 1 << aod::hf_cand_3prong::DecayType::DplusToPiKPi)) {
574+
continue;
575+
}
576+
if (particle.originMcGen() == RecoDecay::OriginType::NonPrompt) {
577+
auto bHadMother = mcGenParticles.rawIteratorAt(particle.idxBhadMotherPart() - mcGenParticles.offset());
578+
flagGenB = getBHadMotherFlag(bHadMother.pdgCode());
579+
ptGenB = bHadMother.pt();
580+
}
581+
fillHistoMCGen(particle);
582+
fillSparseMCGen(particle, ptGenB, flagGenB, cent, occ);
583+
}
544584
}
545-
// if (storeCentrality || storeOccupancy) {
546-
// auto collision = particle.template mcCollision_as<o2::aod::McCollisions>();
547-
// if (storeCentrality && centEstimator != CentralityEstimator::None) {
548-
// cent = getCentrality(collision);
549-
// }
550-
// if (storeOccupancy && occEstimator != OccupancyEstimator::None) {
551-
// occ = getOccupancy(collision);
552-
// }
553-
// }
554-
fillHistoMCGen(particle);
555-
fillSparseMCGen(particle, ptGenB, flagGenB, cent, occ);
556585
}
557586
}
558587

559588
// /// Get the occupancy
560589
// /// \param collision is the collision with the occupancy information
561590
// template <typename Cand>
562-
// float getOccupancy(Cand const& candidate, CollisionsCent const&, McCollisionsCent const&)
591+
// float getOccupancy(Cand const& candidate, CollisionsCent const&, McRecoCollisionsCent const&)
563592
// {
564593
// float occupancy = -999.;
565594
// if constexpr (std::is_same_v<Cand, CandDplusData::iterator>) {
@@ -576,8 +605,8 @@ struct HfTaskDplus {
576605
// occupancy = collision.trackOccupancyInTimeRange();
577606
// break;
578607
// }
579-
// } else if constexpr (std::is_same_v<Cand, McParticles::iterator>) {
580-
// auto collision = candidate.template mcCollision_as<McCollisionsCent>();
608+
// } else if constexpr (std::is_same_v<Cand, CandDplusMcGen::iterator>) {
609+
// auto collision = candidate.template mcCollision_as<McRecoCollisionsCent>();
581610
// switch (occEstimator) {
582611
// case 1:
583612
// occupancy = collision.trackOccupancyInTimeRange();
@@ -591,7 +620,7 @@ struct HfTaskDplus {
591620
// break;
592621
// }
593622
// } else {
594-
// auto collision = candidate.template collision_as<McCollisionsCent>();
623+
// auto collision = candidate.template collision_as<McRecoCollisionsCent>();
595624
// switch (occEstimator) {
596625
// case 1:
597626
// occupancy = collision.trackOccupancyInTimeRange();
@@ -629,6 +658,75 @@ struct HfTaskDplus {
629658
return occupancy;
630659
}
631660

661+
/// \brief Function to get MC collision occupancy
662+
/// \param collSlice collection of reconstructed collisions
663+
/// \return collision occupancy
664+
template <typename CCs>
665+
int getMcGenCollOccupancy(CCs const& collSlice)
666+
{
667+
LOG(info) << "Getting MC gen collision occupancy";
668+
LOG(info) << "Associated reconstructed collisions: " << collSlice.size();
669+
float multiplicity{0.f};
670+
int occupancy = 0;
671+
for (const auto& collision : collSlice) {
672+
float collMult{0.f};
673+
collMult = collision.numContrib();
674+
// LOG(info) << "Slice collision multiplicity: " << collMult;
675+
676+
if (collMult > multiplicity) {
677+
occupancy = getOccupancy(collision);
678+
multiplicity = collMult;
679+
// LOG(info) << "Updating occupancy to: " << occupancy;
680+
// LOG(info) << "Updating multiplicity to: " << collMult;
681+
}
682+
} // end loop over collisions
683+
684+
// LOG(info) << " Occupancy --------->" << occupancy;
685+
return occupancy;
686+
}
687+
688+
/// Get the centrality
689+
/// \param collision is the collision with the centrality information
690+
template <typename Coll>
691+
float getCentrality(Coll const& collision)
692+
{
693+
float cent = -999.;
694+
switch (centEstimator) {
695+
case CentralityEstimator::FT0C:
696+
cent = collision.centFT0C();
697+
break;
698+
case CentralityEstimator::FT0M:
699+
cent = collision.centFT0M();
700+
break;
701+
default:
702+
LOG(warning) << "Centrality estimator not valid. Possible values are FT0C, FT0M. Fallback to FT0C";
703+
cent = collision.centFT0C();
704+
break;
705+
}
706+
return cent;
707+
}
708+
709+
/// Get the centrality
710+
/// \param collision is the collision with the centrality information
711+
float getMcGenCollCentrality(McRecoCollisionsCent const& collSlice) {
712+
LOG(info) << "Getting MC gen collision centrality";
713+
LOG(info) << "Associated reconstructed collisions: " << collSlice.size();
714+
float centrality{-1};
715+
float multiplicity{0.f};
716+
for (const auto& collision : collSlice) {
717+
float collMult = collision.numContrib();
718+
LOG(info) << "Slice collision multiplicity: " << collMult;
719+
if (collMult > multiplicity) {
720+
LOG(info) << "Updating centrality to: " << getCentrality(collision);
721+
LOG(info) << "Updating multiplicity to: " << collMult;
722+
centrality = getCentrality(collision);
723+
multiplicity = collMult;
724+
}
725+
}
726+
LOG(info) << " Centrality ---------> " << centrality;
727+
return centrality;
728+
}
729+
632730
// /// Get the centrality
633731
// /// \param collision is the collision with the centrality information
634732
// template <typename Cand>
@@ -649,8 +747,8 @@ struct HfTaskDplus {
649747
// cent = collision.centFT0C();
650748
// break;
651749
// }
652-
// } else if constexpr (std::is_same_v<Cand, McParticles::iterator>) {
653-
// auto collision = candidate.template mcCollision_as<McCollisionsCent>();
750+
// } else if constexpr (std::is_same_v<Cand, CandDplusMcGen::iterator>) {
751+
// auto collision = candidate.template mcCollision_as<McRecoCollisionsCent>();
654752
// switch (centEstimator) {
655753
// case CentralityEstimator::FT0C:
656754
// cent = collision.centFT0C();
@@ -664,7 +762,7 @@ struct HfTaskDplus {
664762
// break;
665763
// }
666764
// } else {
667-
// auto collision = candidate.template collision_as<McCollisionsCent>();
765+
// auto collision = candidate.template collision_as<McRecoCollisionsCent>();
668766
// switch (centEstimator) {
669767
// case CentralityEstimator::FT0C:
670768
// cent = collision.centFT0C();
@@ -681,27 +779,6 @@ struct HfTaskDplus {
681779
// return cent;
682780
// }
683781

684-
/// Get the centrality
685-
/// \param collision is the collision with the centrality information
686-
template <typename Coll>
687-
float getCentrality(Coll const& collision)
688-
{
689-
float cent = -999.;
690-
switch (centEstimator) {
691-
case CentralityEstimator::FT0C:
692-
cent = collision.centFT0C();
693-
break;
694-
case CentralityEstimator::FT0M:
695-
cent = collision.centFT0M();
696-
break;
697-
default:
698-
LOG(warning) << "Centrality estimator not valid. Possible values are FT0C, FT0M. Fallback to FT0C";
699-
cent = collision.centFT0C();
700-
break;
701-
}
702-
return cent;
703-
}
704-
705782
/// Convert the B hadron mother PDG for non prompt candidates to a flag
706783
/// \param pdg of the b hadron mother
707784
int getBHadMotherFlag(const int& flagBhad)
@@ -734,21 +811,23 @@ struct HfTaskDplus {
734811
}
735812
PROCESS_SWITCH(HfTaskDplus, processDataWithMl, "Process data with ML", false);
736813

737-
void processMc(CandDplusMcReco const& candidates,
738-
McParticles const& mcParticles,
739-
McCollisionsCent const& mcCollisions,
740-
aod::McCollisions const&)
814+
void processMc(CandDplusMcReco const& mcRecoParticles,
815+
CandDplusMcGen const& mcGenParticles,
816+
McRecoCollisionsCent const& mcRecoCollisions,
817+
aod::McCollisions const& mcGenCollisions)
741818
{
742-
runMCAnalysis<false>(candidates, mcParticles, mcCollisions);
819+
runMCRecoAnalysis<false>(mcRecoParticles, mcRecoCollisions);
820+
runMCGenAnalysis(mcGenCollisions, mcRecoCollisions, mcGenParticles);
743821
}
744822
PROCESS_SWITCH(HfTaskDplus, processMc, "Process MC w/o ML", false);
745823

746-
void processMcWithMl(CandDplusMcRecoWithMl const& candidates,
747-
McParticles const& mcParticles,
748-
McCollisionsCent const& mcCollisions,
749-
aod::McCollisions const&)
824+
void processMcWithMl(CandDplusMcRecoWithMl const& mcRecoParticles,
825+
CandDplusMcGen const& mcGenParticles,
826+
McRecoCollisionsCent const& mcRecoCollisions,
827+
aod::McCollisions const& mcGenCollisions)
750828
{
751-
runMCAnalysis<true>(candidates, mcParticles, mcCollisions);
829+
runMCRecoAnalysis<true>(mcRecoParticles, mcRecoCollisions);
830+
runMCGenAnalysis(mcGenCollisions, mcRecoCollisions, mcGenParticles);
752831
}
753832
PROCESS_SWITCH(HfTaskDplus, processMcWithMl, "Process MC with ML", false);
754833
};

‎PWGHF/Utils/utilsEvSelHf.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -320,9 +320,9 @@ struct HfEventSelection : o2::framework::ConfigurableGroup {
320320
hPosYAfterEvSel->Fill(collision.posY());
321321
hPosZAfterEvSel->Fill(posZ);
322322
hNumPvContributorsAfterSel->Fill(collision.numContrib());
323-
LOG(info) << "Fillng centrality: " << centrality;
323+
// LOG(info) << "Fillng centrality: " << centrality;
324324
hSelCollisionsCent->Fill(centrality);
325-
LOG(info) << "Fillng centrality/occupancy: " << centrality << "/" << occupancy;
325+
// LOG(info) << "Fillng centrality/occupancy: " << centrality << "/" << occupancy;
326326
hCollisionsCentOcc->Fill(centrality, occupancy);
327327
}
328328
};

0 commit comments

Comments
 (0)