Skip to content

Commit 9d7a963

Browse files
alcalivaalcaliva
committed
Fix proton double counting in coalescence and add jet/UE spectrum reweighting
1 parent 1a48c05 commit 9d7a963

1 file changed

Lines changed: 105 additions & 15 deletions

File tree

‎PWGLF/Tasks/Nuspex/antinucleiInJets.cxx‎

Lines changed: 105 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -661,6 +661,31 @@ struct AntinucleiInJets {
661661
return (particle.isPhysicalPrimary() || fromHF);
662662
}
663663

664+
// Check if particle is from heavy-flavor hadron decay
665+
bool isFromHF(aod::McParticle const& particle, aod::McParticles const& mcParticles)
666+
{
667+
// Keep only pi, K, p, d, e, mu
668+
int pdg = std::abs(particle.pdgCode());
669+
if (!(pdg == PDG_t::kPiPlus || pdg == PDG_t::kKPlus || pdg == PDG_t::kProton || pdg == o2::constants::physics::Pdg::kDeuteron || pdg == PDG_t::kElectron || pdg == PDG_t::kMuonMinus))
670+
return false;
671+
672+
// Constants for identifying heavy-flavor (charm and bottom) content from PDG codes
673+
static constexpr int CharmQuark = 4;
674+
static constexpr int BottomQuark = 5;
675+
static constexpr int Hundreds = 100;
676+
static constexpr int Thousands = 1000;
677+
678+
// Check if particle is from heavy-flavor decay
679+
bool fromHF = false;
680+
if (particle.has_mothers()) {
681+
auto mother = mcParticles.iteratorAt(particle.mothersIds()[0]);
682+
int motherPdg = std::abs(mother.pdgCode());
683+
fromHF = (motherPdg / Hundreds == CharmQuark || motherPdg / Hundreds == BottomQuark || motherPdg / Thousands == CharmQuark || motherPdg / Thousands == BottomQuark);
684+
}
685+
686+
return fromHF;
687+
}
688+
664689
// Evaluate proton–neutron coalescence for deuteron formation
665690
template <typename ReducedPart>
666691
bool passDeuteronCoalescence(const ReducedPart& p, const ReducedPart& n, double p0, TRandom3& random)
@@ -2932,6 +2957,7 @@ struct AntinucleiInJets {
29322957
// Define per-event particle containers
29332958
std::vector<ReducedParticle> chargedParticles;
29342959
std::vector<ReducedParticle> protons;
2960+
std::vector<ReducedParticle> protonsFromHF;
29352961
std::vector<ReducedParticle> neutrons;
29362962
std::vector<fastjet::PseudoJet> fjParticles;
29372963

@@ -2945,6 +2971,7 @@ struct AntinucleiInJets {
29452971
// Clear containers at the start of the event loop
29462972
chargedParticles.clear();
29472973
protons.clear();
2974+
protonsFromHF.clear();
29482975
neutrons.clear();
29492976
fjParticles.clear();
29502977

@@ -2972,6 +2999,11 @@ struct AntinucleiInJets {
29722999
protons.push_back({particle.px(), particle.py(), particle.pz(), particle.pdgCode(), mcId, false});
29733000
}
29743001

3002+
// Store Protons from HF decays
3003+
if (isFromHF(particle, mcParticles) && std::abs(particle.pdgCode()) == PDG_t::kProton) {
3004+
protonsFromHF.push_back({particle.px(), particle.py(), particle.pz(), particle.pdgCode(), mcId, false});
3005+
}
3006+
29753007
// Store Neutrons
29763008
if (particle.isPhysicalPrimary() && std::abs(particle.pdgCode()) == PDG_t::kNeutron) {
29773009
neutrons.push_back({particle.px(), particle.py(), particle.pz(), particle.pdgCode(), mcId, false});
@@ -2981,17 +3013,12 @@ struct AntinucleiInJets {
29813013
if (!isPhysicalPrimaryOrFromHF(particle, mcParticles))
29823014
continue;
29833015

2984-
// Select particles within acceptance
2985-
if (particle.eta() < minEta || particle.eta() > maxEta || particle.pt() < MinPtParticle)
3016+
// Select particles within acceptance and skip protons (will be added after coalescence)
3017+
if (particle.eta() < minEta || particle.eta() > maxEta || particle.pt() < MinPtParticle || std::abs(particle.pdgCode()) == PDG_t::kProton)
29863018
continue;
29873019
chargedParticles.push_back({particle.px(), particle.py(), particle.pz(), particle.pdgCode(), mcId, false});
29883020
}
29893021

2990-
// Reject empty events
2991-
if (chargedParticles.empty())
2992-
continue;
2993-
registryMC.fill(HIST("genEventsCoalescence"), 2.5);
2994-
29953022
// Build deuterons
29963023
for (int ip = 0; ip < static_cast<int>(protons.size()); ip++) {
29973024
auto& proton = protons[ip];
@@ -3005,6 +3032,10 @@ struct AntinucleiInJets {
30053032

30063033
if (passDeuteronCoalescence(proton, neutron, coalescenceMomentum, mRand)) {
30073034

3035+
// If proton and neutron are used by coalescence they are removed
3036+
neutron.used = true;
3037+
proton.used = true;
3038+
30083039
int sign = (proton.pdgCode > 0) ? +1 : -1;
30093040
int deuteronPdg = sign * o2::constants::physics::Pdg::kDeuteron;
30103041

@@ -3014,18 +3045,40 @@ struct AntinucleiInJets {
30143045
double energyDeut = std::sqrt(pxDeut * pxDeut + pyDeut * pyDeut + pzDeut * pzDeut + massDeut * massDeut);
30153046
LorentzVector pd(pxDeut, pyDeut, pzDeut, energyDeut);
30163047
if (pd.Eta() < minEta || pd.Eta() > maxEta || pd.Pt() < MinPtParticle)
3017-
continue;
3048+
break;
30183049

30193050
// Store Deuteron
30203051
chargedParticles.push_back({pxDeut, pyDeut, pzDeut, deuteronPdg, proton.mcIndex, false});
3021-
3022-
neutron.used = true;
3023-
proton.used = true;
30243052
break;
30253053
}
30263054
}
30273055
}
30283056

3057+
// Add unused protons
3058+
for (int ip = 0; ip < static_cast<int>(protons.size()); ip++) {
3059+
auto& proton = protons[ip];
3060+
if (proton.used)
3061+
continue;
3062+
if (proton.eta() < minEta || proton.eta() > maxEta || proton.pt() < MinPtParticle)
3063+
continue;
3064+
chargedParticles.push_back({proton.px, proton.py, proton.pz, proton.pdgCode, proton.mcIndex, proton.used});
3065+
}
3066+
3067+
// Add protons from HF decays
3068+
for (int ip = 0; ip < static_cast<int>(protonsFromHF.size()); ip++) {
3069+
auto& proton = protonsFromHF[ip];
3070+
if (proton.used)
3071+
continue;
3072+
if (proton.eta() < minEta || proton.eta() > maxEta || proton.pt() < MinPtParticle)
3073+
continue;
3074+
chargedParticles.push_back({proton.px, proton.py, proton.pz, proton.pdgCode, proton.mcIndex, proton.used});
3075+
}
3076+
3077+
// Reject empty events
3078+
if (chargedParticles.empty())
3079+
continue;
3080+
registryMC.fill(HIST("genEventsCoalescence"), 2.5);
3081+
30293082
for (const auto& part : chargedParticles) {
30303083
if (part.used)
30313084
continue;
@@ -3086,10 +3139,29 @@ struct AntinucleiInJets {
30863139
// Analyze jet constituents
30873140
std::vector<fastjet::PseudoJet> jetConstituents = jet.constituents();
30883141
for (const auto& particle : jetConstituents) {
3142+
3143+
// Calculate weight for antiprotons
3144+
double wAntipJet = 1.0;
3145+
if (applyReweighting) {
3146+
if (particle.pt() >= antiprotonsInsideJets->GetXaxis()->GetXmin() && particle.pt() < antiprotonsInsideJets->GetXaxis()->GetXmax()) {
3147+
int ipt = antiprotonsInsideJets->FindBin(particle.pt());
3148+
wAntipJet = antiprotonsInsideJets->GetBinContent(ipt);
3149+
}
3150+
}
3151+
3152+
// Calculate weight for antideuterons
3153+
double wAntidJet = 1.0;
3154+
if (applyReweighting) {
3155+
if (0.5 * particle.pt() >= antiprotonsInsideJets->GetXaxis()->GetXmin() && 0.5 * particle.pt() < antiprotonsInsideJets->GetXaxis()->GetXmax()) {
3156+
int ipt = antiprotonsInsideJets->FindBin(0.5 * particle.pt());
3157+
wAntidJet = antiprotonsInsideJets->GetBinContent(ipt) * antiprotonsInsideJets->GetBinContent(ipt);
3158+
}
3159+
}
3160+
30893161
if (particle.user_index() == PDG_t::kProtonBar)
3090-
registryMC.fill(HIST("antiproton_coal_jet"), particle.pt());
3162+
registryMC.fill(HIST("antiproton_coal_jet"), particle.pt(), wAntipJet);
30913163
if (particle.user_index() == -o2::constants::physics::Pdg::kDeuteron)
3092-
registryMC.fill(HIST("antideuteron_coal_jet"), particle.pt());
3164+
registryMC.fill(HIST("antideuteron_coal_jet"), particle.pt(), wAntidJet);
30933165
}
30943166

30953167
// Set up two perpendicular cone axes for underlying event estimation
@@ -3119,11 +3191,29 @@ struct AntinucleiInJets {
31193191
if (deltaRUe1 > rJet && deltaRUe2 > rJet)
31203192
continue;
31213193

3194+
// Calculate weight for antiprotons
3195+
double wAntipUE = 1.0;
3196+
if (applyReweighting) {
3197+
if (chParticle.pt() >= antiprotonsPerpCone->GetXaxis()->GetXmin() && chParticle.pt() < antiprotonsPerpCone->GetXaxis()->GetXmax()) {
3198+
int ipt = antiprotonsPerpCone->FindBin(chParticle.pt());
3199+
wAntipUE = antiprotonsPerpCone->GetBinContent(ipt);
3200+
}
3201+
}
3202+
3203+
// Calculate weight for antideuterons
3204+
double wAntidUE = 1.0;
3205+
if (applyReweighting) {
3206+
if (0.5 * chParticle.pt() >= antiprotonsPerpCone->GetXaxis()->GetXmin() && 0.5 * chParticle.pt() < antiprotonsPerpCone->GetXaxis()->GetXmax()) {
3207+
int ipt = antiprotonsPerpCone->FindBin(0.5 * chParticle.pt());
3208+
wAntidUE = antiprotonsPerpCone->GetBinContent(ipt) * antiprotonsPerpCone->GetBinContent(ipt);
3209+
}
3210+
}
3211+
31223212
// Fill histograms for UE
31233213
if (chParticle.pdgCode == PDG_t::kProtonBar)
3124-
registryMC.fill(HIST("antiproton_coal_ue"), chParticle.pt());
3214+
registryMC.fill(HIST("antiproton_coal_ue"), chParticle.pt(), wAntipUE);
31253215
if (chParticle.pdgCode == -o2::constants::physics::Pdg::kDeuteron)
3126-
registryMC.fill(HIST("antideuteron_coal_ue"), chParticle.pt());
3216+
registryMC.fill(HIST("antideuteron_coal_ue"), chParticle.pt(), wAntidUE);
31273217
}
31283218
}
31293219
if (isAtLeastOneJetSelected) {

0 commit comments

Comments
 (0)