Skip to content

Commit b9375c5

Browse files
scannitoscannito
authored
[PWGLF] Added process function for closure in MCGen (#15703)
1 parent 6547610 commit b9375c5

File tree

1 file changed

+246
-3
lines changed

1 file changed

+246
-3
lines changed

‎PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx‎

Lines changed: 246 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,7 @@
6363
#include <array>
6464
#include <cmath>
6565
#include <cstdlib>
66+
#include <deque>
6667
#include <memory>
6768
#include <ranges>
6869
#include <string>
@@ -157,7 +158,7 @@ struct PhiStrangenessCorrelation {
157158
Configurable<int> eventSelectionType{"eventSelectionType", 0, "Event selection type: 0 - default selection only, 1 - default + phi meson selection"};
158159

159160
// Configurable for analysis mode
160-
Configurable<int> analysisMode{"analysisMode", kMassvsMass, "Analysis mode: 0 - old method with online normalization, 1 - new method with offline normlization, 2 - deltay vs deltaphi"};
161+
Configurable<int> analysisMode{"analysisMode", kDeltaYvsDeltaPhi, "Analysis mode: 0 - old method with online normalization, 1 - new method with offline normlization, 2 - deltay vs deltaphi"};
161162

162163
// Configurable for event selection
163164
Configurable<float> cutZVertex{"cutZVertex", 10.0f, "Accepted z-vertex range (cm)"}; // TO BE REMOVED
@@ -173,13 +174,13 @@ struct PhiStrangenessCorrelation {
173174

174175
// Configurables for K0s selection
175176
struct : ConfigurableGroup {
176-
Configurable<bool> selectK0sInSigRegion{"selectK0sInSigRegion", false, "Select K0s candidates in signal region"};
177+
Configurable<bool> selectK0sInSigRegion{"selectK0sInSigRegion", true, "Select K0s candidates in signal region"};
177178
Configurable<std::pair<float, float>> rangeMK0sSignal{"rangeMK0sSignal", {0.47f, 0.53f}, "K0S mass range for signal extraction"};
178179
} k0sConfigs;
179180

180181
// Configurables for Pions selection
181182
struct : ConfigurableGroup {
182-
Configurable<bool> selectPionInSigRegion{"selectPionInSigRegion", false, "Select Pion candidates in signal region"};
183+
Configurable<bool> selectPionInSigRegion{"selectPionInSigRegion", true, "Select Pion candidates in signal region"};
183184
Configurable<float> pidTPCMax{"pidTPCMax", 2.0f, "Maximum nSigma TPC"};
184185
Configurable<float> pidTOFMax{"pidTOFMax", 2.0f, "Maximum nSigma TOF"};
185186
Configurable<float> tofPIDThreshold{"tofPIDThreshold", 0.5f, "Minimum pT after which TOF PID is applicable"};
@@ -260,6 +261,24 @@ struct PhiStrangenessCorrelation {
260261

261262
static constexpr std::array<std::string_view, 2> phiMassRegionLabels{"Signal", "Sideband"};
262263
static constexpr std::array<std::string_view, ParticleOfInterestSize> particleOfInterestLabels{"Phi", "K0S", "Pion" /*"PionTPC", "PionTPCTOF"*/};
264+
static constexpr std::array<std::string_view, ParticleOfInterestSize - 1> assocParticleLabels{"K0S", "Pi"};
265+
266+
// Light structures to store only the necessary information for the correlation analysis at MCGen level
267+
struct MiniParticle {
268+
float pt;
269+
float y;
270+
float phi;
271+
};
272+
273+
struct MiniEvent {
274+
float multiplicity;
275+
std::vector<MiniParticle> phiParticles;
276+
std::vector<MiniParticle> k0sParticles;
277+
std::vector<MiniParticle> pionParticles;
278+
};
279+
280+
// Buffer for mixed event, organized as a vector of deques, one for each multiplicity bin, containing the past events with their particles of interest needed for mixing
281+
std::vector<std::deque<MiniEvent>> eventBuffer;
263282

264283
void init(InitContext&)
265284
{
@@ -316,6 +335,12 @@ struct PhiStrangenessCorrelation {
316335
histos.add("pi/h3PiMCGen", "Pion in MC Gen", kTH3F, {binnedmultAxis, binnedpTPiAxis, yAxis});
317336
histos.add("pi/h4PiMCGenAssocReco", "Pion in MC Gen Assoc Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTPiAxis, yAxis});
318337

338+
histos.add("phiK0S/h5PhiK0SClosureMCGen", "Deltay vs deltaphi for Phi and K0Short in MCGen", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTK0SAxis, deltayAxis, deltaphiAxis});
339+
histos.add("phiPi/h5PhiPiClosureMCGen", "Deltay vs deltaphi for Phi and Pion in MCGen", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, deltaphiAxis});
340+
341+
histos.add("phiK0S/h5PhiK0SClosureMCGenME", "Deltay vs deltaphi for Phi and K0Short in MCGen ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTK0SAxis, deltayAxis, deltaphiAxis});
342+
histos.add("phiPi/h5PhiPiClosureMCGenME", "Deltay vs deltaphi for Phi and Pion in MCGen ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, deltaphiAxis});
343+
319344
// Load efficiency maps from CCDB
320345
if (applyEfficiency) {
321346
ccdb->setURL(ccdbUrl);
@@ -327,6 +352,8 @@ struct PhiStrangenessCorrelation {
327352
loadEfficiencyMapFromCCDB(static_cast<ParticleOfInterest>(i));
328353
}
329354
}
355+
356+
eventBuffer.resize(binsMult->size() - 1);
330357
}
331358

332359
void loadEfficiencyMapFromCCDB(ParticleOfInterest poi)
@@ -354,6 +381,15 @@ struct PhiStrangenessCorrelation {
354381
return RecoDecay::constrainAngle(phiTrigger - phiAssociated, -o2::constants::math::PIHalf);
355382
}
356383

384+
int getCentBin(float multiplicity)
385+
{
386+
if (multiplicity < binsMult->front() || multiplicity >= binsMult->back())
387+
return -1;
388+
389+
auto it = std::upper_bound(binsMult->begin(), binsMult->end(), multiplicity);
390+
return std::distance(binsMult->begin(), it) - 1;
391+
}
392+
357393
template <typename TCollision, typename TPhiCands, typename TK0SCands, typename TPionCands>
358394
void processPhiK0SPionSE(TCollision const& collision, TPhiCands const& phiCandidates, TK0SCands const& k0sReduced, TPionCands const& pionTracks)
359395
{
@@ -798,6 +834,213 @@ struct PhiStrangenessCorrelation {
798834
}
799835

800836
PROCESS_SWITCH(PhiStrangenessCorrelation, processParticleEfficiency, "Process function for Efficiency Computation for Particles of Interest", false);
837+
838+
/*void processMCGenClosureSE(MCCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles)
839+
{
840+
float multiplicity = mcCollision.centFT0M();
841+
842+
std::vector<int64_t> phiIndices;
843+
std::vector<int64_t> k0sIndices;
844+
std::vector<int64_t> pionIndices;
845+
std::vector<std::vector<int64_t>> assocIndices;
846+
847+
auto inYAcceptance = [&](const auto& mcParticle) {
848+
return std::abs(mcParticle.y()) <= yConfigs.cfgYAcceptance;
849+
};
850+
851+
for (const auto& mcParticle : mcParticles) {
852+
if (!inYAcceptance(mcParticle))
853+
continue;
854+
855+
switch (std::abs(mcParticle.pdgCode())) {
856+
case o2::constants::physics::Pdg::kPhi:
857+
if (eventSelectionType == 0 && mcParticle.pt() >= minPtMcGenConfigs.minPhiPt)
858+
phiIndices.push_back(mcParticle.globalIndex());
859+
break;
860+
case PDG_t::kK0Short:
861+
if (mcParticle.isPhysicalPrimary() && mcParticle.pt() >= minPtMcGenConfigs.v0SettingMinPt)
862+
k0sIndices.push_back(mcParticle.globalIndex());
863+
break;
864+
case PDG_t::kPiPlus:
865+
if (mcParticle.isPhysicalPrimary() && mcParticle.pt() >= minPtMcGenConfigs.cMinPionPtcut)
866+
pionIndices.push_back(mcParticle.globalIndex());
867+
break;
868+
default:
869+
break;
870+
}
871+
}
872+
873+
assocIndices.push_back(k0sIndices);
874+
assocIndices.push_back(pionIndices);
875+
876+
for (std::size_t iTrigg{0}; iTrigg < phiIndices.size(); ++iTrigg) {
877+
auto& phiParticle = mcParticles.rawIteratorAt(phiIndices[iTrigg]);
878+
879+
static_for<0, assocIndices.size() - 1>([&](auto i_idx) {
880+
constexpr unsigned int i = i_idx.value;
881+
882+
for (std::size_t iAssoc{0}; iAssoc < assocIndices[i].size(); ++iAssoc) {
883+
auto& assocParticle = mcParticles.rawIteratorAt(assocIndices[i][iAssoc]);
884+
885+
histos.fill(HIST("mcGenClosure/h5Phi") + HIST(assocParticleLabels[i]) + HIST("ClosureGenSE"), multiplicity, phiParticle.pt(), assocParticle.pt(), phiParticle.y() - assocParticle.y(), getDeltaPhi(phiParticle.phi(), assocParticle.phi()));
886+
}
887+
});
888+
}
889+
}
890+
891+
PROCESS_SWITCH(PhiStrangenessCorrelation, processMCGenClosureSE, "Process function for MC Gen Closure Test in SE", false);
892+
893+
void processMCGenClosureME(MCCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles)
894+
{
895+
float multiplicity = mcCollision.centFT0M();
896+
897+
std::vector<MiniParticle> phiParticles;
898+
std::vector<MiniParticle> k0sParticles;
899+
std::vector<MiniParticle> pionParticles;
900+
901+
auto inYAcceptance = [&](const auto& mcParticle) {
902+
return std::abs(mcParticle.y()) <= yConfigs.cfgYAcceptance;
903+
};
904+
905+
for (const auto& mcParticle : mcParticles) {
906+
if (!inYAcceptance(mcParticle))
907+
continue;
908+
909+
switch (std::abs(mcParticle.pdgCode())) {
910+
case o2::constants::physics::Pdg::kPhi:
911+
if (eventSelectionType == 0 && mcParticle.pt() >= minPtMcGenConfigs.minPhiPt)
912+
phiParticles.emplace_back(mcParticle.pt(), mcParticle.y(), mcParticle.phi());
913+
break;
914+
case PDG_t::kK0Short:
915+
if (mcParticle.isPhysicalPrimary() && mcParticle.pt() >= minPtMcGenConfigs.v0SettingMinPt)
916+
k0sParticles.emplace_back(mcParticle.pt(), mcParticle.y(), mcParticle.phi());
917+
break;
918+
case PDG_t::kPiPlus:
919+
if (mcParticle.isPhysicalPrimary() && mcParticle.pt() >= minPtMcGenConfigs.cMinPionPtcut)
920+
pionParticles.emplace_back(mcParticle.pt(), mcParticle.y(), mcParticle.phi());
921+
break;
922+
default:
923+
break;
924+
}
925+
}
926+
927+
if (phiParticles.empty() && k0sParticles.empty() && pionParticles.empty())
928+
return;
929+
930+
int multBin = getCentBin(multiplicity);
931+
932+
// Loop over past events in the same multiplicity bin and fill histograms with all combinations of current phi particles and past K0S and pion particles
933+
for (const auto& pastEvent : eventBuffer[multBin]) {
934+
for (const auto& phiParticle : phiParticles) {
935+
for (const auto& k0sParticle : pastEvent.k0sParticles) {
936+
histos.fill(HIST("mcGenClosure/h5PhiK0SClosureGenME"), multiplicity, phiParticle.pt, k0sParticle.pt, phiParticle.y - k0sParticle.y, getDeltaPhi(phiParticle.phi, k0sParticle.phi));
937+
}
938+
for (const auto& pionParticle : pastEvent.pionParticles) {
939+
histos.fill(HIST("mcGenClosure/h5PhiPiClosureGenME"), multiplicity, phiParticle.pt, pionParticle.pt, phiParticle.y - pionParticle.y, getDeltaPhi(phiParticle.phi, pionParticle.phi));
940+
}
941+
}
942+
}
943+
944+
// Add current event to buffer
945+
MiniEvent currentEvent;
946+
currentEvent.multiplicity = multiplicity;
947+
currentEvent.phiParticles = std::move(phiParticles);
948+
currentEvent.k0sParticles = std::move(k0sParticles);
949+
currentEvent.pionParticles = std::move(pionParticles);
950+
951+
eventBuffer[multBin].push_front(std::move(currentEvent));
952+
if (eventBuffer[multBin].size() > static_cast<std::size_t>(cfgNoMixedEvents.value))
953+
eventBuffer[multBin].pop_back();
954+
}
955+
956+
PROCESS_SWITCH(PhiStrangenessCorrelation, processMCGenClosureME, "Process function for MC Gen Closure Test in ME", false);*/
957+
958+
void processMCGenClosure(MCCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles)
959+
{
960+
float multiplicity = mcCollision.centFT0M();
961+
962+
std::vector<MiniParticle> phiParticles;
963+
std::vector<MiniParticle> k0sParticles;
964+
std::vector<MiniParticle> pionParticles;
965+
966+
auto inYAcceptance = [&](const auto& mcParticle) {
967+
return std::abs(mcParticle.y()) <= yConfigs.cfgYAcceptance;
968+
};
969+
970+
for (const auto& mcParticle : mcParticles) {
971+
if (!inYAcceptance(mcParticle))
972+
continue;
973+
974+
switch (std::abs(mcParticle.pdgCode())) {
975+
case o2::constants::physics::Pdg::kPhi:
976+
if (eventSelectionType == 0 && mcParticle.pt() >= minPtMcGenConfigs.minPhiPt)
977+
phiParticles.emplace_back(mcParticle.pt(), mcParticle.y(), mcParticle.phi());
978+
break;
979+
case PDG_t::kK0Short:
980+
if (mcParticle.isPhysicalPrimary() && mcParticle.pt() >= minPtMcGenConfigs.v0SettingMinPt)
981+
k0sParticles.emplace_back(mcParticle.pt(), mcParticle.y(), mcParticle.phi());
982+
break;
983+
case PDG_t::kPiPlus:
984+
if (mcParticle.isPhysicalPrimary() && mcParticle.pt() >= minPtMcGenConfigs.cMinPionPtcut)
985+
pionParticles.emplace_back(mcParticle.pt(), mcParticle.y(), mcParticle.phi());
986+
break;
987+
default:
988+
break;
989+
}
990+
}
991+
992+
if (phiParticles.empty() && k0sParticles.empty() && pionParticles.empty())
993+
return;
994+
995+
int multBin = getCentBin(multiplicity);
996+
if (multBin < 0)
997+
return;
998+
999+
std::vector<MiniParticle>* currentAssocParticles[] = {&k0sParticles, &pionParticles};
1000+
1001+
for (const auto& phiParticle : phiParticles) {
1002+
static_for<0, assocParticleLabels.size() - 1>([&](auto i_idx) {
1003+
constexpr unsigned int i = i_idx.value;
1004+
1005+
for (const auto& assocParticle : *(currentAssocParticles[i])) {
1006+
histos.fill(HIST("phi") + HIST(assocParticleLabels[i]) + HIST("/h5Phi") + HIST(assocParticleLabels[i]) + HIST("ClosureMCGen"),
1007+
multiplicity, phiParticle.pt, assocParticle.pt,
1008+
phiParticle.y - assocParticle.y,
1009+
getDeltaPhi(phiParticle.phi, assocParticle.phi));
1010+
}
1011+
});
1012+
}
1013+
1014+
for (const auto& pastEvent : eventBuffer[multBin]) {
1015+
const std::vector<MiniParticle>* pastAssocParticles[] = {&pastEvent.k0sParticles, &pastEvent.pionParticles};
1016+
1017+
// Loop over past events in the same multiplicity bin and fill histograms with all combinations of current phi particles and past associated particles
1018+
for (const auto& phiParticle : phiParticles) {
1019+
static_for<0, assocParticleLabels.size() - 1>([&](auto i_idx) {
1020+
constexpr unsigned int i = i_idx.value;
1021+
1022+
for (const auto& assocParticle : *(pastAssocParticles[i])) {
1023+
histos.fill(HIST("phi") + HIST(assocParticleLabels[i]) + HIST("/h5Phi") + HIST(assocParticleLabels[i]) + HIST("ClosureMCGenME"),
1024+
multiplicity, phiParticle.pt, assocParticle.pt,
1025+
phiParticle.y - assocParticle.y,
1026+
getDeltaPhi(phiParticle.phi, assocParticle.phi));
1027+
}
1028+
});
1029+
}
1030+
}
1031+
1032+
MiniEvent currentEvent;
1033+
currentEvent.multiplicity = multiplicity;
1034+
currentEvent.phiParticles = std::move(phiParticles);
1035+
currentEvent.k0sParticles = std::move(k0sParticles);
1036+
currentEvent.pionParticles = std::move(pionParticles);
1037+
1038+
eventBuffer[multBin].push_front(std::move(currentEvent));
1039+
if (eventBuffer[multBin].size() > static_cast<std::size_t>(cfgNoMixedEvents.value))
1040+
eventBuffer[multBin].pop_back();
1041+
}
1042+
1043+
PROCESS_SWITCH(PhiStrangenessCorrelation, processMCGenClosure, "Process function for MC Gen Closure Test in SE and ME", false);
8011044
};
8021045

8031046
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)