Skip to content

[PWGCF] FemtoUniverse -- Adding correction on secondary contaminations for the p-phi analysis #11955

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
zchochul merged 2 commits into from
Jul 7, 2025
Merged

[PWGCF] FemtoUniverse -- Adding correction on secondary contaminations for the p-phi analysis #11955

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ef8d9ee
adding secondary contamination calculations for p-phi analysis
Jul 7, 2025
754df18
removing commited sections
Jul 7, 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
5 changes: 4 additions & 1 deletion PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@
framework::Configurable<std::string> confEffCorCCDBPath{"confEffCorCCDBPath", "", "[Efficiency Correction] CCDB path to histograms"};
framework::Configurable<std::vector<std::string>> confEffCorCCDBTimestamps{"confEffCorCCDBTimestamps", {}, "[Efficiency Correction] Timestamps of histograms in CCDB (0 can be used as a placeholder, e.g. when running subwagons)"};
framework::Configurable<std::string> confEffCorVariables{"confEffCorVariables", "pt", "[Efficiency Correction] Variables for efficiency correction histogram dimensions (available: 'pt'; 'pt,eta'; 'pt,mult'; 'pt,eta,mult')"};
framework::Configurable<bool> confEffCorSetMultToConst{"confEffCorSetMultToConst", false, "[Efficiency Correction] Multiplicity for the histograms set to the constant value"};
};

class EfficiencyCorrection
Expand All @@ -62,6 +63,7 @@
auto init(framework::HistogramRegistry* registry, std::vector<framework::AxisSpec> axisSpecs) -> void
{
shouldFillHistograms = config->confEffCorFillHist;
shouldSetMultToConst = config->confEffCorSetMultToConst;

histRegistry = registry;
if (shouldFillHistograms) {
Expand Down Expand Up @@ -126,7 +128,7 @@
histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hMCTruth"),
particle.pt(),
particle.eta(),
particle.template fdCollision_as<CollisionType>().multV0M());
shouldSetMultToConst ? 100 : particle.template fdCollision_as<CollisionType>().multV0M());
}

template <uint8_t N, typename CollisionType = aod::FdCollisions>
Expand Down Expand Up @@ -277,13 +279,14 @@

bool shouldApplyCorrection{false};
bool shouldFillHistograms{false};
bool shouldSetMultToConst{false};

o2::ccdb::BasicCCDBManager& ccdb{o2::ccdb::BasicCCDBManager::instance()};
std::array<TH1*, 2> hLoaded{nullptr, nullptr};

framework::HistogramRegistry* histRegistry{};
static constexpr std::string_view histDirectory{"EfficiencyCorrection"};

Check failure on line 288 in PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/constexpr-constant] 

Use UpperCamelCase for names of constexpr constants. Names of special constants may be prefixed with "k".
static constexpr std::string_view histSuffix[2]{"one", "two"};

Check failure on line 289 in PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/constexpr-constant] 

Use UpperCamelCase for names of constexpr constants. Names of special constants may be prefixed with "k".
};

} // namespace o2::analysis::femto_universe::efficiency_correction
Expand Down
115 changes: 70 additions & 45 deletions PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackPhi.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,29 +16,31 @@
/// \author Anton Riedel, TU München, anton.riedel@tum.de
/// \author Zuzanna Chochulska, WUT Warsaw & CTU Prague, zchochul@cern.ch

#include <vector>
#include <string>
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseContainer.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseDetaDphiStar.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseEventHisto.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniversePairCleaner.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseParticleHisto.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseTrackSelection.h"
#include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h"

#include "Framework/ASoAHelpers.h"
#include "Framework/AnalysisTask.h"
#include "Framework/runDataProcessing.h"
#include "Framework/HistogramRegistry.h"
#include "Framework/ASoAHelpers.h"
#include "Framework/O2DatabasePDGPlugin.h"
#include "Framework/runDataProcessing.h"
#include "ReconstructionDataFormats/PID.h"

#include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseParticleHisto.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseEventHisto.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniversePairCleaner.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseContainer.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseDetaDphiStar.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseTrackSelection.h"
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCalculator.h"
#include <TFile.h>
#include <TH1.h>

#include <string>
#include <vector>

using namespace o2;
using namespace o2::analysis::femto_universe;
using namespace o2::analysis::femto_universe::efficiency;
using namespace o2::analysis::femto_universe::efficiency_correction;
using namespace o2::framework;
using namespace o2::framework::expressions;
using namespace o2::soa;
Expand Down Expand Up @@ -187,8 +189,10 @@

ColumnBinningPolicy<aod::collision::PosZ, aod::femtouniversecollision::MultNtr> colBinning{{ConfBinsVtx, ConfBinsMult}, true};

