Skip to content

Commit 9fde296

Browse files
kcui1kcui1
authored
[PWGLF] Add INELgtONE selection in hStrangeCorrelation.cxx (#17658)
1 parent 0247fa2 commit 9fde296

1 file changed

Lines changed: 135 additions & 16 deletions

File tree

‎PWGLF/Tasks/Strangeness/hStrangeCorrelation.cxx‎

Lines changed: 135 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -126,6 +126,7 @@ struct HStrangeCorrelation {
126126
Configurable<bool> doCorrelationPion{"doCorrelationPion", false, "do Pion correlation"};
127127
Configurable<bool> doGenEventSelection{"doGenEventSelection", true, "use event selections when performing closure test for the gen events"};
128128
Configurable<bool> selectINELgtZERO{"selectINELgtZERO", true, "select INEL>0 events"};
129+
Configurable<bool> selectINELgtONE{"selectINELgtONE", false, "select INEL>1 events (at least 2 charged particles in |eta| < 1)"};
129130
Configurable<float> zVertexCut{"zVertexCut", 10, "Cut on PV position"};
130131
Configurable<bool> requireAllGoodITSLayers{"requireAllGoodITSLayers", false, " require that in the event all ITS are good"};
131132
Configurable<bool> requireGoodTriggerTVX{"requireGoodTriggerTVX", false, " require acceptable FT0C-FT0A time difference"};
@@ -612,6 +613,18 @@ struct HStrangeCorrelation {
612613
bool physicalPrimary = false;
613614
};
614615

616+
// One object of a GenStudy h-K0 pair: generated kinematics plus whether it has a
617+
// reconstructed counterpart, in exactly the sense the GenStudy single-particle
618+
// folders use.
619+
struct GenStudyPairObject {
620+
float pt = 0.0f;
621+
float eta = 0.0f;
622+
float phi = 0.0f;
623+
int64_t globalIndex = -1;
624+
int64_t motherIndex = -1;
625+
bool reconstructed = false;
626+
};
627+
615628
struct PairLossTruthK0Info {
616629
int64_t globalIndex = -1;
617630
float pt = 0.0f;
@@ -723,20 +736,37 @@ struct HStrangeCorrelation {
723736
return localDensity;
724737
}
725738

726-
/// Generated-level counterpart: counts the primary charged particles in the associated pt range inside the same cone, skipping the reference particle
727-
template <typename TMcParticles>
728-
int computeLocalDensityGen(float etaRef, float phiRef, TMcParticles const& mcParticles, int64_t skipId)
739+
/// Collects the MC index of a particle and of its decay products, so that a particle never contributes to its own local density
740+
template <typename TMcParticle>
741+
void collectDescendantIds(TMcParticle const& mcParticle, std::vector<int64_t>& ids, int depth = 0)
729742
{
743+
ids.push_back(mcParticle.globalIndex());
744+
if (depth >= 3 || !mcParticle.has_daughters()) {
745+
return;
746+
}
747+
for (auto const& daughter : mcParticle.template daughters_as<aod::McParticles>()) {
748+
collectDescendantIds(daughter, ids, depth + 1);
749+
}
750+
}
751+
752+
/// Generated-level counterpart: the density the reconstruction would have measured around the generated
753+
/// direction, i.e. the same associated-quality tracks of the same collision, so that the axis means the
754+
/// same thing here as in the reconstructed histograms and can be used to correct data binned in it
755+
template <typename TMcParticle, typename TTracks>
756+
int computeLocalDensityGen(TMcParticle const& mcParticle, TTracks const& tracks)
757+
{
758+
std::vector<int64_t> skipIds;
759+
collectDescendantIds(mcParticle, skipIds);
730760
int localDensity = 0;
731-
for (auto const& mcParticle : mcParticles) {
732-
if (mcParticle.globalIndex() == skipId || !mcParticle.isPhysicalPrimary() || !isPairLossTriggerPdg(mcParticle.pdgCode())) {
761+
for (auto const& track : tracks) {
762+
if (!isValidAssocHadron(track)) {
733763
continue;
734764
}
735-
if (mcParticle.pt() < axisRanges[2][0] || mcParticle.pt() > axisRanges[2][1]) {
765+
if (track.has_mcParticle() && std::find(skipIds.begin(), skipIds.end(), track.mcParticleId()) != skipIds.end()) {
736766
continue;
737767
}
738-
double deltaEta = mcParticle.eta() - etaRef;
739-
double deltaPhi = RecoDecay::constrainAngle(mcParticle.phi() - phiRef, -PI);
768+
double deltaEta = track.eta() - mcParticle.eta();
769+
double deltaPhi = RecoDecay::constrainAngle(track.phi() - mcParticle.phi(), -PI);
740770
if (std::hypot(deltaEta, deltaPhi) < masterConfigurations.localDensityConeRadius) {
741771
localDensity++;
742772
}
@@ -2795,8 +2825,18 @@ struct HStrangeCorrelation {
27952825
// bin by bin and NotReconstructed/Gen reads directly as the loss.
27962826
histos.add("PairLossK0/GenStudy/Gen/hTrigger", "generated triggers;#it{p}_{T}^{gen} (GeV/#it{c});#eta^{gen};#varphi^{gen};#it{N}_{ch}^{gen}", kTHnF, {axesConfigurations.axisPtQA, axesConfigurations.axisEta, axesConfigurations.axisPhi, axisGenStudyNch});
27972827
histos.add("PairLossK0/GenStudy/Gen/hK0Short", "generated K0s;#it{p}_{T}^{gen} (GeV/#it{c});#eta^{gen};#varphi^{gen};#it{N}_{ch}^{gen};findable", kTHnF, {axesConfigurations.axisPtQA, axesConfigurations.axisEta, axesConfigurations.axisPhi, axisGenStudyNch, axisGenStudyFindable});
2828+
// h-K0 correlations of the very same objects. Gen/ is every generated pair and the
2829+
// four exclusive classes below split it by which of the two objects was
2830+
// reconstructed, so Reconstructed + OnlyTriggerReconstructed + OnlyK0Reconstructed
2831+
// + NotReconstructed equals Gen bin by bin.
2832+
histos.add("PairLossK0/GenStudy/Gen/hCorrelation", "generated h-K0s pairs;#Delta#eta;#Delta#varphi;#it{p}_{T}^{trigger} (GeV/#it{c});#it{p}_{T}^{K^{0}_{S}} (GeV/#it{c});#it{N}_{ch}^{gen}", kTHnF, {axisDeltaEtaNDim, axisDeltaPhiNDim, axisPtTriggerNDim, axisPtAssocNDim, axisGenStudyNch});
27982833
histos.addClone("PairLossK0/GenStudy/Gen/", "PairLossK0/GenStudy/Reconstructed/");
27992834
histos.addClone("PairLossK0/GenStudy/Gen/", "PairLossK0/GenStudy/NotReconstructed/");
2835+
// Only the correlation exists for the two mixed classes -- a single particle is
2836+
// either reconstructed or not, so cloning the single-particle folders here would
2837+
// only produce histograms with no meaning.
2838+
histos.add("PairLossK0/GenStudy/OnlyTriggerReconstructed/hCorrelation", "h-K0s pairs with only the trigger reconstructed;#Delta#eta;#Delta#varphi;#it{p}_{T}^{trigger} (GeV/#it{c});#it{p}_{T}^{K^{0}_{S}} (GeV/#it{c});#it{N}_{ch}^{gen}", kTHnF, {axisDeltaEtaNDim, axisDeltaPhiNDim, axisPtTriggerNDim, axisPtAssocNDim, axisGenStudyNch});
2839+
histos.add("PairLossK0/GenStudy/OnlyK0Reconstructed/hCorrelation", "h-K0s pairs with only the K0s reconstructed;#Delta#eta;#Delta#varphi;#it{p}_{T}^{trigger} (GeV/#it{c});#it{p}_{T}^{K^{0}_{S}} (GeV/#it{c});#it{N}_{ch}^{gen}", kTHnF, {axisDeltaEtaNDim, axisDeltaPhiNDim, axisPtTriggerNDim, axisPtAssocNDim, axisGenStudyNch});
28002840

28012841
for (auto const& histogram : {histos.get<THn>(HIST("PairLossK0/GenStudy/Gen/hK0Short")),
28022842
histos.get<THn>(HIST("PairLossK0/GenStudy/Reconstructed/hK0Short")),
@@ -3018,6 +3058,8 @@ struct HStrangeCorrelation {
30183058
}
30193059
histos.addClone("Generated/", "GeneratedWithPV/");
30203060

3061+
// The density axis is the reconstructed one on both sides: these generated histograms are counted
3062+
// from the tracks of the best collision, exactly like their reconstructed counterparts.
30213063
if (masterConfigurations.doLocalDensityStudy && masterConfigurations.doPPAnalysis) {
30223064
histos.add("GeneratedWithPV/hTriggerLocalDensity", "", kTH3F, {axesConfigurations.axisPtQA, axesConfigurations.axisEta, axesConfigurations.axisLocalDensity});
30233065
for (int i = 0; i < AssocParticleTypesNoHadron; i++) {
@@ -3194,6 +3236,9 @@ struct HStrangeCorrelation {
31943236
if (!collision.isInelGt0() && masterConfigurations.selectINELgtZERO) {
31953237
return false;
31963238
}
3239+
if (!collision.isInelGt1() && masterConfigurations.selectINELgtONE) {
3240+
return false;
3241+
}
31973242
if (!collision.selection_bit(aod::evsel::kIsGoodITSLayersAll) && masterConfigurations.requireAllGoodITSLayers) {
31983243
return false;
31993244
}
@@ -4095,7 +4140,7 @@ struct HStrangeCorrelation {
40954140
}
40964141
}
40974142

4098-
void processMCGenerated(aod::McCollision const& /*mcCollision*/, soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSels, aod::CentFT0Ms, aod::CentFT0Cs, aod::PVMults>> const& collisions, aod::McParticles const& mcParticles)
4143+
void processMCGenerated(aod::McCollision const& /*mcCollision*/, soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSels, aod::CentFT0Ms, aod::CentFT0Cs, aod::PVMults>> const& collisions, aod::McParticles const& mcParticles, TracksCompleteMC const& tracks)
40994144
{
41004145
histos.fill(HIST("hClosureTestEventCounter"), 2.5f);
41014146

@@ -4143,17 +4188,20 @@ struct HStrangeCorrelation {
41434188

41444189
// determine best collision properties
41454190
int biggestNContribs = -1;
4191+
int64_t bestCollisionId = -1;
41464192
float bestCollisionFT0Mpercentile = -1;
41474193
float bestCollisionFT0Cpercentile = -1;
41484194
float bestCollisionVtxZ = 0.0f;
41494195
bool bestCollisionSel8 = false;
41504196
bool bestCollisionINELgtZERO = false;
4197+
bool bestCollisionINELgtONE = false;
41514198
bool isCollisionSelect = false;
41524199
uint32_t bestCollisionTriggerPresenceMap = 0;
41534200

41544201
for (auto const& collision : collisions) {
41554202
if (biggestNContribs < collision.numContrib()) {
41564203
biggestNContribs = collision.numContrib();
4204+
bestCollisionId = collision.globalIndex();
41574205
bestCollisionFT0Mpercentile = collision.centFT0M();
41584206
bestCollisionFT0Cpercentile = collision.centFT0C();
41594207
if (masterConfigurations.applyNewMCSelection) {
@@ -4162,6 +4210,7 @@ struct HStrangeCorrelation {
41624210
bestCollisionSel8 = collision.sel8();
41634211
bestCollisionVtxZ = collision.posZ();
41644212
bestCollisionINELgtZERO = collision.isInelGt0();
4213+
bestCollisionINELgtONE = collision.isInelGt1();
41654214
}
41664215
if (triggerPresenceMap.size() > 0) {
41674216
bestCollisionTriggerPresenceMap = triggerPresenceMap[collision.globalIndex()];
@@ -4210,6 +4259,9 @@ struct HStrangeCorrelation {
42104259
if (!bestCollisionINELgtZERO) {
42114260
return;
42124261
}
4262+
if (masterConfigurations.selectINELgtONE && !bestCollisionINELgtONE) {
4263+
return;
4264+
}
42134265
}
42144266

42154267
histos.fill(HIST("hClosureTestEventCounter"), 3.5f);
@@ -4233,6 +4285,10 @@ struct HStrangeCorrelation {
42334285
}
42344286
}
42354287

4288+
// The local density of a generated particle is counted from the tracks of the best collision,
4289+
// the same collision whose centrality already labels these generated histograms.
4290+
const auto bestCollisionTracks = tracks.sliceBy(pairLossTracksPerCollision, bestCollisionId);
4291+
42364292
for (auto const& mcParticle : mcParticles) {
42374293
if (doAssocPhysicalPrimaryInGen && !mcParticle.isPhysicalPrimary()) {
42384294
continue;
@@ -4246,7 +4302,7 @@ struct HStrangeCorrelation {
42464302
histos.fill(HIST("GeneratedWithPV/hTrigger"), gpt, geta, bestCollisionFT0Mpercentile);
42474303
}
42484304
if (masterConfigurations.doLocalDensityStudy && masterConfigurations.doPPAnalysis) {
4249-
histos.fill(HIST("GeneratedWithPV/hTriggerLocalDensity"), gpt, geta, computeLocalDensityGen(geta, mcParticle.phi(), mcParticles, mcParticle.globalIndex()));
4305+
histos.fill(HIST("GeneratedWithPV/hTriggerLocalDensity"), gpt, geta, computeLocalDensityGen(mcParticle, bestCollisionTracks));
42504306
}
42514307
if (mcParticle.pdgCode() > 0) {
42524308
histos.fill(HIST("GeneratedWithPV/hPositiveTrigger"), gpt, geta, bestCollisionFT0Mpercentile);
@@ -4314,7 +4370,7 @@ struct HStrangeCorrelation {
43144370
histos.fill(HIST("GeneratedWithPV/h") + HIST(Particlenames[Index]) + HIST("_MidYVsMult"), gpt, bestCollisionFT0Mpercentile);
43154371
}
43164372
if (masterConfigurations.doLocalDensityStudy && masterConfigurations.doPPAnalysis) {
4317-
histos.fill(HIST("GeneratedWithPV/h") + HIST(Particlenames[Index]) + HIST("LocalDensity"), gpt, geta, computeLocalDensityGen(geta, mcParticle.phi(), mcParticles, mcParticle.globalIndex()));
4373+
histos.fill(HIST("GeneratedWithPV/h") + HIST(Particlenames[Index]) + HIST("LocalDensity"), gpt, geta, computeLocalDensityGen(mcParticle, bestCollisionTracks));
43184374
}
43194375
}
43204376
});
@@ -4356,8 +4412,15 @@ struct HStrangeCorrelation {
43564412
histos.fill(HIST("PairLossK0/GenStudy/hEventCounter"), 0.0f);
43574413

43584414
// Generated-level event selection. No reconstructed variable is used.
4359-
if (masterConfigurations.selectINELgtZERO && !o2::pwglf::isINELgt0mc(mcParticles, pdgDB)) {
4360-
return;
4415+
// INEL>1 implies INEL>0, so only the tighter enabled selection has to be evaluated
4416+
if (masterConfigurations.selectINELgtONE) {
4417+
if (!o2::pwglf::isINELgt1mc(mcParticles, pdgDB)) {
4418+
return;
4419+
}
4420+
} else if (masterConfigurations.selectINELgtZERO) {
4421+
if (!o2::pwglf::isINELgt0mc(mcParticles, pdgDB)) {
4422+
return;
4423+
}
43614424
}
43624425
histos.fill(HIST("PairLossK0/GenStudy/hEventCounter"), 1.0f);
43634426
if (std::abs(mcCollision.posZ()) > masterConfigurations.zVertexCut) {
@@ -4392,6 +4455,8 @@ struct HStrangeCorrelation {
43924455
// them.
43934456
std::unordered_set<int64_t> reconstructedTrackMcIds;
43944457
std::unordered_set<int64_t> reconstructedV0McIds;
4458+
std::vector<GenStudyPairObject> genStudyTriggers;
4459+
std::vector<GenStudyPairObject> genStudyK0s;
43954460
for (auto const& collision : recCollisions) {
43964461
const auto trackSlice = tracks.sliceBy(pairLossTracksPerCollision, collision.globalIndex());
43974462
for (auto const& track : trackSlice) {
@@ -4433,6 +4498,13 @@ struct HStrangeCorrelation {
44334498
} else {
44344499
histos.fill(HIST("PairLossK0/GenStudy/NotReconstructed/hTrigger"), genPt, genEta, genPhi, generatedNch);
44354500
}
4501+
genStudyTriggers.push_back(GenStudyPairObject{
4502+
.pt = genPt,
4503+
.eta = genEta,
4504+
.phi = genPhi,
4505+
.globalIndex = static_cast<int64_t>(mcParticle.globalIndex()),
4506+
.motherIndex = mcParticle.has_mothers() ? static_cast<int64_t>(mcParticle.mothers_first_as<aod::McParticles>().globalIndex()) : -1,
4507+
.reconstructed = reconstructedTrackMcIds.count(mcParticle.globalIndex()) > 0});
44364508
}
44374509
}
44384510

@@ -4470,6 +4542,38 @@ struct HStrangeCorrelation {
44704542
} else {
44714543
histos.fill(HIST("PairLossK0/GenStudy/NotReconstructed/hK0Short"), genPt, genEta, genPhi, generatedNch, k0Findable);
44724544
}
4545+
genStudyK0s.push_back(GenStudyPairObject{
4546+
.pt = genPt,
4547+
.eta = genEta,
4548+
.phi = genPhi,
4549+
.globalIndex = static_cast<int64_t>(mcParticle.globalIndex()),
4550+
.motherIndex = -1,
4551+
.reconstructed = reconstructedV0McIds.count(mcParticle.globalIndex()) > 0});
4552+
}
4553+
}
4554+
4555+
// h-K0 correlations of the objects collected above, in generated coordinates.
4556+
// Same delta-phi / delta-eta convention as every other correlation in this task
4557+
// (trigger minus associated), and the same autocorrelation rejection: a trigger
4558+
// that is a decay product of the K0 it would be paired with is skipped.
4559+
// Every pair goes into Gen/ and into exactly one of the four exclusive classes.
4560+
for (auto const& trigger : genStudyTriggers) {
4561+
for (auto const& k0 : genStudyK0s) {
4562+
if (trigger.globalIndex == k0.globalIndex || trigger.motherIndex == k0.globalIndex) {
4563+
continue;
4564+
}
4565+
const float deltaPhi = computeDeltaPhi(trigger.phi, k0.phi);
4566+
const float deltaEta = trigger.eta - k0.eta;
4567+
histos.fill(HIST("PairLossK0/GenStudy/Gen/hCorrelation"), deltaEta, deltaPhi, trigger.pt, k0.pt, generatedNch);
4568+
if (trigger.reconstructed && k0.reconstructed) {
4569+
histos.fill(HIST("PairLossK0/GenStudy/Reconstructed/hCorrelation"), deltaEta, deltaPhi, trigger.pt, k0.pt, generatedNch);
4570+
} else if (trigger.reconstructed) {
4571+
histos.fill(HIST("PairLossK0/GenStudy/OnlyTriggerReconstructed/hCorrelation"), deltaEta, deltaPhi, trigger.pt, k0.pt, generatedNch);
4572+
} else if (k0.reconstructed) {
4573+
histos.fill(HIST("PairLossK0/GenStudy/OnlyK0Reconstructed/hCorrelation"), deltaEta, deltaPhi, trigger.pt, k0.pt, generatedNch);
4574+
} else {
4575+
histos.fill(HIST("PairLossK0/GenStudy/NotReconstructed/hCorrelation"), deltaEta, deltaPhi, trigger.pt, k0.pt, generatedNch);
4576+
}
44734577
}
44744578
}
44754579
};
@@ -4601,6 +4705,7 @@ struct HStrangeCorrelation {
46014705
float genBestCollisionVtxZ = 0.0f;
46024706
bool genBestCollisionSel8 = false;
46034707
bool genBestCollisionINELgtZERO = false;
4708+
bool genBestCollisionINELgtONE = false;
46044709
bool genCollisionSelected = false;
46054710
int genLargestNContributors = -1;
46064711
uint32_t genBestCollisionTriggerPresenceMap = 0;
@@ -4617,6 +4722,7 @@ struct HStrangeCorrelation {
46174722
genBestCollisionSel8 = recCollision.sel8();
46184723
genBestCollisionVtxZ = recCollision.posZ();
46194724
genBestCollisionINELgtZERO = recCollision.isInelGt0();
4725+
genBestCollisionINELgtONE = recCollision.isInelGt1();
46204726
}
46214727
if (triggerPresenceMap.size() > 0) {
46224728
genBestCollisionTriggerPresenceMap = triggerPresenceMap[recCollision.globalIndex()];
@@ -4627,7 +4733,8 @@ struct HStrangeCorrelation {
46274733
genEventSelected = genEventSelected && genCollisionSelected;
46284734
} else if (masterConfigurations.doGenEventSelection) {
46294735
genEventSelected = genEventSelected && genBestCollisionSel8 && std::abs(genBestCollisionVtxZ) <= masterConfigurations.zVertexCut &&
4630-
genBestCollisionINELgtZERO && genBestCollisionMultiplicity >= axisRanges[5][0] && genBestCollisionMultiplicity <= axisRanges[5][1];
4736+
genBestCollisionINELgtZERO && (!masterConfigurations.selectINELgtONE || genBestCollisionINELgtONE) &&
4737+
genBestCollisionMultiplicity >= axisRanges[5][0] && genBestCollisionMultiplicity <= axisRanges[5][1];
46314738
}
46324739

46334740
if (genEventSelected && masterConfigurations.doCorrelationK0Short && TESTBIT(doCorrelation, IndexK0)) {
@@ -5852,6 +5959,7 @@ struct HStrangeCorrelation {
58525959
float bestCollisionVtxZ = 0.0f;
58535960
bool bestCollisionSel8 = false;
58545961
bool bestCollisionINELgtZERO = false;
5962+
bool bestCollisionINELgtONE = false;
58555963
bool isCollisionSelect = false;
58565964
int biggestNContribs = -1;
58575965
uint32_t bestCollisionTriggerPresenceMap = 0;
@@ -5866,6 +5974,7 @@ struct HStrangeCorrelation {
58665974
bestCollisionSel8 = recCollision.sel8();
58675975
bestCollisionVtxZ = recCollision.posZ();
58685976
bestCollisionINELgtZERO = recCollision.isInelGt0();
5977+
bestCollisionINELgtONE = recCollision.isInelGt1();
58695978
}
58705979
if (triggerPresenceMap.size() > 0) {
58715980
bestCollisionTriggerPresenceMap = triggerPresenceMap[recCollision.globalIndex()];
@@ -5893,6 +6002,9 @@ struct HStrangeCorrelation {
58936002
if (!bestCollisionINELgtZERO) {
58946003
return;
58956004
}
6005+
if (masterConfigurations.selectINELgtONE && !bestCollisionINELgtONE) {
6006+
return;
6007+
}
58966008
if (bestCollisionCentpercentile > axisRanges[5][1] || bestCollisionCentpercentile < axisRanges[5][0]) {
58976009
return;
58986010
}
@@ -6222,8 +6334,15 @@ struct HStrangeCorrelation {
62226334
float multEta08 = -1;
62236335
float multEta05 = -1;
62246336
histos.fill(HIST("Prediction/hEventSelection"), 0.5);
6225-
if (masterConfigurations.selectINELgtZERO && !o2::pwglf::isINELgt0mc(mcParticles, pdgDB)) {
6226-
return;
6337+
// INEL>1 implies INEL>0, so only the tighter enabled selection has to be evaluated
6338+
if (masterConfigurations.selectINELgtONE) {
6339+
if (!o2::pwglf::isINELgt1mc(mcParticles, pdgDB)) {
6340+
return;
6341+
}
6342+
} else if (masterConfigurations.selectINELgtZERO) {
6343+
if (!o2::pwglf::isINELgt0mc(mcParticles, pdgDB)) {
6344+
return;
6345+
}
62276346
}
62286347
histos.fill(HIST("Prediction/hEventSelection"), 1.5);
62296348
if (std::abs(mcCollision.posZ()) > masterConfigurations.zVertexCut) {

0 commit comments

Comments
 (0)