Skip to content

Commit 163770b

Browse files
BanajitBarmanBanajitBarman
committed
For centrality calibration in MC
1 parent 1baf73c commit 163770b

1 file changed

Lines changed: 92 additions & 67 deletions

File tree

‎PWGLF/Tasks/Nuspex/spectraTOF.cxx‎

Lines changed: 92 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -632,6 +632,8 @@ struct tofSpectra {
632632
hh->GetXaxis()->SetBinLabel(4, "INEL>1");
633633
hh->GetXaxis()->SetBinLabel(5, "hasParticleInFT0C && hasParticleInFT0A");
634634
histos.add("MC/MultiplicityRecoEv", "MC multiplicity", kTH1D, {multAxis});
635+
histos.add("MC/CentFT0CVsTrueMult", "MC Centrality vs True Multiplicity", kTH2D, {multAxis, {150, 0, 150, "Track multiplicity"}});
636+
635637
histos.add("MC/Multiplicity", "MC multiplicity", kTH1D, {multAxis});
636638
histos.add("MC/MultiplicityMCINELgt0", "MC multiplicity", kTH1D, {multAxis});
637639
histos.add("MC/MultiplicityMCINELgt1", "MC multiplicity", kTH1D, {multAxis});
@@ -645,7 +647,7 @@ struct tofSpectra {
645647
}
646648
}
647649

648-
hMultiplicityvsPercentile = histos.add<TH2>("Mult/vsPercentile", "Multiplicity vs percentile", HistType::kTH2D, {{150, 0, 150}, {100, 0, 100, "Track multiplicity"}});
650+
hMultiplicityvsPercentile = histos.add<TH2>("Mult/vsPercentile", "Multiplicity vs percentile", HistType::kTH2D, {{multAxis}, {150, 0, 150, "Track multiplicity"}});
649651