EfficiencyConfigurableGroup effConfGroup;
EfficiencyCalculator<TH1> efficiencyCalculator{&effConfGroup};
HistogramRegistry effCorrRegistry{"EfficiencyCorrection", {}, OutputObjHandlingPolicy::AnalysisObject};

EffCorConfigurableGroup effCorConfGroup;
EfficiencyCorrection effCorrection{&effCorConfGroup};

float weight = 1;

Expand Down Expand Up @@ -372,6 +376,12 @@
}
}

/// @returns 1 if positive, -1 if negative, 0 if zero
auto sign(auto number) -> int8_t
{
return (number > 0) - (number < 0);
}

void init(InitContext&)
{
if (ConfIsMC) {
Expand All @@ -383,7 +393,14 @@
registryMCpT.add("MCReco/C_p_pT", "; #it{p_T} (GeV/#it{c}); Counts", kTH1F, {{100, 0, 10}});
registryMCpT.add("MCReco/NC_p_pT", "; #it{p_T} (GeV/#it{c}); Counts", kTH1F, {{100, 0, 10}});
}
efficiencyCalculator.init();

effCorrection.init(
&effCorrRegistry,
{
static_cast<framework::AxisSpec>(ConfBinsTempFitVarpT),
{ConfBinsEta, -1, 1},
ConfBinsMult,
});

eventHisto.init(&qaRegistry);
qaRegistry.add("PhiDaugh_pos/nSigmaTPC", "; #it{p} (GeV/#it{c}); n#sigma_{TPC}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}});
Expand Down Expand Up @@ -467,7 +484,7 @@
}
}

template <typename PartitionType, typename PartType, typename MCParticles = std::nullptr_t>
template <bool isMC, typename PartitionType, typename PartType, typename MCParticles = std::nullptr_t>
void doSameEvent(PartitionType groupPartsTrack, PartitionType groupPartsPhi, PartType parts, float magFieldTesla, int multCol, [[maybe_unused]] MCParticles mcParts = nullptr)
{
for (auto const& phicandidate : groupPartsPhi) {
Expand All @@ -493,6 +510,17 @@
qaRegistry.fill(HIST("PhiDaugh_neg/hDCAxy"), negChild.p(), negChild.tempFitVar());

trackHistoPartPhi.fillQA<false, false>(phicandidate);
if constexpr (isMC) {
// reco
effCorrection.fillRecoHist<ParticleNo::ONE, FilteredFDCollisions>(phicandidate, 333);
// truth
auto mcPartId1 = phicandidate.fdMCParticleId();
auto const& mcpart1 = mcParts.iteratorAt(mcPartId1);
if (mcpart1.pdgMCTruth() != 333) {

Check failure on line 519 in PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackPhi.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
continue;
}
effCorrection.fillTruthHist<ParticleNo::ONE, FilteredFDCollisions>(phicandidate);
}
}

for (auto const& track : groupPartsTrack) {
Expand Down Expand Up @@ -529,6 +557,18 @@
qaRegistry.fill(HIST("Hadron_neg/nSigmaTOFPr"), track.p(), tofNSigmaPr);
}
trackHistoPartTrack.fillQA<false, false>(track);

if constexpr (isMC) {
effCorrection.fillRecoHist<ParticleNo::TWO, FilteredFDCollisions>(track, ConfTrackPDGCode);

// truth
auto mcPartId2 = track.fdMCParticleId();
auto const& mcpart2 = mcParts.iteratorAt(mcPartId2);
if (mcpart2.pdgMCTruth() != ConfTrackPDGCode) {
continue;
}
effCorrection.fillTruthHist<ParticleNo::TWO, FilteredFDCollisions>(track);
}
}

/// Now build the combinations
Expand Down Expand Up @@ -556,13 +596,8 @@
if (!pairCleaner.isCleanPair(track, phicandidate, parts)) {
continue;
}

weight = efficiencyCalculator.getWeight(ParticleNo::ONE, phicandidate.pt()) * efficiencyCalculator.getWeight(ParticleNo::TWO, track.pt());

if constexpr (std::is_same<PartType, FemtoRecoParticles>::value)
sameEventCont.setPair<true>(track, phicandidate, multCol, ConfUse3D, weight);
else
sameEventCont.setPair<false>(track, phicandidate, multCol, ConfUse3D, weight);
weight = effCorrection.getWeight<FilteredFDCollisions>(ParticleNo::ONE, phicandidate) * effCorrection.getWeight<FilteredFDCollisions>(ParticleNo::TWO, track);
sameEventCont.setPair<isMC>(track, phicandidate, multCol, ConfUse3D, weight);
}

// // Used for better fitting of invariant mass background.
Expand Down Expand Up @@ -592,7 +627,7 @@
// }
}

