Skip to content

[PWGLF] generatorsID in UPC task #10537

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ChiaraDeMartin95 merged 4 commits into from
Mar 18, 2025
Merged

[PWGLF] generatorsID in UPC task #10537

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6df0893
generatorsID
nepeivodaRS Mar 17, 2025
cc287dc
fix pdg codes
nepeivodaRS Mar 17, 2025
49967a6
fix pdg codes 2
nepeivodaRS Mar 17, 2025
7e2007f
remove unused vars
nepeivodaRS Mar 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
210 changes: 176 additions & 34 deletions PWGLF/Tasks/Strangeness/derivedupcanalysis.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,8 @@ struct Derivedupcanalysis {
Configurable<bool> analyseOmega{"analyseOmega", true, "process Omega-like candidates"};
Configurable<bool> analyseAntiOmega{"analyseAntiOmega", true, "process AntiOmega-like candidates"};

Configurable<std::vector<int>> generatorIds{"generatorIds", std::vector<int>{-1}, "MC generatorIds to process"};

// Event selections
Configurable<bool> rejectITSROFBorder{"rejectITSROFBorder", true, "reject events at ITS ROF border"};
Configurable<bool> rejectTFBorder{"rejectTFBorder", true, "reject events at TF border"};
Expand Down Expand Up @@ -861,6 +863,7 @@ struct Derivedupcanalysis {
histos.add("eventQA/hTracksPVeta1VsTracksGlobal", "hTracksPVeta1VsTracksGlobal", kTH3F, {axisNTracksPVeta1, axisNTracksGlobal, axisSelGap});
histos.add("eventQA/hCentralityVsTracksGlobal", "hCentralityVsTracksGlobal", kTH3F, {axisFT0Cqa, axisNTracksGlobal, axisSelGap});
histos.add("eventQA/hGapSide", "Gap side; Entries", kTH1F, {{5, -0.5, 4.5}});
histos.add("eventQA/hRawGapSide", "Raw Gap side; Entries", kTH1F, {{6, -1.5, 4.5}});
histos.add("eventQA/hSelGapSide", "Selected gap side; Entries", kTH1F, {axisSelGap});
histos.add("eventQA/hPosX", "Vertex position in x", kTH2F, {{100, -0.1, 0.1}, axisSelGap});
histos.add("eventQA/hPosY", "Vertex position in y", kTH2F, {{100, -0.1, 0.1}, axisSelGap});
Expand Down Expand Up @@ -977,14 +980,16 @@ struct Derivedupcanalysis {
}

template <typename TCollision>
int getGapSide(TCollision const& collision)
int getGapSide(TCollision const& collision, bool fillQA)
{
int selGapSide = sgSelector.trueGap(collision, upcCuts.fv0a, upcCuts.ft0a, upcCuts.ft0c, upcCuts.zdc);
histos.fill(HIST("eventQA/hGapSide"), collision.gapSide());
histos.fill(HIST("eventQA/hSelGapSide"), selGapSide);
histos.fill(HIST("eventQA/hFT0"), collision.totalFT0AmplitudeA(), collision.totalFT0AmplitudeC(), selGapSide);
histos.fill(HIST("eventQA/hFDD"), collision.totalFDDAmplitudeA(), collision.totalFDDAmplitudeC(), selGapSide);
histos.fill(HIST("eventQA/hZN"), collision.energyCommonZNA(), collision.energyCommonZNC(), selGapSide);
if (fillQA) {
histos.fill(HIST("eventQA/hGapSide"), collision.gapSide());
histos.fill(HIST("eventQA/hSelGapSide"), selGapSide);
histos.fill(HIST("eventQA/hFT0"), collision.totalFT0AmplitudeA(), collision.totalFT0AmplitudeC(), selGapSide);
histos.fill(HIST("eventQA/hFDD"), collision.totalFDDAmplitudeA(), collision.totalFDDAmplitudeC(), selGapSide);
histos.fill(HIST("eventQA/hZN"), collision.energyCommonZNA(), collision.energyCommonZNC(), selGapSide);
}
return selGapSide;
}

Expand Down Expand Up @@ -1014,7 +1019,7 @@ struct Derivedupcanalysis {
float qaBin;
};

const std::array<SelectionCheck, 15> checks = {{
const std::array<SelectionCheck, 14> checks = {{
{true, true, 0.0}, // All collisions
{requireIsTriggerTVX, collision.selection_bit(aod::evsel::kIsTriggerTVX), 1.0}, // Triggered by FT0M
{true, std::fabs(collision.posZ()) < maxZVtxPosition, 2.0}, // Vertex-Z selected
Expand All @@ -1029,7 +1034,6 @@ struct Derivedupcanalysis {
{requireNoCollInTimeRangeNarrow, collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeNarrow), 11.0}, // No collision within +-4 µs
{minOccupancy >= 0, collision.trackOccupancyInTimeRange() >= minOccupancy, 12.0}, // Above min occupancy
{maxOccupancy > 0, collision.trackOccupancyInTimeRange() < maxOccupancy, 13.0}, // Below max occupancy
{studyUPConly, collision.isUPC(), 14.0}, // Study UPC collisions only
}};

for (const auto& check : checks) {
Expand All @@ -1041,6 +1045,12 @@ struct Derivedupcanalysis {
}
}

if (studyUPConly && !collision.isUPC()) {
return false;
} else if (collision.isUPC()) {
histos.fill(HIST("eventQA/hEventSelection"), 14.0); // is UPC compatible
}

// Additional check for UPC collision flag
if (useUPCflag && collision.flags() < 1) {
return false;
Expand Down Expand Up @@ -1511,27 +1521,27 @@ struct Derivedupcanalysis {
const int pdgV0 = v0.pdgCode();
const bool isPhysPrim = v0.isPhysicalPrimary();

const bool isPositiveProton = (pdgPos == 2212);
const bool isPositivePion = (pdgPos == 211) || (doTreatPiToMuon && pdgPos == -13);
const bool isNegativeProton = (pdgNeg == -2212);
const bool isNegativePion = (pdgNeg == -211) || (doTreatPiToMuon && pdgNeg == 13);
const bool isPositiveProton = (pdgPos == PDG_t::kProton);
const bool isPositivePion = (pdgPos == PDG_t::kPiPlus) || (doTreatPiToMuon && pdgPos == PDG_t::kMuonPlus);
const bool isNegativeProton = (pdgNeg == kProtonBar);
const bool isNegativePion = (pdgNeg == PDG_t::kPiMinus) || (doTreatPiToMuon && pdgNeg == PDG_t::kMuonMinus);

switch (pdgV0) {
case 310: // K0Short
case PDG_t::kK0Short: // K0Short
if (isPositivePion && isNegativePion) {
bitMap.set(selConsiderK0Short);
if (isPhysPrim)
bitMap.set(selPhysPrimK0Short);
}
break;
case 3122: // Lambda
case PDG_t::kLambda0: // Lambda
if (isPositiveProton && isNegativePion) {
bitMap.set(selConsiderLambda);
if (isPhysPrim)
bitMap.set(selPhysPrimLambda);
}
break;
case -3122: // AntiLambda
case PDG_t::kLambda0Bar: // AntiLambda
if (isPositivePion && isNegativeProton) {
bitMap.set(selConsiderAntiLambda);
if (isPhysPrim)
Expand All @@ -1541,6 +1551,60 @@ struct Derivedupcanalysis {
}
}

template <typename TCasc>
void computeCascadeMCAssociation(const TCasc& casc, std::bitset<kSelNum>& bitMap)
{
const int pdgPos = casc.pdgCodePositive();
const int pdgNeg = casc.pdgCodeNegative();
const int pdgBach = casc.pdgCodeBachelor();
const int pdgCasc = casc.pdgCode();
const bool isPhysPrim = casc.isPhysicalPrimary();

const bool isPositiveProton = (pdgPos == PDG_t::kProton);

const bool isPositivePion = (pdgPos == PDG_t::kPiPlus);
const bool isBachelorPositivePion = (pdgBach == PDG_t::kPiPlus);

const bool isNegativeProton = (pdgNeg == kProtonBar);

const bool isNegativePion = (pdgNeg == PDG_t::kPiMinus);
const bool isBachelorNegativePion = (pdgBach == PDG_t::kPiMinus);

const bool isBachelorPositiveKaon = (pdgBach == PDG_t::kKPlus);
const bool isBachelorNegativeKaon = (pdgBach == PDG_t::kKMinus);

switch (pdgCasc) {
case PDG_t::kXiMinus: // Xi
if (isPositiveProton && isNegativePion && isBachelorNegativePion) {
bitMap.set(selConsiderXi);
if (isPhysPrim)
bitMap.set(selPhysPrimXi);
}
break;
case PDG_t::kXiPlusBar: // Anti-Xi
if (isNegativeProton && isPositivePion && isBachelorPositivePion) {
bitMap.set(selConsiderAntiXi);
if (isPhysPrim)
bitMap.set(selPhysPrimAntiXi);
}
break;
case PDG_t::kOmegaMinus: // Omega
if (isPositiveProton && isNegativePion && isBachelorNegativeKaon) {
bitMap.set(selConsiderOmega);
if (isPhysPrim)
bitMap.set(selPhysPrimOmega);
}
break;
case PDG_t::kOmegaPlusBar: // Anti-Omega
if (isNegativeProton && isPositivePion && isBachelorPositiveKaon) {
bitMap.set(selConsiderAntiOmega);
if (isPhysPrim)
bitMap.set(selPhysPrimAntiOmega);
}
break;
}
}

template <typename TV0, typename TCollision>
void analyseV0Candidate(TV0 const& v0, TCollision const& coll, int const& gap, std::bitset<kSelNum> const& selMap)
{
Expand Down Expand Up @@ -1585,6 +1649,9 @@ struct Derivedupcanalysis {
std::vector<int> listBestCollisionIds(mcCollisions.size(), -1);

for (auto const& mcCollision : mcCollisions) {
if (std::find(generatorIds->begin(), generatorIds->end(), mcCollision.generatorsID()) == generatorIds->end()) {
continue;
}
auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex());
// Find the collision with the biggest nbr of PV contributors
// Follows what was done here: https://github.com/AliceO2Group/O2Physics/blob/master/Common/TableProducer/mcCollsExtra.cxx#L93
Expand All @@ -1595,7 +1662,7 @@ struct Derivedupcanalysis {
continue;
}

int selGapSide = collision.isUPC() ? getGapSide(collision) : -1;
int selGapSide = collision.isUPC() ? getGapSide(collision, false) : -1;
if (studyUPConly && (selGapSide < -0.5))
continue;

Expand All @@ -1613,6 +1680,11 @@ struct Derivedupcanalysis {
void fillGenMCHistogramsQA(StraMCCollisionsFull const& mcCollisions, StraCollisonsFullMC const& collisions)
{
for (auto const& mcCollision : mcCollisions) {
// LOGF(info, "Generator ID is %i", mcCollision.generatorsID());
if (std::find(generatorIds->begin(), generatorIds->end(), mcCollision.generatorsID()) == generatorIds->end()) {
continue;
}

histos.fill(HIST("eventQA/mc/hEventSelectionMC"), 0.0);
histos.fill(HIST("eventQA/mc/hMCNParticlesEta10"), mcCollision.multMCNParticlesEta10(), 0 /* all gen. events*/);

Expand All @@ -1636,7 +1708,7 @@ struct Derivedupcanalysis {
continue;
}

int selGapSide = collision.isUPC() ? getGapSide(collision) : -1;
int selGapSide = collision.isUPC() ? getGapSide(collision, false) : -1;
if (studyUPConly && (selGapSide < -0.5))
continue;

Expand Down Expand Up @@ -1670,7 +1742,9 @@ struct Derivedupcanalysis {
return;
} // event is accepted

int selGapSide = collision.isUPC() ? getGapSide(collision) : -1;
histos.fill(HIST("eventQA/hRawGapSide"), collision.gapSide());

int selGapSide = collision.isUPC() ? getGapSide(collision, true) : -1;
if (studyUPConly && (selGapSide < -0.5))
return;

Expand All @@ -1697,15 +1771,22 @@ struct Derivedupcanalysis {
StraMCCollisionsFull const&,
V0MCCoresFull const&)
{
if (!collision.has_straMCCollision()) {
histos.fill(HIST("eventQA/mc/hFakeEvents"), 0); // no assoc. MC collisions
} else {
const auto& mcCollision = collision.straMCCollision_as<StraMCCollisionsFull>();
if (std::find(generatorIds->begin(), generatorIds->end(), mcCollision.generatorsID()) == generatorIds->end()) {
return;
}
}

if (!acceptEvent(collision, true)) {
return;
} // event is accepted

if (!collision.has_straMCCollision()) {
histos.fill(HIST("eventQA/mc/hFakeEvents"), 0); // no assoc. MC collisions
}
histos.fill(HIST("eventQA/hRawGapSide"), collision.gapSide());

int selGapSide = collision.isUPC() ? getGapSide(collision) : -1;
int selGapSide = collision.isUPC() ? getGapSide(collision, true) : -1;
if (studyUPConly && (selGapSide < -0.5))
return;

Expand Down Expand Up @@ -1748,7 +1829,9 @@ struct Derivedupcanalysis {
return;
} // event is accepted

int selGapSide = collision.isUPC() ? getGapSide(collision) : -1;
histos.fill(HIST("eventQA/hRawGapSide"), collision.gapSide());

int selGapSide = collision.isUPC() ? getGapSide(collision, true) : -1;
if (studyUPConly && (selGapSide < -0.5))
return;

Expand All @@ -1764,6 +1847,60 @@ struct Derivedupcanalysis {
} // end casc loop
}

void processCascadesMC(StraCollisonFullMC const& collision,
CascadeCandidatesMC const& fullCascades,
DauTracks const&,
aod::MotherMCParts const&,
StraMCCollisionsFull const&,
CascMCCoresFull const&)
{
if (!collision.has_straMCCollision()) {
histos.fill(HIST("eventQA/mc/hFakeEvents"), 0); // no assoc. MC collisions
} else {
const auto& mcCollision = collision.straMCCollision_as<StraMCCollisionsFull>();
if (std::find(generatorIds->begin(), generatorIds->end(), mcCollision.generatorsID()) == generatorIds->end()) {
return;
}
}

if (!acceptEvent(collision, true)) {
return;
} // event is accepted

histos.fill(HIST("eventQA/hRawGapSide"), collision.gapSide());

int selGapSide = collision.isUPC() ? getGapSide(collision, true) : -1;
if (studyUPConly && (selGapSide < -0.5))
return;

fillHistogramsQA(collision, selGapSide);

if (collision.has_straMCCollision()) {
const auto& mcCollision = collision.straMCCollision_as<StraMCCollisionsFull>();
histos.fill(HIST("eventQA/mc/hNTracksGlobalvsMCNParticlesEta10rec"), collision.multNTracksGlobal(), mcCollision.multMCNParticlesEta10());
histos.fill(HIST("eventQA/mc/hNTracksPVeta1vsMCNParticlesEta10rec"), collision.multNTracksPVeta1(), mcCollision.multMCNParticlesEta10());
histos.fill(HIST("eventQA/mc/hNTracksGlobalvstotalMultMCParticles"), collision.multNTracksGlobal(), mcCollision.totalMultMCParticles());
histos.fill(HIST("eventQA/mc/hNTracksPVeta1vstotalMultMCParticles"), collision.multNTracksPVeta1(), mcCollision.totalMultMCParticles());
}

for (const auto& casc : fullCascades) {
std::bitset<kSelNum> selMap = computeBitmapCascade(casc, collision);

if (doMCAssociation) {
if (casc.has_cascMCCore()) {
const auto& cascMC = casc.cascMCCore_as<CascMCCoresFull>();
computeCascadeMCAssociation(cascMC, selMap);
}
} else {
// the candidate may belong to any particle species
setBits(selMap, {selConsiderXi, selConsiderAntiXi, selConsiderOmega, selConsiderAntiOmega,
selPhysPrimXi, selPhysPrimAntiXi, selPhysPrimOmega, selPhysPrimAntiOmega});
}

analyseCascCandidate(casc, collision, selGapSide, selMap);
} // end casc loop
}

void processGenerated(StraMCCollisionsFull const& mcCollisions,
V0MCCoresFull const& V0MCCores,
CascMCCoresFull const& CascMCCores,
Expand All @@ -1780,9 +1917,9 @@ struct Derivedupcanalysis {
// Kinematics (|y| < rapidityCut)
float pTmc = v0MC.ptMC();
float ymc = 1e3;
if (v0MC.pdgCode() == 310)
if (v0MC.pdgCode() == PDG_t::kK0Short)
ymc = v0MC.rapidityMC(0);
else if (std::abs(v0MC.pdgCode()) == 3122)
else if (std::abs(v0MC.pdgCode()) == PDG_t::kLambda0)
ymc = v0MC.rapidityMC(1);
if (std::abs(ymc) > rapidityCut)
continue;
Expand All @@ -1792,6 +1929,11 @@ struct Derivedupcanalysis {
if (std::abs(mcCollision.posZ()) > maxZVtxPosition)
continue;

// Collision is of the proccess of interest
if (std::find(generatorIds->begin(), generatorIds->end(), mcCollision.generatorsID()) == generatorIds->end()) {
continue;
}

float centrality = -1.f;
int nTracksGlobal = -1;

Expand All @@ -1802,13 +1944,13 @@ struct Derivedupcanalysis {
}

// Fill histograms
if (v0MC.pdgCode() == 310) {
if (v0MC.pdgCode() == PDG_t::kK0Short) {
histos.fill(HIST(kParticlenames[0]) + HIST("/mc/h6dGen"), centrality, nTracksGlobal, mcCollision.multMCNParticlesEta10(), pTmc, static_cast<int>(upcCuts.genGapSide), ymc);
}
if (v0MC.pdgCode() == 3122) {
if (v0MC.pdgCode() == PDG_t::kLambda0) {
histos.fill(HIST(kParticlenames[1]) + HIST("/mc/h6dGen"), centrality, nTracksGlobal, mcCollision.multMCNParticlesEta10(), pTmc, static_cast<int>(upcCuts.genGapSide), ymc);
}
if (v0MC.pdgCode() == -3122) {
if (v0MC.pdgCode() == PDG_t::kLambda0Bar) {
histos.fill(HIST(kParticlenames[2]) + HIST("/mc/h6dGen"), centrality, nTracksGlobal, mcCollision.multMCNParticlesEta10(), pTmc, static_cast<int>(upcCuts.genGapSide), ymc);
}
} // V0 end
Expand All @@ -1821,9 +1963,9 @@ struct Derivedupcanalysis {
// Kinematics (|y| < rapidityCut)
float pTmc = cascMC.ptMC();
float ymc = 1e3;
if (std::abs(cascMC.pdgCode()) == 3312)
if (std::abs(cascMC.pdgCode()) == PDG_t::kXiMinus)
ymc = RecoDecay::y(std::array{cascMC.pxMC(), cascMC.pyMC(), cascMC.pzMC()}, o2::constants::physics::MassXiMinus);
else if (std::abs(cascMC.pdgCode()) == 3334)
else if (std::abs(cascMC.pdgCode()) == PDG_t::kOmegaMinus)
ymc = RecoDecay::y(std::array{cascMC.pxMC(), cascMC.pyMC(), cascMC.pzMC()}, o2::constants::physics::MassOmegaMinus);
if (std::abs(ymc) > rapidityCut)
continue;
Expand All @@ -1842,16 +1984,16 @@ struct Derivedupcanalysis {
}

// Fill histograms
if (cascMC.pdgCode() == 3312) {
if (cascMC.pdgCode() == PDG_t::kXiMinus) {
histos.fill(HIST(kParticlenames[3]) + HIST("/mc/h6dGen"), centrality, nTracksGlobal, mcCollision.multMCNParticlesEta10(), pTmc, static_cast<int>(upcCuts.genGapSide), ymc);
}
if (cascMC.pdgCode() == -3312) {
if (cascMC.pdgCode() == PDG_t::kXiPlusBar) {
histos.fill(HIST(kParticlenames[4]) + HIST("/mc/h6dGen"), centrality, nTracksGlobal, mcCollision.multMCNParticlesEta10(), pTmc, static_cast<int>(upcCuts.genGapSide), ymc);
}
if (cascMC.pdgCode() == 3334) {
if (cascMC.pdgCode() == PDG_t::kOmegaMinus) {
histos.fill(HIST(kParticlenames[5]) + HIST("/mc/h6dGen"), centrality, nTracksGlobal, mcCollision.multMCNParticlesEta10(), pTmc, static_cast<int>(upcCuts.genGapSide), ymc);
}
if (cascMC.pdgCode() == -3334) {
if (cascMC.pdgCode() == PDG_t::kOmegaPlusBar) {
histos.fill(HIST(kParticlenames[6]) + HIST("/mc/h6dGen"), centrality, nTracksGlobal, mcCollision.multMCNParticlesEta10(), pTmc, static_cast<int>(upcCuts.genGapSide), ymc);
}
} // Cascade end
Expand Down
Loading