650652
for (int i = 0; i < NpCharge; i++) {
651653
switch (i) {
@@ -2487,26 +2489,27 @@ struct tofSpectra {
24872489
}
24882490
}
24892491

2492+
// ==============================================================================
2493+
// HELPER 1: For Reconstructed Events
2494+
// Accepts 'eventMultiplicity' (which will be the actual Reco Centrality)
2495+
// ==============================================================================
24902496
template <std::size_t i, typename ParticleType>
2491-
void fillParticleHistogramsMCRecoEvs(ParticleType const& mcParticle, RecoMCCollisions::iterator const& collision)
2497+
void fillParticleHistogramsMCRecoEvs(ParticleType const& mcParticle, RecoMCCollisions::iterator const& collision, float eventMultiplicity)
24922498
{
2493-
if (!isParticleEnabled<i>()) { // Check if the particle is enabled
2499+
if (!isParticleEnabled<i>()) {
24942500
return;
24952501
}
24962502

24972503
if (mcParticle.pdgCode() != PDGs[i]) {
24982504
return;
24992505
}
25002506

2501-
const auto& mcCollision = collision.mcCollision_as<GenMCCollisions>();
2502-
const float multiplicity = getMultiplicityMC(mcCollision);
2503-
25042507
if (mcParticle.isPhysicalPrimary()) {
25052508
if (isEventSelected<false, false>(collision)) {
25062509
if (includeCentralityMC) {
2507-
histos.fill(HIST(hpt_den_prm_goodev[i]), mcParticle.pt(), multiplicity, mcParticle.eta());
2510+
histos.fill(HIST(hpt_den_prm_goodev[i]), mcParticle.pt(), eventMultiplicity, mcParticle.eta());
25082511
} else {
2509-
histos.fill(HIST(hpt_den_prm_goodev[i]), mcParticle.pt(), multiplicity);
2512+
histos.fill(HIST(hpt_den_prm_goodev[i]), mcParticle.pt(), eventMultiplicity);
25102513
}
25112514
} else {
25122515
bool isSelected = collision.sel8();
@@ -2529,40 +2532,43 @@ struct tofSpectra {
25292532
}
25302533
if (isSelected) {
25312534
if (includeCentralityMC) {
2532-
histos.fill(HIST(hpt_den_prm_evsel[i]), mcParticle.pt(), multiplicity, mcParticle.eta());
2535+
histos.fill(HIST(hpt_den_prm_evsel[i]), mcParticle.pt(), eventMultiplicity, mcParticle.eta());
25332536
} else {
2534-
histos.fill(HIST(hpt_den_prm_evsel[i]), mcParticle.pt(), multiplicity);
2537+
histos.fill(HIST(hpt_den_prm_evsel[i]), mcParticle.pt(), eventMultiplicity);
25352538
}
25362539
}
25372540
} else {
25382541
if (includeCentralityMC) {
2539-
histos.fill(HIST(hpt_den_prm_recoev[i]), mcParticle.pt(), multiplicity, mcParticle.eta());
2542+
histos.fill(HIST(hpt_den_prm_recoev[i]), mcParticle.pt(), eventMultiplicity, mcParticle.eta());
25402543
} else {
2541-
histos.fill(HIST(hpt_den_prm_recoev[i]), mcParticle.pt(), multiplicity);
2544+
histos.fill(HIST(hpt_den_prm_recoev[i]), mcParticle.pt(), eventMultiplicity);
25422545
}
25432546
}
25442547
}
25452548
}
25462549
}
25472550

2551+
// ==============================================================================
2552+
// HELPER 2: For All Generated Events (Signal Loss Denominator)
2553+
// Accepts 'eventMultiplicity' (which will be the True Multiplicity |eta| < 0.5)
2554+
// ==============================================================================
25482555
template <std::size_t i, typename ParticleType>
2549-
void fillParticleHistogramsMCGenEvs(ParticleType const& mcParticle, GenMCCollisions::iterator const& mcCollision)
2556+
void fillParticleHistogramsMCGenEvs(ParticleType const& mcParticle, float eventMultiplicity)
25502557
{
2551-
2552-
if (!isParticleEnabled<i>()) { // Check if the particle is enabled
2558+
if (!isParticleEnabled<i>()) {
25532559
return;
25542560
}
25552561

25562562
if (mcParticle.pdgCode() != PDGs[i]) {
25572563
return;
25582564
}
25592565

2560-
const float multiplicity = getMultiplicityMC(mcCollision);
25612566
if (mcParticle.isPhysicalPrimary()) {
2562-
if (std::abs(mcCollision.posZ()) < evselOptions.cfgCutVertex) {
2563-
histos.fill(HIST(hpt_den_prm_mcgoodev[i]), mcParticle.pt(), multiplicity);
2567+
// Ensure this matches the array name for your generated denominator histograms!
2568+
if (includeCentralityMC) {
2569+
histos.fill(HIST(hpt_den_prm_mcgoodev[i]), mcParticle.pt(), eventMultiplicity);
25642570
} else {
2565-
histos.fill(HIST(hpt_den_prm_mcbadev[i]), mcParticle.pt(), multiplicity);
2571+
histos.fill(HIST(hpt_den_prm_mcbadev[i]), mcParticle.pt(), eventMultiplicity);
25662572
}
25672573
}
25682574
}
@@ -2584,6 +2590,9 @@ struct tofSpectra {
25842590
histos.fill(HIST("MC/GenRecoCollisions"), 1.f, mcCollisions.size());
25852591
histos.fill(HIST("MC/GenRecoCollisions"), 2.f, collisions.size());
25862592

2593+
// ==============================================================================
2594+
// 1. RECONSTRUCTED TRACKS LOOP
2595+
// ==============================================================================
25872596
for (const auto& track : tracks) {
25882597
if (!track.has_collision()) {
25892598
if (track.sign() > 0) {
@@ -2613,107 +2622,123 @@ struct tofSpectra {
26132622
fillTrackHistogramsMC<i>(track, mcParticle, track.collision_as<RecoMCCollisions>(), mcParticles);
26142623
});
26152624
}
2625+
2626+
// (Keeping your original includeCentralityMC legacy block intact to avoid breaking other parts of your code)
26162627
if (includeCentralityMC) {
26172628
for (const auto& collision : collisions) {
2618-
if (!collision.has_mcCollision()) {
2619-
continue;
2620-
}
2629+
if (!collision.has_mcCollision()) continue;
26212630
const auto& mcCollision = collision.mcCollision_as<GenMCCollisions>();
26222631
const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache);
26232632
const float multiplicity = getMultiplicity(collision);
26242633
for (const auto& mcParticle : particlesInCollision) {
2625-
2626-
if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) {
2627-
continue;
2628-
}
2629-
static_for<0, 17>([&](auto i) {
2630-
fillParticleHistogramsMC<i>(multiplicity, mcParticle);
2631-
});
2634+
if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) continue;
2635+
static_for<0, 17>([&](auto i) { fillParticleHistogramsMC<i>(multiplicity, mcParticle); });
26322636
}
26332637
}
26342638
} else {
26352639
for (const auto& mcParticle : mcParticles) {
2636-
if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) {
2637-
continue;
2638-
}
2639-
2640+
if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) continue;
26402641
const auto& mcCollision = mcParticle.mcCollision_as<GenMCCollisions>();
26412642
const float multiplicity = getMultiplicityMC(mcCollision);
2642-
2643-
static_for<0, 17>([&](auto i) {
2644-
fillParticleHistogramsMC<i>(multiplicity, mcParticle);
2645-
});
2643+
static_for<0, 17>([&](auto i) { fillParticleHistogramsMC<i>(multiplicity, mcParticle); });
26462644
}
26472645
}
2648-
// Loop on reconstructed collisions
2646+
2647+
// ==============================================================================
2648+
// 2. RECONSTRUCTED COLLISIONS LOOP (Builds the 2D Calibration Map)
2649+
// ==============================================================================
26492650
for (const auto& collision : collisions) {
26502651
if (!collision.has_mcCollision()) {
26512652
continue;
26522653
}
26532654
const auto& mcCollision = collision.mcCollision_as<GenMCCollisions>();
26542655
const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache);
2656+
26552657
if (evselOptions.cfgINELCut.value == 1) {
2656-
if (!o2::pwglf::isINELgt0mc(particlesInCollision, pdgDB)) {
2657-
continue;
2658-
}
2658+
if (!o2::pwglf::isINELgt0mc(particlesInCollision, pdgDB)) continue;
26592659
}
26602660
if (evselOptions.cfgINELCut.value == 2) {
2661-
if (!o2::pwglf::isINELgt1mc(particlesInCollision, pdgDB)) {
2662-
continue;
2661+
if (!o2::pwglf::isINELgt1mc(particlesInCollision, pdgDB)) continue;
2662+
}
2663+
2664+
// -- Get Actual Reconstructed Centrality --
2665+
float recoCentrality = collision.centFT0C();
2666+
2667+
// -- Calculate True Multiplicity in |eta| < 0.5 --
2668+
int trueMultEta05 = 0;
2669+
for (const auto& mcParticle : particlesInCollision) {
2670+
if (!mcParticle.isPhysicalPrimary()) continue;
2671+
auto* pdgParticle = pdgDB->GetParticle(mcParticle.pdgCode());
2672+
if (pdgParticle != nullptr && std::abs(pdgParticle->Charge()) > 0.01) {
2673+
if (std::abs(mcParticle.eta()) < 0.5f) {
2674+
trueMultEta05++;
2675+
}
26632676
}
26642677
}
26652678

26662679
if (isEventSelected<false, false>(collision)) {
2667-
histos.fill(HIST("MC/MultiplicityRecoEv"), getMultiplicityMC(mcCollision));
2680+
histos.fill(HIST("MC/MultiplicityRecoEv"), recoCentrality);
2681+
// -- Fill the 2D map to extract the 1D calibration profile later! --
2682+
histos.fill(HIST("MC/CentFT0CVsTrueMult"), recoCentrality, trueMultEta05);
26682683
}
2684+
26692685
for (const auto& mcParticle : particlesInCollision) {
2670-
if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) {
2671-
continue;
2672-
}
2686+
if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) continue;
26732687
static_for<0, 17>([&](auto i) {
2674-
fillParticleHistogramsMCRecoEvs<i>(mcParticle, collision);
2688+
// Pass the actual Reco Centrality to the helper function
2689+
fillParticleHistogramsMCRecoEvs<i>(mcParticle, collision, recoCentrality);
26752690
});
26762691
}
26772692
}
26782693