template <typename PartitionType, typename PartType, typename MCParticles = std::nullptr_t>
template <bool isMC, typename PartitionType, typename PartType, typename MCParticles = std::nullptr_t>
void doMixedEvent(PartitionType groupPartsTrack, PartitionType groupPartsPhi, PartType parts, float magFieldTesla, int multCol, [[maybe_unused]] MCParticles mcParts = nullptr)
{
for (auto const& [track, phicandidate] : combinations(CombinationsFullIndexPolicy(groupPartsTrack, groupPartsPhi))) {
Expand All @@ -612,31 +647,23 @@
continue;
}
}

weight = efficiencyCalculator.getWeight(ParticleNo::ONE, phicandidate.pt()) * efficiencyCalculator.getWeight(ParticleNo::TWO, track.pt());

if constexpr (std::is_same<PartType, FemtoRecoParticles>::value)
mixedEventCont.setPair<true>(track, phicandidate, multCol, ConfUse3D, weight);
else
mixedEventCont.setPair<false>(track, phicandidate, multCol, ConfUse3D, weight);
weight = effCorrection.getWeight<FilteredFDCollisions>(ParticleNo::ONE, phicandidate) * effCorrection.getWeight<FilteredFDCollisions>(ParticleNo::TWO, track);
mixedEventCont.setPair<isMC>(track, phicandidate, multCol, ConfUse3D, weight);
}
}

void processSameEvent(FilteredFDCollision const& col, FilteredFemtoFullParticles const& parts)
{
auto thegroupPartsTrack = partsTrack->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
auto thegroupPartsPhi = partsPhi->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
// auto thegroupPartsPhiDaugh = partsPhiDaugh->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
// auto thegroupPartsKaons = partsKaons->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
eventHisto.fillQA(col);
doSameEvent(thegroupPartsTrack, thegroupPartsPhi, parts, col.magField(), col.multNtr());
doSameEvent<false>(thegroupPartsTrack, thegroupPartsPhi, parts, col.magField(), col.multNtr());
}
PROCESS_SWITCH(FemtoUniversePairTaskTrackPhi, processSameEvent, "Enable processing same event", true);

void processMixedEvent(FilteredFDCollisions const& cols, FilteredFemtoFullParticles const& parts)
{
for (auto const& [collision1, collision2] : soa::selfCombinations(colBinning, 5, -1, cols, cols)) {

const int multiplicityCol = collision1.multNtr();
mixQaRegistry.fill(HIST("MixingQA/hMECollisionBins"), colBinning.getBin({collision1.posZ(), multiplicityCol}));

Expand All @@ -650,20 +677,19 @@
continue;
}

doMixedEvent(groupPartsTrack, groupPartsPhi, parts, magFieldTesla1, multiplicityCol);
doMixedEvent<false>(groupPartsTrack, groupPartsPhi, parts, magFieldTesla1, multiplicityCol);
}
}
PROCESS_SWITCH(FemtoUniversePairTaskTrackPhi, processMixedEvent, "Enable processing mixed events", true);

///--------------------------------------------MC-------------------------------------------------///
void processSameEventMCReco(FilteredFDCollision const& col, FemtoRecoParticles const& parts, aod::FdMCParticles const& mcparts)
{
eventHisto.fillQA(col);
// Reco
auto thegroupPartsTrack = partsTrackMCReco->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
auto thegroupPartsPhi = partsPhiMCReco->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
// auto thegroupPartsPhiDaugh = partsPhiDaugh->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
// auto thegroupPartsKaons = partsKaons->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
eventHisto.fillQA(col);
doSameEvent(thegroupPartsTrack, thegroupPartsPhi, parts, col.magField(), col.multNtr(), mcparts);
doSameEvent<true>(thegroupPartsTrack, thegroupPartsPhi, parts, col.magField(), col.multNtr(), mcparts);
}
PROCESS_SWITCH(FemtoUniversePairTaskTrackPhi, processSameEventMCReco, "Enable processing same event for MC Reco", true);

Expand All @@ -684,7 +710,7 @@
auto groupPartsTrack = partsTrackMCReco->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache);
auto groupPartsPhi = partsPhiMCReco->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache);

