diff --git a/PWGCF/FemtoDream/Tasks/femtoDreamPairEfficiency.cxx b/PWGCF/FemtoDream/Tasks/femtoDreamPairEfficiency.cxx index 44fb6b8cba0..620c7b33839 100644 --- a/PWGCF/FemtoDream/Tasks/femtoDreamPairEfficiency.cxx +++ b/PWGCF/FemtoDream/Tasks/femtoDreamPairEfficiency.cxx @@ -83,6 +83,12 @@ struct FemtoDreamPairEfficiency { SliceCache cache; Preslice perColMc = mcparticle::mcCollisionId; + struct : ConfigurableGroup { + std::string prefix = std::string("Mode"); + Configurable countPairs{"countPairs", true, "Count Pairs"}; + Configurable countTriplets{"countTriplets", false, "Count Triplets"}; + } mode; + // Event Selections struct : ConfigurableGroup { std::string prefix = std::string("EventSelection"); @@ -90,6 +96,7 @@ struct FemtoDreamPairEfficiency { Configurable offlineCheck{"offlineCheck", true, "Check for Sel8"}; Configurable etaAbsMax{"etaAbsMax", 0.8f, "Common eta cut for particles in dNdEta distribution"}; Configurable kstarMax{"kstarMax", 999.f, "Cut on kstar"}; + Configurable q3Max{"q3Max", 999.f, "Cut on Q3"}; } eventSelection; struct : ConfigurableGroup { @@ -146,6 +153,27 @@ struct FemtoDreamPairEfficiency { Configurable tpctofNsigmaMax{"tpctofNsigmaMax", 3.f, "TPCTOC nsigma max"}; } trackCuts2; + struct : ConfigurableGroup { + std::string prefix = std::string("SelectionTrack3"); + Configurable sign{"sign", 1, "Sign of charge"}; + Configurable ptMin{"ptMin", 0.0f, "pt min"}; + Configurable ptMax{"ptMax", 3.0f, "pt max"}; + Configurable etaAbsMax{"etaAbsMax", 0.8f, "|eta| max"}; + Configurable dcazAbsMax{"dcazAbsMax", 0.1f, "|dca_z| max"}; + Configurable useDcaxyPtDepCut{"useDcaxyPtDepCut", true, "|dca_z| max"}; + Configurable tpcClusterMin{"tpcClusterMin", 80.f, "TPC clusters min"}; + Configurable tpcCrossedOverClusterMin{"tpcCrossedOverClusterMin", 0.83f, "TPC clusters/TPC crossed rows min"}; + Configurable tpcCrossedMin{"tpcCrossedMin", 70.f, "TPC crossed rows min"}; + Configurable tpcSharedMax{"tpcSharedMax", 160.f, "TPC shared clusters max"}; + Configurable itsClusterMin{"itsClusterMin", 0.f, "ITS clusters min"}; + Configurable itsIbClusterMin{"itsIbClusterMin", 0.f, "ITS inner barrle min"}; + Configurable pdgCode{"pdgCode", 2212, "PDG code"}; + Configurable pidThreshold{"pidThreshold", 0.75f, "Momentum threshold for PID"}; + Configurable itsNsigmaMax{"itsNsigmaMax", 99.f, "its nsigma max"}; + Configurable tpcNsigmaMax{"tpcNsigmaMax", 3.f, "TPC nsigma max"}; + Configurable tpctofNsigmaMax{"tpctofNsigmaMax", 3.f, "TPCTOC nsigma max"}; + } trackCuts3; + HistogramRegistry dataEventHist{"dataEventHist", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; HistogramRegistry mcEventHist{"mcEventHist", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; @@ -153,6 +181,11 @@ struct FemtoDreamPairEfficiency { 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(HIST("hDataEventSelection"))->GetXaxis()->SetBinLabel(1, "All collisions"); dataEventHist.get(HIST("hDataEventSelection"))->GetXaxis()->SetBinLabel(2, "sel8 cut"); @@ -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(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(1, "All collisions"); mcEventHist.get(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(2, "Sel8 cut"); mcEventHist.get(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(3, "posZ cut"); mcEventHist.get(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(4, "INEL>0 cut"); mcEventHist.get(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(5, "With at least a gen coll"); - mcEventHist.get(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(6, "With at least a pair"); - mcEventHist.get(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(7, "With at least a lowkstar pair"); + mcEventHist.get(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(6, "With at least a pair/triplet"); + mcEventHist.get(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(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(1, "All collisions"); mcEventHist.get(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(2, "posZ cut"); mcEventHist.get(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(3, "INEL>0 cut"); - mcEventHist.get(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(4, "With at least a pair"); - mcEventHist.get(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(5, "With at least a lowkstar pair"); + mcEventHist.get(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(4, "With at least a pair/triplet"); + mcEventHist.get(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(5, "With at least a lowkstar pair/lowq3 triplet"); mcEventHist.get(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(6, "With at least a reco coll"); // MC Event information for Rec and Gen @@ -236,6 +283,9 @@ struct FemtoDreamPairEfficiency { int pairsFound = 0; std::vector vecKstar{}; + int tripletsFound = 0; + std::vector vecq3{}; + template bool checkRecoTrackSelections(const T1& track, const T2& sel) { @@ -417,6 +467,72 @@ struct FemtoDreamPairEfficiency { return pairs; } + template + 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 int countGenPairs(const T& tracks) { @@ -444,6 +560,42 @@ struct FemtoDreamPairEfficiency { return pairs; } + template + 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(collision, true)) @@ -452,16 +604,33 @@ struct FemtoDreamPairEfficiency { auto tracksWithItsPid = soa::Attach(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 @@ -491,17 +660,35 @@ struct FemtoDreamPairEfficiency { const auto& mcCollision = collision.mcCollision_as(); 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(); @@ -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());