2679-
// Loop on generated collisions
2694+
// ==============================================================================
2695+
// 3. GENERATED COLLISIONS LOOP (Builds the Signal Loss Denominator)
2696+
// ==============================================================================
26802697
for (const auto& mcCollision : mcCollisions) {
26812698
if (std::abs(mcCollision.posZ()) > evselOptions.cfgCutVertex) {
26822699
continue;
26832700
}
2684-
histos.fill(HIST("MC/Multiplicity"), getMultiplicityMC(mcCollision));
2701+
26852702
const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache);
2686-
bool hasParticleInFT0C = false;
2687-
bool hasParticleInFT0A = false;
2688-
if (evselOptions.cfgINELCut.value == EvSelInelGt0Cut) {
2689-
if (!o2::pwglf::isINELgt0mc(particlesInCollision, pdgDB)) {
2690-
continue;
2703+
2704+
// -- Calculate True Multiplicity in |eta| < 0.5 --
2705+
int trueMultEta05 = 0;
2706+
for (const auto& mcParticle : particlesInCollision) {
2707+
if (!mcParticle.isPhysicalPrimary()) continue;
2708+
auto* pdgParticle = pdgDB->GetParticle(mcParticle.pdgCode());
2709+
if (pdgParticle != nullptr && std::abs(pdgParticle->Charge()) > 0.01) {
2710+
if (std::abs(mcParticle.eta()) < 0.5f) {
2711+
trueMultEta05++;
2712+
}
26912713
}
26922714
}
2693-
histos.fill(HIST("MC/MultiplicityMCINELgt0"), getMultiplicityMC(mcCollision));
2715+
2716+
histos.fill(HIST("MC/Multiplicity"), trueMultEta05);
2717+
2718+
if (evselOptions.cfgINELCut.value == EvSelInelGt0Cut) {
2719+
if (!o2::pwglf::isINELgt0mc(particlesInCollision, pdgDB)) continue;
2720+
}
2721+
histos.fill(HIST("MC/MultiplicityMCINELgt0"), trueMultEta05);
2722+
26942723
if (evselOptions.cfgINELCut.value == EvSelInelGt1Cut) {
2695-
if (!o2::pwglf::isINELgt1mc(particlesInCollision, pdgDB)) {
2696-
continue;
2697-
}
2724+
if (!o2::pwglf::isINELgt1mc(particlesInCollision, pdgDB)) continue;
26982725
}
2699-
histos.fill(HIST("MC/MultiplicityMCINELgt1"), getMultiplicityMC(mcCollision));
2726+
histos.fill(HIST("MC/MultiplicityMCINELgt1"), trueMultEta05);
2727+
27002728
for (const auto& mcParticle : particlesInCollision) {
2701-
if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) {
2702-
continue;
2703-
}
2729+
if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) continue;
27042730
static_for<0, 17>([&](auto i) {
2705-
fillParticleHistogramsMCGenEvs<i>(mcParticle, mcCollision);
2731+
// Pass True Multiplicity directly so it plots generated particles against it
2732+
fillParticleHistogramsMCGenEvs<i>(mcParticle, trueMultEta05);
27062733
});
27072734
}
2735+
27082736
if (mcCollision.isInelGt0()) {
27092737
histos.fill(HIST("MC/GenRecoCollisions"), 3.f);
27102738
}
27112739
if (mcCollision.isInelGt1()) {
27122740
histos.fill(HIST("MC/GenRecoCollisions"), 4.f);
27132741
}
2714-
if (hasParticleInFT0C && hasParticleInFT0A) {
2715-
histos.fill(HIST("MC/GenRecoCollisions"), 5.f);
2716-
}
27172742
}
27182743
}
27192744
PROCESS_SWITCH(tofSpectra, processMC, "Process MC", false);

0 commit comments

Comments
 (0)