Skip to content

[PWGCF] Update pair efficiency task #13728

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
ariedel-cern merged 1 commit into from
Nov 8, 2025
Merged

[PWGCF] Update pair efficiency task #13728

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
257 changes: 232 additions & 25 deletions PWGCF/FemtoDream/Tasks/femtoDreamPairEfficiency.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,13 +83,20 @@ struct FemtoDreamPairEfficiency {
SliceCache cache;
Preslice<McParticles> perColMc = mcparticle::mcCollisionId;

struct : ConfigurableGroup {
std::string prefix = std::string("Mode");
Configurable<bool> countPairs{"countPairs", true, "Count Pairs"};
Configurable<bool> countTriplets{"countTriplets", false, "Count Triplets"};
} mode;

// Event Selections
struct : ConfigurableGroup {
std::string prefix = std::string("EventSelection");
Configurable<float> zvtxAbsMax{"zvtxAbsMax", 10.f, "|z-Vertex| max"};
Configurable<bool> offlineCheck{"offlineCheck", true, "Check for Sel8"};
Configurable<float> etaAbsMax{"etaAbsMax", 0.8f, "Common eta cut for particles in dNdEta distribution"};
Configurable<float> kstarMax{"kstarMax", 999.f, "Cut on kstar"};
Configurable<float> q3Max{"q3Max", 999.f, "Cut on Q3"};
} eventSelection;

struct : ConfigurableGroup {
Expand Down Expand Up @@ -146,13 +153,39 @@ struct FemtoDreamPairEfficiency {
Configurable<float> tpctofNsigmaMax{"tpctofNsigmaMax", 3.f, "TPCTOC nsigma max"};
} trackCuts2;

struct : ConfigurableGroup {
std::string prefix = std::string("SelectionTrack3");
Configurable<int> sign{"sign", 1, "Sign of charge"};
Configurable<float> ptMin{"ptMin", 0.0f, "pt min"};
Configurable<float> ptMax{"ptMax", 3.0f, "pt max"};
Configurable<float> etaAbsMax{"etaAbsMax", 0.8f, "|eta| max"};
Configurable<float> dcazAbsMax{"dcazAbsMax", 0.1f, "|dca_z| max"};
Configurable<bool> useDcaxyPtDepCut{"useDcaxyPtDepCut", true, "|dca_z| max"};
Configurable<float> tpcClusterMin{"tpcClusterMin", 80.f, "TPC clusters min"};
Configurable<float> tpcCrossedOverClusterMin{"tpcCrossedOverClusterMin", 0.83f, "TPC clusters/TPC crossed rows min"};
Configurable<float> tpcCrossedMin{"tpcCrossedMin", 70.f, "TPC crossed rows min"};
Configurable<float> tpcSharedMax{"tpcSharedMax", 160.f, "TPC shared clusters max"};
Configurable<float> itsClusterMin{"itsClusterMin", 0.f, "ITS clusters min"};
Configurable<float> itsIbClusterMin{"itsIbClusterMin", 0.f, "ITS inner barrle min"};
Configurable<int> pdgCode{"pdgCode", 2212, "PDG code"};
Configurable<float> pidThreshold{"pidThreshold", 0.75f, "Momentum threshold for PID"};
Configurable<float> itsNsigmaMax{"itsNsigmaMax", 99.f, "its nsigma max"};
Configurable<float> tpcNsigmaMax{"tpcNsigmaMax", 3.f, "TPC nsigma max"};
Configurable<float> tpctofNsigmaMax{"tpctofNsigmaMax", 3.f, "TPCTOC nsigma max"};
} trackCuts3;

HistogramRegistry dataEventHist{"dataEventHist", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
HistogramRegistry mcEventHist{"mcEventHist", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};

HistogramRegistry dataTrackHist{"dataTrackHist", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};

void init(InitContext&)
{

if (mode.countPairs.value && mode.countTriplets.value) {
LOG(fatal) << "We cannot count pairs and triplets at the same time. Turn one of the off.";
}

dataEventHist.add("hDataEventSelection", "hDataEventSelection", kTH1F, {{6, -0.5f, 5.5f}});
dataEventHist.get<TH1>(HIST("hDataEventSelection"))->GetXaxis()->SetBinLabel(1, "All collisions");
dataEventHist.get<TH1>(HIST("hDataEventSelection"))->GetXaxis()->SetBinLabel(2, "sel8 cut");
Expand Down Expand Up @@ -193,21 +226,35 @@ struct FemtoDreamPairEfficiency {
dataTrackHist.add("Track2/tpcNsigma", "Track 2 nsigma tpc", kTH2F, {{600, 0, 6}, {1000, -5, 5}});
dataTrackHist.add("Track2/tpctofNsigma", "Track 2 nsigma tpctof", kTH2F, {{600, 0, 6}, {500, 0, 5}});

if (mode.countTriplets) {
dataTrackHist.add("Track3/pt", "Track 3 pt", kTH1F, {{600, 0, 6}});
dataTrackHist.add("Track3/eta", "Track 3 eta", kTH1F, {{200, -1, 1}});
dataTrackHist.add("Track3/phi", "Track 3 phi", kTH1F, {{720, 0, o2::constants::math::TwoPI}});
dataTrackHist.add("Track3/tpcCluster", "Track 3 cluster", kTH1F, {{160, 0, 160}});
dataTrackHist.add("Track3/tpcCrossed", "Track 3 crossed rows", kTH1F, {{160, 0, 160}});
dataTrackHist.add("Track3/tpcShared", "Track 3 cluster shared", kTH1F, {{160, 0, 160}});
dataTrackHist.add("Track3/itsCluster", "Track 3 its cluster", kTH1F, {{0, 0, 7}});
dataTrackHist.add("Track3/itsIbCluster", "Track 3 its cluster inner barrel", kTH1F, {{3, 0, 3}});
dataTrackHist.add("Track3/itsNsigma", "Track 3 nsigma its", kTH2F, {{600, 0, 6}, {1000, -5, 5}});
dataTrackHist.add("Track3/tpcNsigma", "Track 3 nsigma tpc", kTH2F, {{600, 0, 6}, {1000, -5, 5}});
dataTrackHist.add("Track3/tpctofNsigma", "Track 3 nsigma tpctof", kTH2F, {{600, 0, 6}, {500, 0, 5}});
}

mcEventHist.add("hRecoEventSelection", "hRecoEventSelection", kTH1F, {{7, -0.5f, 6.5f}});
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(1, "All collisions");
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(2, "Sel8 cut");
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(3, "posZ cut");
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(4, "INEL>0 cut");
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(5, "With at least a gen coll");
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(6, "With at least a pair");
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(7, "With at least a lowkstar pair");
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(6, "With at least a pair/triplet");
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(7, "With at least a lowkstar pair/lowq3 triplet");

mcEventHist.add("hGenEventSelection", "hGenEventSelection", kTH1F, {{6, -0.5f, 5.5f}});
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(1, "All collisions");
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(2, "posZ cut");
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(3, "INEL>0 cut");
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(4, "With at least a pair");
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(5, "With at least a lowkstar pair");
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(4, "With at least a pair/triplet");
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(5, "With at least a lowkstar pair/lowq3 triplet");
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(6, "With at least a reco coll");

// MC Event information for Rec and Gen
Expand Down Expand Up @@ -236,6 +283,9 @@ struct FemtoDreamPairEfficiency {
int pairsFound = 0;
std::vector<float> vecKstar{};

int tripletsFound = 0;
std::vector<float> vecq3{};

template <typename T1, typename T2>
bool checkRecoTrackSelections(const T1& track, const T2& sel)
{
Expand Down Expand Up @@ -417,6 +467,72 @@ struct FemtoDreamPairEfficiency {
return pairs;
}

template <typename T>
int countRecoTriplets(const T& tracks)
{
int triplets = 0;
for (auto track1 = tracks.begin(); track1 != tracks.end(); track1++) {
if (!checkRecoTrackSelections(track1, trackCuts1) || !checkRecoTrackPidSelections(track1, trackCuts1)) {
continue;
}
for (auto track2 = track1 + 1; track2 != tracks.end(); track2++) {
if (!checkRecoTrackSelections(track2, trackCuts2) || !checkRecoTrackPidSelections(track2, trackCuts2)) {
continue;
}

for (auto track3 = track2 + 1; track3 != tracks.end(); track3++) {
if (!checkRecoTrackSelections(track3, trackCuts3) || !checkRecoTrackPidSelections(track3, trackCuts3)) {
continue;
}

auto nsigma1 = getNSigmaValues(track1, trackCuts1.pdgCode.value);
dataTrackHist.fill(HIST("Track1/pt"), track1.pt());
dataTrackHist.fill(HIST("Track1/eta"), track1.eta());
dataTrackHist.fill(HIST("Track1/phi"), track1.phi());
dataTrackHist.fill(HIST("Track1/tpcCluster"), track1.tpcNClsFound());
dataTrackHist.fill(HIST("Track1/tpcCrossed"), track1.tpcNClsCrossedRows());
dataTrackHist.fill(HIST("Track1/tpcShared"), track1.tpcNClsShared());
dataTrackHist.fill(HIST("Track1/itsCluster"), track1.itsNCls());
dataTrackHist.fill(HIST("Track1/itsIbCluster"), track1.itsNClsInnerBarrel());
dataTrackHist.fill(HIST("Track1/itsNsigma"), track1.p(), nsigma1[0]);
dataTrackHist.fill(HIST("Track1/tpcNsigma"), track1.p(), nsigma1[1]);
dataTrackHist.fill(HIST("Track1/tpctofNsigma"), track1.p(), std::hypot(nsigma1[1], nsigma1[2]));

auto nsigma2 = getNSigmaValues(track2, trackCuts2.pdgCode.value);
dataTrackHist.fill(HIST("Track2/pt"), track2.pt());
dataTrackHist.fill(HIST("Track2/eta"), track2.eta());
dataTrackHist.fill(HIST("Track2/phi"), track2.phi());
dataTrackHist.fill(HIST("Track2/tpcCluster"), track2.tpcNClsFound());
dataTrackHist.fill(HIST("Track2/tpcCrossed"), track2.tpcNClsCrossedRows());
dataTrackHist.fill(HIST("Track2/tpcShared"), track2.tpcNClsShared());
dataTrackHist.fill(HIST("Track2/itsCluster"), track2.itsNCls());
dataTrackHist.fill(HIST("Track2/itsIbCluster"), track2.itsNClsInnerBarrel());
dataTrackHist.fill(HIST("Track2/itsNsigma"), track2.p(), nsigma2[0]);
dataTrackHist.fill(HIST("Track2/tpcNsigma"), track2.p(), nsigma2[1]);
dataTrackHist.fill(HIST("Track2/tpctofNsigma"), track2.p(), std::hypot(nsigma2[1], nsigma2[2]));

auto nsigma3 = getNSigmaValues(track3, trackCuts3.pdgCode.value);
dataTrackHist.fill(HIST("Track3/pt"), track3.pt());
dataTrackHist.fill(HIST("Track3/eta"), track3.eta());
dataTrackHist.fill(HIST("Track3/phi"), track3.phi());
dataTrackHist.fill(HIST("Track3/tpcCluster"), track3.tpcNClsFound());
dataTrackHist.fill(HIST("Track3/tpcCrossed"), track3.tpcNClsCrossedRows());
dataTrackHist.fill(HIST("Track3/tpcShared"), track3.tpcNClsShared());
dataTrackHist.fill(HIST("Track3/itsCluster"), track3.itsNCls());
dataTrackHist.fill(HIST("Track3/itsIbCluster"), track3.itsNClsInnerBarrel());
dataTrackHist.fill(HIST("Track3/itsNsigma"), track3.p(), nsigma2[0]);
dataTrackHist.fill(HIST("Track3/tpcNsigma"), track3.p(), nsigma2[1]);
dataTrackHist.fill(HIST("Track3/tpctofNsigma"), track3.p(), std::hypot(nsigma3[1], nsigma3[2]));

triplets++;

vecq3.push_back(femtoDream::FemtoDreamMath::getQ3(track1, femtoDream::getMass(trackCuts1.pdgCode.value), track2, femtoDream::getMass(trackCuts2.pdgCode.value), track3, femtoDream::getMass(trackCuts3.pdgCode.value)));
}
}
}
return triplets;
}

template <typename T>
int countGenPairs(const T& tracks)
{
Expand Down Expand Up @@ -444,6 +560,42 @@ struct FemtoDreamPairEfficiency {
return pairs;
}

template <typename T>
int countGenTriplets(const T& tracks)
{
int triplets = 0;
for (auto track1 = tracks.begin(); track1 != tracks.end(); track1++) {
if (!track1.isPhysicalPrimary() ||
track1.pdgCode() != trackCuts1.pdgCode.value ||
track1.pt() < trackCuts1.ptMin.value ||
track1.pt() > trackCuts1.ptMax.value ||
std::fabs(track1.eta()) > trackCuts1.etaAbsMax.value) {
continue;
}
for (auto track2 = track1 + 1; track2 != tracks.end(); track2++) {
if (!track2.isPhysicalPrimary() ||
track2.pdgCode() != trackCuts2.pdgCode.value ||
track2.pt() < trackCuts2.ptMin.value ||
track2.pt() > trackCuts2.ptMax.value ||
std::fabs(track2.eta()) > trackCuts2.etaAbsMax.value) {
continue;
}
for (auto track3 = track2 + 1; track3 != tracks.end(); track3++) {
if (!track3.isPhysicalPrimary() ||
track3.pdgCode() != trackCuts3.pdgCode.value ||
track3.pt() < trackCuts3.ptMin.value ||
track3.pt() > trackCuts3.ptMax.value ||
std::fabs(track3.eta()) > trackCuts3.etaAbsMax.value) {
continue;
}
vecq3.push_back(femtoDream::FemtoDreamMath::getQ3(track1, femtoDream::getMass(trackCuts1.pdgCode.value), track2, femtoDream::getMass(trackCuts2.pdgCode.value), track3, femtoDream::getMass(trackCuts3.pdgCode.value)));
triplets++;
}
}
}
return triplets;
}

void processdNdetaData(Collisions::iterator const& collision, FilteredFullTracks const& tracks)
{
if (!checkEventSelections<false>(collision, true))
Expand All @@ -452,16 +604,33 @@ struct FemtoDreamPairEfficiency {
auto tracksWithItsPid = soa::Attach<FilteredFullTracks, aod::pidits::ITSNSigmaEl, aod::pidits::ITSNSigmaPi, aod::pidits::ITSNSigmaKa,
aod::pidits::ITSNSigmaPr, aod::pidits::ITSNSigmaDe, aod::pidits::ITSNSigmaTr, aod::pidits::ITSNSigmaHe>(tracks);

vecKstar.clear();
pairsFound = countRecoPairs(tracksWithItsPid);

if (pairsFound == 0) {
return;
if (mode.countPairs.value) {
vecKstar.clear();
pairsFound = countRecoPairs(tracksWithItsPid);
if (pairsFound == 0) {
return;
}
}
if (mode.countTriplets.value) {
vecq3.clear();
tripletsFound = countRecoTriplets(tracksWithItsPid);
if (tripletsFound == 0) {
return;
}
}

dataEventHist.fill(HIST("hDataEventSelection"), 4); // found a pair

if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
return;
if (mode.countPairs.value) {
if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
return;
}
}

if (mode.countTriplets.value) {
if (*std::min_element(vecq3.begin(), vecq3.end()) > eventSelection.q3Max) {
return;
}
}

dataEventHist.fill(HIST("hDataEventSelection"), 5); // found a lowkstar pair
Expand Down Expand Up @@ -491,17 +660,35 @@ struct FemtoDreamPairEfficiency {
const auto& mcCollision = collision.mcCollision_as<GenCollisions>();
auto mcParticlesThisColl = mcParticles.sliceByCached(mcparticle::mcCollisionId, mcCollision.globalIndex(), cache);

vecKstar.clear();
pairsFound = countGenPairs(mcParticlesThisColl);

if (pairsFound == 0) {
return;
if (mode.countPairs.value) {
vecKstar.clear();
pairsFound = countGenPairs(mcParticlesThisColl);
if (pairsFound == 0) {
return;
}
}
if (mode.countTriplets.value) {
vecq3.clear();
tripletsFound = countGenTriplets(mcParticlesThisColl);
if (tripletsFound == 0) {
return;
}
}

mcEventHist.fill(HIST("hRecoEventSelection"), 5);

if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
return;
if (mode.countPairs.value) {
if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
return;
}
}

if (mode.countTriplets.value) {
if (*std::min_element(vecq3.begin(), vecq3.end()) > eventSelection.q3Max) {
return;
}
}

mcEventHist.fill(HIST("hRecoEventSelection"), 6);

float genCentrality = mcCollision.centFT0M();
Expand Down Expand Up @@ -551,17 +738,37 @@ struct FemtoDreamPairEfficiency {
return;
}
mcEventHist.fill(HIST("hGenEventSelection"), 2);
vecKstar.clear();
pairsFound = countGenPairs(mcParticles);
if (pairsFound == 0) {
return;

if (mode.countPairs.value) {
vecKstar.clear();
pairsFound = countGenPairs(mcParticles);
if (pairsFound == 0) {
return;
}
}
if (mode.countTriplets.value) {
vecq3.clear();
tripletsFound = countGenTriplets(mcParticles);
if (tripletsFound == 0) {
return;
}
}

mcEventHist.fill(HIST("hGenEventSelection"), 3);
if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
return;

if (mode.countPairs.value) {
if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
return;
}
}

if (mode.countTriplets.value) {
if (*std::min_element(vecq3.begin(), vecq3.end()) > eventSelection.q3Max) {
return;
}
}
mcEventHist.fill(HIST("hGenEventSelection"), 4);

mcEventHist.fill(HIST("hGenEventSelection"), 4);
mcEventHist.fill(HIST("hGenMcVertexZ"), mcCollision.posZ());
mcEventHist.fill(HIST("hGenMcMultiplicity"), mcCollision.multMCNParticlesEta08());
mcEventHist.fill(HIST("hGenMcCentrality"), mcCollision.centFT0M());
Expand Down
Loading