Skip to content

Commit c6f333e

Browse files
skundu692skundu692
authored
[PWGLF] Add new process function for doublephi and modify charge pairing for f1p correlation (#12957)
1 parent e750600 commit c6f333e

File tree

2 files changed

+229
-19
lines changed

2 files changed

+229
-19
lines changed

‎PWGLF/Tasks/Resonances/doublephimeson.cxx‎

Lines changed: 217 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -13,27 +13,32 @@
1313
/// \author sourav kundu
1414
/// \since 02/11/2023
1515

16+
#include "PWGLF/DataModel/ReducedDoublePhiTables.h"
17+
18+
#include "Common/Core/trackUtilities.h"
19+
20+
#include "CommonConstants/PhysicsConstants.h"
21+
#include "Framework/ASoAHelpers.h"
22+
#include "Framework/AnalysisDataModel.h"
23+
#include "Framework/AnalysisTask.h"
24+
#include "Framework/StepTHn.h"
25+
#include "Framework/runDataProcessing.h"
1626
#include <Framework/Configurable.h>
17-
#include <TLorentzVector.h>
27+
1828
#include <Math/GenVector/Boost.h>
19-
#include <Math/Vector4D.h>
2029
#include <Math/Vector3D.h>
30+
#include <Math/Vector4D.h>
31+
#include <TLorentzVector.h>
2132
#include <TMath.h>
33+
#include <TVector2.h>
34+
2235
#include <fairlogger/Logger.h>
36+
2337
#include <iostream>
2438
#include <iterator>
2539
#include <string>
2640
#include <vector>
2741

28-
#include "Framework/AnalysisTask.h"
29-
#include "Framework/ASoAHelpers.h"
30-
#include "Framework/runDataProcessing.h"
31-
#include "Framework/AnalysisDataModel.h"
32-
#include "Framework/StepTHn.h"
33-
#include "Common/Core/trackUtilities.h"
34-
#include "PWGLF/DataModel/ReducedDoublePhiTables.h"
35-
#include "CommonConstants/PhysicsConstants.h"
36-
3742
using namespace o2;
3843
using namespace o2::framework;
3944
using namespace o2::framework::expressions;
@@ -631,6 +636,207 @@ struct doublephimeson {
631636
}
632637
PROCESS_SWITCH(doublephimeson, processopti, "Process Optimized same event", false);
633638

639+
void processopti2(aod::RedPhiEvents::iterator const& collision, aod::PhiTracks const& phitracks)
640+
{
641+
if (additionalEvsel && (collision.numPos() < 2 || collision.numNeg() < 2)) {
642+
return;
643+
}
644+
645+
// === Helpers ===
646+
// PDG phi mass (centralized constant, easy to vary for systematics)
647+
constexpr double mPhiPDG = 1.019461; // GeV/c^2
648+
649+
// ΔR with proper Δφ wrapping (−π, π]
650+
const auto deltaR = [](double phi1, double eta1, double phi2, double eta2) {
651+
const double dphi = TVector2::Phi_mpi_pi(phi1 - phi2);
652+
const double deta = eta1 - eta2;
653+
return std::sqrt(dphi * dphi + deta * deta);
654+
};
655+
656+
// Combined mass offset from PDG (tie-breaker metric)
657+
const auto deltam = [=](double m1, double m2) {
658+
return std::sqrt((m1 - mPhiPDG) * (m1 - mPhiPDG) + (m2 - mPhiPDG) * (m2 - mPhiPDG));
659+
};
660+
661+
// Anti-merging/duplication veto:
662+
// Require same-sign kaons from the two φ's to be separated in (η,φ).
663+
// This protects against split/merged tracks and near-duplicate topologies.
664+
const auto passDaughterDR = [&](const ROOT::Math::PtEtaPhiMVector& kplusA,
665+
const ROOT::Math::PtEtaPhiMVector& kplusB,
666+
const ROOT::Math::PtEtaPhiMVector& kminusA,
667+
const ROOT::Math::PtEtaPhiMVector& kminusB) {
668+
const double dRkplus = deltaR(kplusA.Phi(), kplusA.Eta(), kplusB.Phi(), kplusB.Eta());
669+
const double dRkminus = deltaR(kminusA.Phi(), kminusA.Eta(), kminusB.Phi(), kminusB.Eta());
670+
return (dRkplus > daughterDeltaR) && (dRkminus > daughterDeltaR);
671+
// If later needed, make pT-aware:
672+
// const double thr = std::max(0.01, daughterDeltaR - 0.002 * std::min(10.0, exoticPt));
673+
// return (dRkplus > thr) && (dRkminus > thr);
674+
};
675+
676+
// === Phi multiplicity (uses PID1 for d1, PID2 for d2) ===
677+
int phimult = 0;
678+
for (auto const& t : phitracks) {
679+
if (t.phiMass() < minPhiMass1 || t.phiMass() > maxPhiMass1)
680+
continue;
681+
const double kpluspt = std::hypot(t.phid1Px(), t.phid1Py());
682+
const double kminuspt = std::hypot(t.phid2Px(), t.phid2Py());
683+
if (kpluspt > maxKaonPt)
684+
continue;
685+
if (kminuspt > maxKaonPt)
686+
continue;
687+
if (!selectionPID(t.phid1TPC(), t.phid1TOF(), t.phid1TOFHit(), strategyPID1, kpluspt))
688+
continue; // d1 → PID1
689+
if (!selectionPID(t.phid2TPC(), t.phid2TOF(), t.phid2TOFHit(), strategyPID2, kminuspt))
690+
continue; // d2 → PID2
691+
++phimult;
692+
}
693+
694+
// === Collect candidates first ===
695+
std::vector<ROOT::Math::PtEtaPhiMVector> exoticres, phi1v, phi2v, kplus1, kplus2, kminus1, kminus2;
696+
std::vector<int> d1id, d2id, d3id, d4id;
697+
698+
const auto n = phitracks.size();
699+
exoticres.reserve(n);
700+
phi1v.reserve(n);
701+
phi2v.reserve(n);
702+
kplus1.reserve(n);
703+
kplus2.reserve(n);
704+
kminus1.reserve(n);
705+
kminus2.reserve(n);
706+
d1id.reserve(n);
707+
d2id.reserve(n);
708+
d3id.reserve(n);
709+
d4id.reserve(n);
710+
711+
for (auto const& t1 : phitracks) {
712+
// Per-φ selection for the first φ
713+
const double kplus1pt = std::hypot(t1.phid1Px(), t1.phid1Py());
714+
const double kminus1pt = std::hypot(t1.phid2Px(), t1.phid2Py());
715+
if (kplus1pt > maxKaonPt)
716+
continue;
717+
if (kminus1pt > maxKaonPt)
718+
continue;
719+
if (!selectionPID(t1.phid1TPC(), t1.phid1TOF(), t1.phid1TOFHit(), strategyPID1, kplus1pt))
720+
continue;
721+
if (!selectionPID(t1.phid2TPC(), t1.phid2TOF(), t1.phid2TOFHit(), strategyPID2, kminus1pt))
722+
continue;
723+
724+
// Set vectors BEFORE QA fills
725+
Phid1.SetXYZM(t1.phiPx(), t1.phiPy(), t1.phiPz(), t1.phiMass());
726+
Phi1kaonplus.SetXYZM(t1.phid1Px(), t1.phid1Py(), t1.phid1Pz(), 0.493);
727+
Phi1kaonminus.SetXYZM(t1.phid2Px(), t1.phid2Py(), t1.phid2Pz(), 0.493);
728+
729+
// PID QA
730+
histos.fill(HIST("hnsigmaTPCTOFKaon"), t1.phid1TPC(), t1.phid1TOF(), kplus1pt);
731+
histos.fill(HIST("hnsigmaTPCKaonPlus"), t1.phid1TPC(), kplus1pt);
732+
histos.fill(HIST("hnsigmaTPCKaonMinus"), t1.phid2TPC(), kminus1pt);
733+
histos.fill(HIST("hPhiMass"), Phid1.M(), Phid1.Pt());
734+
735+
const auto id1 = t1.index();
736+
737+
for (auto const& t2 : phitracks) {
738+
const auto id2 = t2.index();
739+
if (id2 <= id1)
740+
continue; // unique unordered pairs
741+
742+
// Per-φ selection for the second φ
743+
const double kplus2pt = std::hypot(t2.phid1Px(), t2.phid1Py());
744+
const double kminus2pt = std::hypot(t2.phid2Px(), t2.phid2Py());
745+
if (kplus2pt > maxKaonPt)
746+
continue;
747+
if (kminus2pt > maxKaonPt)
748+
continue;
749+
if (!selectionPID(t2.phid1TPC(), t2.phid1TOF(), t2.phid1TOFHit(), strategyPID1, kplus2pt))
750+
continue;
751+
if (!selectionPID(t2.phid2TPC(), t2.phid2TOF(), t2.phid2TOFHit(), strategyPID2, kminus2pt))
752+
continue;
753+
754+
// Disallow shared same-sign daughters between the two φ's
755+
if ((t1.phid1Index() == t2.phid1Index()) || (t1.phid2Index() == t2.phid2Index()))
756+
continue;
757+
758+
Phid2.SetXYZM(t2.phiPx(), t2.phiPy(), t2.phiPz(), t2.phiMass());
759+
Phi2kaonplus.SetXYZM(t2.phid1Px(), t2.phid1Py(), t2.phid1Pz(), 0.493);
760+
Phi2kaonminus.SetXYZM(t2.phid2Px(), t2.phid2Py(), t2.phid2Pz(), 0.493);
761+
762+
// Mass windows
763+
if (t1.phiMass() < minPhiMass1 || t1.phiMass() > maxPhiMass1)
764+
continue;
765+
if (t2.phiMass() < minPhiMass2 || t2.phiMass() > maxPhiMass2)
766+
continue;
767+
768+
exotic = Phid1 + Phid2;
769+
if (exotic.M() < minExoticMass || exotic.M() > maxExoticMass)
770+
continue;
771+
772+
// Store candidate and bookkeeping
773+
exoticres.emplace_back(exotic.Pt(), exotic.Eta(), exotic.Phi(), exotic.M());
774+
phi1v.emplace_back(Phid1.Pt(), Phid1.Eta(), Phid1.Phi(), Phid1.M());
775+
phi2v.emplace_back(Phid2.Pt(), Phid2.Eta(), Phid2.Phi(), Phid2.M());
776+
777+
kplus1.emplace_back(Phi1kaonplus.Pt(), Phi1kaonplus.Eta(), Phi1kaonplus.Phi(), 0.493);
778+
kminus1.emplace_back(Phi1kaonminus.Pt(), Phi1kaonminus.Eta(), Phi1kaonminus.Phi(), 0.493);
779+
kplus2.emplace_back(Phi2kaonplus.Pt(), Phi2kaonplus.Eta(), Phi2kaonplus.Phi(), 0.493);
780+
kminus2.emplace_back(Phi2kaonminus.Pt(), Phi2kaonminus.Eta(), Phi2kaonminus.Phi(), 0.493);
781+
782+
d1id.push_back(t1.phid1Index());
783+
d2id.push_back(t2.phid1Index());
784+
d3id.push_back(t1.phid2Index());
785+
d4id.push_back(t2.phid2Index());
786+
}
787+
}
788+
789+
if (exoticres.empty())
790+
return;
791+
792+
// === Special handling if exactly two candidates were built ===
793+
if (exoticres.size() == 2) {
794+
const int i0 = 0, i1 = 1;
795+
796+
const bool sameFour =
797+
((d1id[i0] == d1id[i1] || d1id[i0] == d2id[i1]) &&
798+
(d2id[i0] == d1id[i1] || d2id[i0] == d2id[i1]) &&
799+
(d3id[i0] == d3id[i1] || d3id[i0] == d4id[i1]) &&
800+
(d4id[i0] == d3id[i1] || d4id[i0] == d4id[i1]));
801+
802+
const double deltam0 = deltam(phi1v[i0].M(), phi2v[i0].M());
803+
const double deltam1 = deltam(phi1v[i1].M(), phi2v[i1].M());
804+
const bool dr0 = passDaughterDR(kplus1[i0], kplus2[i0], kminus1[i0], kminus2[i0]);
805+
const bool dr1 = passDaughterDR(kplus1[i1], kplus2[i1], kminus1[i1], kminus2[i1]);
806+
807+
if (sameFour) {
808+
// Choose the candidate closer to mφ (if it passes ΔR)
809+
int keep = (deltam1 < deltam0 && dr1) ? i1 : (dr0 ? i0 : -1);
810+
if (keep >= 0) {
811+
const double dR = deltaR(phi1v[keep].Phi(), phi1v[keep].Eta(), phi2v[keep].Phi(), phi2v[keep].Eta());
812+
const double dm = (keep == i0) ? deltam0 : deltam1;
813+
histos.fill(HIST("SEMassUnlike"), exoticres[keep].M(), exoticres[keep].Pt(), dR, dm, phimult);
814+
}
815+
} else {
816+
// Independent candidates → fill both (respect ΔR)
817+
if (dr0) {
818+
const double dR0 = deltaR(phi1v[i0].Phi(), phi1v[i0].Eta(), phi2v[i0].Phi(), phi2v[i0].Eta());
819+
histos.fill(HIST("SEMassUnlike"), exoticres[i0].M(), exoticres[i0].Pt(), dR0, deltam0, phimult);
820+
}
821+
if (dr1) {
822+
const double dR1 = deltaR(phi1v[i1].Phi(), phi1v[i1].Eta(), phi2v[i1].Phi(), phi2v[i1].Eta());
823+
histos.fill(HIST("SEMassUnlike"), exoticres[i1].M(), exoticres[i1].Pt(), dR1, deltam1, phimult);
824+
}
825+
}
826+
return;
827+
}
828+
829+
// === General case: fill each candidate once (respect ΔR) ===
830+
for (size_t i = 0; i < exoticres.size(); ++i) {
831+
if (!passDaughterDR(kplus1[i], kplus2[i], kminus1[i], kminus2[i]))
832+
continue;
833+
const double dR = deltaR(phi1v[i].Phi(), phi1v[i].Eta(), phi2v[i].Phi(), phi2v[i].Eta());
834+
const double dm = deltam(phi1v[i].M(), phi2v[i].M());
835+
histos.fill(HIST("SEMassUnlike"), exoticres[i].M(), exoticres[i].Pt(), dR, dm, phimult);
836+
}
837+
}
838+
PROCESS_SWITCH(doublephimeson, processopti2, "Process Optimized same event", false);
839+
634840
SliceCache cache;
635841
using BinningTypeVertexContributor = ColumnBinningPolicy<aod::collision::PosZ, aod::collision::NumContrib>;
636842
void processMixedEvent(aod::RedPhiEvents& collisions, aod::PhiTracks& phitracks)

‎PWGLF/Tasks/Resonances/f1protoncorrelation.cxx‎

Lines changed: 12 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -303,13 +303,15 @@ struct f1protoncorrelation {
303303

304304
if (f1track.f1SignalStat() > 0) {
305305
// check charge
306+
float pairCharge = f1track.f1SignalStat() * protontrack.protonCharge();
306307
int f1Charge = f1track.f1SignalStat();
308+
int pionCharge = -1;
309+
int kaonCharge = 1;
307310
if (f1Charge == 2) {
308-
f1Charge = -1;
311+
pionCharge = 1;
312+
kaonCharge = -1;
309313
}
310-
int pionCharge = -1.0 * f1Charge;
311-
float pairCharge = f1Charge * protontrack.protonCharge();
312-
histos.fill(HIST("hPhaseSpaceProtonKaonSame"), Proton.Eta() - Kaon.Eta(), PhiAtSpecificRadiiTPC(Proton, Kaon, protontrack.protonCharge(), f1Charge, bz, bz)); // Phase Space Proton kaon
314+
histos.fill(HIST("hPhaseSpaceProtonKaonSame"), Proton.Eta() - Kaon.Eta(), PhiAtSpecificRadiiTPC(Proton, Kaon, protontrack.protonCharge(), kaonCharge, bz, bz)); // Phase Space Proton kaon
313315
histos.fill(HIST("hPhaseSpaceProtonPionSame"), Proton.Eta() - Kaon.Eta(), PhiAtSpecificRadiiTPC(Proton, Pion, protontrack.protonCharge(), pionCharge, bz, bz)); // Phase Space Proton Pion
314316
histos.fill(HIST("h2SameEventInvariantMassUnlike_mass"), relative_momentum, F1.Pt(), F1.M(), pairCharge); // F1 sign = 1 unlike, F1 sign = -1 like
315317
if (fillSparse) {
@@ -533,14 +535,16 @@ struct f1protoncorrelation {
533535
}
534536
auto relative_momentum = getkstar(F1, Proton);
535537
if (t1.f1SignalStat() > 0) {
538+
float pairCharge = t1.f1SignalStat() * t2.protonCharge();
536539
int f1Charge = t1.f1SignalStat();
540+
int pionCharge = -1;
541+
int kaonCharge = 1;
537542
if (f1Charge == 2) {
538-
f1Charge = -1;
543+
pionCharge = 1;
544+
kaonCharge = -1;
539545
}
540-
int pionCharge = -1.0 * f1Charge;
541-
float pairCharge = f1Charge * t2.protonCharge();
542546
histos.fill(HIST("h2MixEventInvariantMassUnlike_mass"), relative_momentum, F1.Pt(), F1.M(), pairCharge); // F1 sign = 1 unlike, F1 sign = -1 like
543-
histos.fill(HIST("hPhaseSpaceProtonKaonMix"), Proton.Eta() - Kaon.Eta(), PhiAtSpecificRadiiTPC(Proton, Kaon, t2.protonCharge(), f1Charge, bz, bz2)); // Phase Space Proton kaon
547+
histos.fill(HIST("hPhaseSpaceProtonKaonMix"), Proton.Eta() - Kaon.Eta(), PhiAtSpecificRadiiTPC(Proton, Kaon, t2.protonCharge(), kaonCharge, bz, bz2)); // Phase Space Proton kaon
544548
histos.fill(HIST("hPhaseSpaceProtonPionMix"), Proton.Eta() - Kaon.Eta(), PhiAtSpecificRadiiTPC(Proton, Pion, t2.protonCharge(), pionCharge, bz, bz2)); // Phase Space Proton Pion
545549
if (fillSparse) {
546550
histos.fill(HIST("MEMassUnlike"), F1.M(), F1.Pt(), Proton.Pt(), relative_momentum, combinedTPC, pairCharge);

0 commit comments

Comments
 (0)