Skip to content

Commit 5756773

Browse files
scannitoscannito
authored
[PWGLF] More protections for getPhiPurity function calling (#10183)
1 parent a5f7ca2 commit 5756773

File tree

1 file changed

+56
-22
lines changed

1 file changed

+56
-22
lines changed

‎PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx‎

Lines changed: 56 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -111,6 +111,7 @@ struct Phik0shortanalysis {
111111
// Configurables on phi pT bins
112112
Configurable<std::vector<double>> binspTPhi{"binspTPhi", {0.4, 0.8, 1.4, 2.0, 2.8, 4.0, 6.0, 10.0}, "pT bin limits for Phi"};
113113
Configurable<float> minPhiPt{"minPhiPt", 0.4f, "Minimum pT for Phi"};
114+
Configurable<float> maxPhiPt{"maxPhiPt", 10.0f, "Maximum pT for Phi"};
114115

115116
// Configurables on phi mass
116117
Configurable<int> nBinsMPhi{"nBinsMPhi", 13, "N bins in cfgmassPhiaxis"};
@@ -139,7 +140,7 @@ struct Phik0shortanalysis {
139140
Configurable<float> upMK0S{"upMK0S", 0.52f, "Upper limit on K0Short mass"};
140141

141142
// Configurable on K0S pT bins
142-
Configurable<std::vector<double>> binspTK0S{"binspTK0S", {0.1, 0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0, 6.0}, "pT bin limits for K0S"};
143+
Configurable<std::vector<double>> binspTK0S{"binspTK0S", {0.1, 0.5, 0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0, 6.0}, "pT bin limits for K0S"};
143144

144145
// Configurable on pion pT bins
145146
Configurable<std::vector<double>> binspTPi{"binspTPi", {0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.5, 2.0, 3.0}, "pT bin limits for pions"};
@@ -740,16 +741,27 @@ struct Phik0shortanalysis {
740741
// Get the phi purity choosing the correct purity function according to the multiplicity and pt of the phi
741742
double getPhiPurity(float multiplicity, const ROOT::Math::PxPyPzMVector& Phi)
742743
{
743-
// Find multiplicity bin using lower_bound
744-
auto multIt = std::lower_bound(binsMult->begin(), binsMult->end(), multiplicity);
745-
auto multIdx = multIt != binsMult->end() ? std::distance(binsMult->begin(), multIt) - 1 : -1;
744+
// Check if multiplicity is out of range
745+
if (multiplicity < binsMult->front() || multiplicity >= binsMult->back()) {
746+
LOG(info) << "Multiplicity out of range: " << multiplicity;
747+
return 0;
748+
}
749+
750+
// Find the multiplicity bin using upper_bound which finds the first element strictly greater than 'multiplicity'
751+
// Subtract 1 to get the correct bin index
752+
auto multIt = std::upper_bound(binsMult->begin(), binsMult->end(), multiplicity);
753+
int multIdx = std::distance(binsMult->begin(), multIt) - 1;
746754

747-
// Find phi-pT bin using lower_bound
748-
auto pTIt = std::lower_bound(binspTPhi->begin(), binspTPhi->end(), Phi.Pt());
749-
auto pTIdx = pTIt != binspTPhi->end() ? std::distance(binspTPhi->begin(), pTIt) - 1 : -1;
755+
// Check if pT is out of range
756+
if (Phi.Pt() < binspTPhi->front() || Phi.Pt() >= binspTPhi->back()) {
757+
LOG(info) << "pT out of range: " << Phi.Pt();
758+
return 0;
759+
}
750760

751-
if (multIdx == -1 || pTIdx == -1)
752-
LOG(fatal) << "Problem computing phi purity!";
761+
// Find the pT bin using upper_bound
762+
// The logic is the same as for multiplicity
763+
auto pTIt = std::upper_bound(binspTPhi->begin(), binspTPhi->end(), Phi.Pt());
764+
int pTIdx = std::distance(binspTPhi->begin(), pTIt) - 1;
753765

754766
return phiPurityFunctions[multIdx][pTIdx]->Eval(Phi.M());
755767
}
@@ -854,6 +866,8 @@ struct Phik0shortanalysis {
854866
bool isCountedPhi = false;
855867
bool isFilledhV0 = false;
856868

869+
double weight{1.0};
870+
857871
// Loop over all positive tracks
858872
for (const auto& track1 : posThisColl) {
859873
if (!selectionTrackResonance<false>(track1, true) || !selectionPIDKaonpTdependent(track1))
@@ -875,7 +889,7 @@ struct Phik0shortanalysis {
875889
continue; // condition to avoid double counting of pair
876890

877891
ROOT::Math::PxPyPzMVector recPhi = recMother(track1, track2, massKa, massKa);
878-
if (recPhi.Pt() < minPhiPt)
892+
if (recPhi.Pt() < minPhiPt || recPhi.Pt() > maxPhiPt)
879893
continue;
880894
if (std::abs(recPhi.Rapidity()) > cfgYAcceptance)
881895
continue;
@@ -886,6 +900,9 @@ struct Phik0shortanalysis {
886900
isCountedPhi = true;
887901
}
888902

903+
if (fillMethodSingleWeight)
904+
weight *= (1 - getPhiPurity(multiplicity, recPhi));
905+
889906
dataPhiHist.fill(HIST("h3PhipurInvMass"), multiplicity, recPhi.Pt(), recPhi.M());
890907

891908
std::array<bool, 3> isCountedK0S{false, false, false};
@@ -962,6 +979,9 @@ struct Phik0shortanalysis {
962979
}
963980
}
964981
}
982+
983+
weight = 1 - weight;
984+
dataEventHist.fill(HIST("hEventSelection"), 5, weight); // at least a Phi in the event
965985
}
966986

967987
PROCESS_SWITCH(Phik0shortanalysis, processQAPurity, "Process for QA and Phi Purities", true);
@@ -1016,7 +1036,7 @@ struct Phik0shortanalysis {
10161036
continue; // condition to avoid double counting of pair
10171037

10181038
ROOT::Math::PxPyPzMVector recPhi = recMother(track1, track2, massKa, massKa);
1019-
if (recPhi.Pt() < minPhiPt)
1039+
if (recPhi.Pt() < minPhiPt || recPhi.Pt() > maxPhiPt)
10201040
continue;
10211041
if (recPhi.M() < lowMPhi || recPhi.M() > upMPhi)
10221042
continue;
@@ -1105,7 +1125,7 @@ struct Phik0shortanalysis {
11051125
continue; // condition to avoid double counting of pair
11061126

11071127
ROOT::Math::PxPyPzMVector recPhi = recMother(track1, track2, massKa, massKa);
1108-
if (recPhi.Pt() < minPhiPt)
1128+
if (recPhi.Pt() < minPhiPt || recPhi.Pt() > maxPhiPt)
11091129
continue;
11101130
if (recPhi.M() < lowMPhi || recPhi.M() > upMPhi)
11111131
continue;
@@ -1213,7 +1233,7 @@ struct Phik0shortanalysis {
12131233
continue;
12141234

12151235
ROOT::Math::PxPyPzMVector recPhi = recMother(track1, track2, massKa, massKa);
1216-
if (recPhi.Pt() < minPhiPt)
1236+
if (recPhi.Pt() < minPhiPt || recPhi.Pt() > maxPhiPt)
12171237
continue;
12181238

12191239
mcPhiHist.fill(HIST("h3PhiRapiditySmearing"), genmultiplicity, recPhi.Rapidity(), mcMotherPhi.y());
@@ -1355,7 +1375,9 @@ struct Phik0shortanalysis {
13551375
for (const auto& mcParticle : mcParticlesThisColl) {
13561376
if (mcParticle.pdgCode() != 333)
13571377
continue;
1358-
if (std::abs(mcParticle.y()) > cfgYAcceptance || mcParticle.pt() < minPhiPt)
1378+
if (mcParticle.pt() < minPhiPt || mcParticle.pt() > maxPhiPt)
1379+
continue;
1380+
if (std::abs(mcParticle.y()) > cfgYAcceptance)
13591381
continue;
13601382

13611383
if (!isCountedMCPhi.at(0)) {
@@ -1436,7 +1458,9 @@ struct Phik0shortanalysis {
14361458
for (const auto& mcParticle : mcParticlesThisColl) {
14371459
if (mcParticle.pdgCode() != 333)
14381460
continue;
1439-
if (std::abs(mcParticle.y()) > cfgYAcceptance || mcParticle.pt() < minPhiPt)
1461+
if (mcParticle.pt() < minPhiPt || mcParticle.pt() > maxPhiPt)
1462+
continue;
1463+
if (std::abs(mcParticle.y()) > cfgYAcceptance)
14401464
continue;
14411465

14421466
if (!isCountedMCPhi.at(0)) {
@@ -1467,7 +1491,9 @@ struct Phik0shortanalysis {
14671491
for (const auto& mcParticle : mcParticlesThisColl) {
14681492
if (mcParticle.pdgCode() != 333)
14691493
continue;
1470-
if (std::abs(mcParticle.y()) > cfgYAcceptance || mcParticle.pt() < minPhiPt)
1494+
if (mcParticle.pt() < minPhiPt || mcParticle.pt() > maxPhiPt)
1495+
continue;
1496+
if (std::abs(mcParticle.y()) > cfgYAcceptance)
14711497
continue;
14721498

14731499
if (!isCountedMCPhi.at(0)) {
@@ -1512,6 +1538,8 @@ struct Phik0shortanalysis {
15121538

15131539
bool isCountedPhi = false;
15141540

1541+
double weight{1.0};
1542+
15151543
// Loop over all positive tracks
15161544
for (const auto& track1 : posThisColl) {
15171545
if (!selectionTrackResonance<true>(track1, true) || !selectionPIDKaonpTdependent(track1))
@@ -1529,17 +1557,20 @@ struct Phik0shortanalysis {
15291557
continue; // condition to avoid double counting of pair
15301558

15311559
ROOT::Math::PxPyPzMVector recPhi = recMother(track1, track2, massKa, massKa);
1532-
if (recPhi.Pt() < minPhiPt)
1560+
if (recPhi.Pt() < minPhiPt || recPhi.Pt() > maxPhiPt)
15331561
continue;
15341562
if (std::abs(recPhi.Rapidity()) > cfgYAcceptance)
15351563
continue;
15361564

15371565
if (!isCountedPhi) {
1538-
mcEventHist.fill(HIST("hRecMCEventSelection"), 7); // at least a Phi in the event
1566+
mcEventHist.fill(HIST("hRecMCEventSelection"), 7); // at least a Phi candidate in the event
15391567
mcEventHist.fill(HIST("hRecMCGenMultiplicityPercentWithPhi"), genmultiplicity);
15401568
isCountedPhi = true;
15411569
}
15421570

1571+
if (fillMethodSingleWeight)
1572+
weight *= (1 - getPhiPurity(genmultiplicity, recPhi));
1573+
15431574
closureMCPhiHist.fill(HIST("h3MCPhipurInvMass"), genmultiplicity, recPhi.Pt(), recPhi.M());
15441575

15451576
std::array<bool, 3> isCountedK0S{false, false, false};
@@ -1620,6 +1651,9 @@ struct Phik0shortanalysis {
16201651
}
16211652
}
16221653
}
1654+
1655+
weight = 1 - weight;
1656+
mcEventHist.fill(HIST("hRecMCEventSelection"), 8, weight); // at least a Phi in the event
16231657
}
16241658

16251659
PROCESS_SWITCH(Phik0shortanalysis, processRecMCClosurePhiQA, "Process for ReCMCQA and Phi in RecMCClosure", false);
@@ -1709,7 +1743,7 @@ struct Phik0shortanalysis {
17091743
}
17101744

17111745
ROOT::Math::PxPyPzMVector recPhi = recMother(track1, track2, massKa, massKa);
1712-
if (recPhi.Pt() < minPhiPt)
1746+
if (recPhi.Pt() < minPhiPt || recPhi.Pt() > maxPhiPt)
17131747
continue;
17141748
if (recPhi.M() < lowMPhi || recPhi.M() > upMPhi)
17151749
continue;
@@ -1833,7 +1867,7 @@ struct Phik0shortanalysis {
18331867
}
18341868

18351869
ROOT::Math::PxPyPzMVector recPhi = recMother(track1, track2, massKa, massKa);
1836-
if (recPhi.Pt() < minPhiPt)
1870+
if (recPhi.Pt() < minPhiPt || recPhi.Pt() > maxPhiPt)
18371871
continue;
18381872
if (recPhi.M() < lowMPhi || recPhi.M() > upMPhi)
18391873
continue;
@@ -2053,7 +2087,7 @@ struct Phik0shortanalysis {
20532087
if (!isPosKaon || !isNegKaon)
20542088
continue;
20552089
}
2056-
if (mcParticle2.pt() < minPhiPt)
2090+
if (mcParticle2.pt() < minPhiPt || mcParticle2.pt() > maxPhiPt)
20572091
continue;
20582092

20592093
if (std::abs(mcParticle2.y()) > cfgYAcceptance)
@@ -2138,7 +2172,7 @@ struct Phik0shortanalysis {
21382172
if (!isPosKaon || !isNegKaon)
21392173
continue;
21402174
}
2141-
if (mcParticle2.pt() < minPhiPt)
2175+
if (mcParticle2.pt() < minPhiPt || mcParticle2.pt() > maxPhiPt)
21422176
continue;
21432177

21442178
if (std::abs(mcParticle2.y()) > cfgYAcceptance)

0 commit comments

Comments
 (0)