doMixedEvent(groupPartsTrack, groupPartsPhi, parts, magFieldTesla1, multiplicityCol, mcparts);
doMixedEvent<true>(groupPartsTrack, groupPartsPhi, parts, magFieldTesla1, multiplicityCol, mcparts);
}
}
PROCESS_SWITCH(FemtoUniversePairTaskTrackPhi, processMixedEventMCReco, "Enable processing mixed events for MC Reco", false);
Expand All @@ -705,18 +731,18 @@
// charge +
if (pdgParticle->Charge() > 0.0) {
registryMCtruth.fill(HIST("MCtruthAllPositivePt"), part.pt());
if (pdgCode == 2212) {

Check failure on line 734 in PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackPhi.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registryMCtruth.fill(HIST("MCtruthPpos"), part.pt(), part.eta());
registryMCtruth.fill(HIST("MCtruthPposPt"), part.pt());
continue;
} else if (pdgCode == 321) {

Check failure on line 738 in PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackPhi.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registryMCtruth.fill(HIST("MCtruthKp"), part.pt(), part.eta());
registryMCtruth.fill(HIST("MCtruthKpPt"), part.pt());
continue;
}
}
// charge 0
if (pdgCode == 333) {

Check failure on line 745 in PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackPhi.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registryMCtruth.fill(HIST("MCtruthPhi"), part.pt(), part.eta());
registryMCtruth.fill(HIST("MCtruthPhiPt"), part.pt());
continue;
Expand All @@ -726,11 +752,11 @@
if (pdgParticle->Charge() < 0.0) {
registryMCtruth.fill(HIST("MCtruthAllNegativePt"), part.pt());

if (pdgCode == -321) {

Check failure on line 755 in PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackPhi.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registryMCtruth.fill(HIST("MCtruthKm"), part.pt(), part.eta());
registryMCtruth.fill(HIST("MCtruthKmPt"), part.pt());
continue;
} else if (pdgCode == -2212) {

Check failure on line 759 in PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackPhi.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registryMCtruth.fill(HIST("MCtruthPneg"), part.pt(), part.eta());
registryMCtruth.fill(HIST("MCtruthPnegPt"), part.pt());
continue;
Expand All @@ -750,18 +776,18 @@

if (mcpart.pdgMCTruth() == ConfTrackPDGCode && (part.pt() > ConfTrackPtLow) && (part.pt() < ConfTrackPtHigh) && isParticleNSigmaAccepted(part.p(), trackCuts.getNsigmaTPC(part, o2::track::PID::Proton), trackCuts.getNsigmaTOF(part, o2::track::PID::Proton), trackCuts.getNsigmaTPC(part, o2::track::PID::Pion), trackCuts.getNsigmaTOF(part, o2::track::PID::Pion), trackCuts.getNsigmaTPC(part, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(part, o2::track::PID::Kaon))) {
registryMCpT.fill(HIST("MCReco/NC_p_pT"), part.pt());
float weightTrack = efficiencyCalculator.getWeight(ParticleNo::TWO, part.pt());
float weightTrack = effCorrection.getWeight<FilteredFDCollisions>(ParticleNo::TWO, part);
registryMCpT.fill(HIST("MCReco/C_p_pT"), part.pt(), weightTrack);
}
if ((mcpart.pdgMCTruth() == 333) && (part.partType() == aod::femtouniverseparticle::ParticleType::kPhi) && (part.pt() > ConfPhiPtLow) && (part.pt() < ConfPhiPtHigh)) {

Check failure on line 782 in PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackPhi.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registryMCpT.fill(HIST("MCReco/NC_phi_pT"), part.pt());
float weightPhi = efficiencyCalculator.getWeight(ParticleNo::ONE, part.pt());
float weightPhi = effCorrection.getWeight<FilteredFDCollisions>(ParticleNo::ONE, part);
registryMCpT.fill(HIST("MCReco/C_phi_pT"), part.pt(), weightPhi);
}

if (isParticleNSigmaAccepted(part.p(), trackCuts.getNsigmaTPC(part, o2::track::PID::Proton), trackCuts.getNsigmaTOF(part, o2::track::PID::Proton), trackCuts.getNsigmaTPC(part, o2::track::PID::Pion), trackCuts.getNsigmaTOF(part, o2::track::PID::Pion), trackCuts.getNsigmaTPC(part, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(part, o2::track::PID::Kaon)))
hTrackDCA.fillQA<true, true>(part);
if ((part.partType() == aod::femtouniverseparticle::ParticleType::kPhi) && (mcpart.pdgMCTruth() == 333) && (mcpart.partOriginMCTruth() == aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kPrimary)) {

Check failure on line 790 in PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackPhi.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registryMCreco.fill(HIST("MCrecoPhi"), mcpart.pt(), mcpart.eta()); // phi
registryMCreco.fill(HIST("MCrecoPhiPt"), mcpart.pt());
} else if (part.partType() == aod::femtouniverseparticle::ParticleType::kTrack) {
Expand All @@ -778,7 +804,6 @@
registryMCreco.fill(HIST("MCrecoPnegPt"), mcpart.pt());
}
}

} // partType kTrack
}
}
Expand Down
Loading