From 97949204c48339e8b6d634390e032e9fe4caec02 Mon Sep 17 00:00:00 2001 From: sawan Date: Tue, 22 Sep 2026 22:26:37 +0530 Subject: [PATCH 1/4] merged conflicts resolved --- PWGLF/Tasks/Resonances/kstarqa.cxx | 72 +++++++++++++++++++----------- 1 file changed, 46 insertions(+), 26 deletions(-) diff --git a/PWGLF/Tasks/Resonances/kstarqa.cxx b/PWGLF/Tasks/Resonances/kstarqa.cxx index 042672d9c41..28269b2b91a 100644 --- a/PWGLF/Tasks/Resonances/kstarqa.cxx +++ b/PWGLF/Tasks/Resonances/kstarqa.cxx @@ -224,8 +224,8 @@ struct Kstarqa { rEventSelection.add("hEventCut", "No. of event after cuts", kTH1D, {{20, 0, 20}}); std::shared_ptr hEventSelection = rEventSelection.get(HIST("hEventCut")); - rEventSelection.add("hTrackCut", "No. of tracks after cuts", kTH1D, {{20, 0, 20}}); - std::shared_ptr hTrackSelection = rEventSelection.get(HIST("hTrackCut")); + rTrackSelection.add("hTrackCut", "No. of tracks after cuts", kTH1D, {{20, 0, 20}}); + std::shared_ptr hTrackSelection = rTrackSelection.get(HIST("hTrackCut")); auto check = [](bool enabled) { return enabled ? "" : " #otimes"; }; // check if a cut is enabled and put #otimes beside that label if not enabled @@ -249,8 +249,8 @@ struct Kstarqa { std::vector trackCutLabels = { "All Tracks", std::string("GlobalTracks") + check(configGp.isGlobalTracks.value), - std::string("pT > ") + std::to_string(configGp.cfgCutPT), - std::string("|#eta| < ") + std::to_string(configGp.cfgCutEtaMax), + std::string("pT > ") + std::to_string(configGp.cfgCutPT.value), + std::string("|#eta| < ") + std::to_string(configGp.cfgCutEtaMax.value), std::string("DCAxy < ") + std::to_string(configGp.cfgCutDCAxyMax.value) + check(!configGp.isApplyPtDepDCACut.value), std::string("DCAz < ") + std::to_string(configGp.cfgCutDCAz.value), std::string("ITS clusters > ") + std::to_string(configGp.cfgITScluster.value), @@ -506,24 +506,27 @@ struct Kstarqa { template bool selectionTrack(const T& candidate) { - - rEventSelection.fill(HIST("hTrackCut"), 0); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 0); if (configGp.isGlobalTracks) { if (!candidate.isGlobalTrack()) return false; - rEventSelection.fill(HIST("hTrackCut"), 1); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 1); if (std::abs(candidate.pt()) < configGp.cfgCutPT) return false; - rEventSelection.fill(HIST("hTrackCut"), 2); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 2); if (std::abs(candidate.eta()) > configGp.cfgCutEtaMax) return false; - rEventSelection.fill(HIST("hTrackCut"), 3); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 3); if (!configGp.isApplyPtDepDCACut) { if (std::abs(candidate.dcaXY()) > configGp.cfgCutDCAxyMax) @@ -534,7 +537,8 @@ struct Kstarqa { return false; } - rEventSelection.fill(HIST("hTrackCut"), 4); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 4); if (!configGp.isApplyPtDepDCACut) { if (std::abs(candidate.dcaZ()) > configGp.cfgCutDCAz) @@ -544,22 +548,26 @@ struct Kstarqa { return false; } - rEventSelection.fill(HIST("hTrackCut"), 5); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 5); if (candidate.itsNCls() < configGp.cfgITScluster) return false; - rEventSelection.fill(HIST("hTrackCut"), 6); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 6); if (candidate.tpcNClsFound() < configGp.cfgTPCcluster) return false; - rEventSelection.fill(HIST("hTrackCut"), 7); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 7); if (configGp.hasITS && !candidate.hasITS()) return false; - rEventSelection.fill(HIST("hTrackCut"), 8); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 8); if (configGp.isITSTPCRefit) { if (!(candidate.flags() & o2::aod::track::ITSrefit) || @@ -568,31 +576,37 @@ struct Kstarqa { } } - rEventSelection.fill(HIST("hTrackCut"), 9); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 9); if (configGp.cfgPVContributor && !candidate.isPVContributor()) return false; - rEventSelection.fill(HIST("hTrackCut"), 10); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 10); } else if (!configGp.isGlobalTracks) { - rEventSelection.fill(HIST("hTrackCut"), 0); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 0); if (!candidate.isGlobalTrackWoDCA()) return false; - rEventSelection.fill(HIST("hTrackCut"), 1); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 1); if (std::abs(candidate.pt()) < configGp.cfgCutPT) return false; - rEventSelection.fill(HIST("hTrackCut"), 2); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 2); if (std::abs(candidate.eta()) > configGp.cfgCutEtaMax) return false; - rEventSelection.fill(HIST("hTrackCut"), 3); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 3); if (!configGp.isApplyPtDepDCACut) { if (std::abs(candidate.dcaXY()) > configGp.cfgCutDCAxyMax) @@ -603,7 +617,8 @@ struct Kstarqa { return false; } - rEventSelection.fill(HIST("hTrackCut"), 4); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 4); if (!configGp.isApplyPtDepDCACut) { if (std::abs(candidate.dcaZ()) > configGp.cfgCutDCAz) @@ -613,12 +628,14 @@ struct Kstarqa { return false; } - rEventSelection.fill(HIST("hTrackCut"), 5); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 5); if (candidate.itsNCls() < configGp.cfgITScluster) return false; - rEventSelection.fill(HIST("hTrackCut"), 6); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 6); if (candidate.tpcNClsFound() < configGp.cfgTPCcluster) return false; @@ -626,7 +643,8 @@ struct Kstarqa { if (configGp.hasITS && !candidate.hasITS()) return false; - rEventSelection.fill(HIST("hTrackCut"), 7); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 7); if (configGp.isITSTPCRefit) { if (!(candidate.flags() & o2::aod::track::ITSrefit) || @@ -635,12 +653,14 @@ struct Kstarqa { } } - rEventSelection.fill(HIST("hTrackCut"), 8); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 8); if (configGp.cfgPVContributor && !candidate.isPVContributor()) return false; - rEventSelection.fill(HIST("hTrackCut"), 9); + if (fillHist) + rEventSelection.fill(HIST("hTrackCut"), 9); } return true; From f4b8fc8dab828a2f5ab9dba6ab3786f025093a26 Mon Sep 17 00:00:00 2001 From: sawan Date: Fri, 18 Sep 2026 15:49:01 +0530 Subject: [PATCH 2/4] fixed errors --- PWGLF/Tasks/Resonances/kstarqa.cxx | 72 +++++++++++------------------- 1 file changed, 26 insertions(+), 46 deletions(-) diff --git a/PWGLF/Tasks/Resonances/kstarqa.cxx b/PWGLF/Tasks/Resonances/kstarqa.cxx index 28269b2b91a..042672d9c41 100644 --- a/PWGLF/Tasks/Resonances/kstarqa.cxx +++ b/PWGLF/Tasks/Resonances/kstarqa.cxx @@ -224,8 +224,8 @@ struct Kstarqa { rEventSelection.add("hEventCut", "No. of event after cuts", kTH1D, {{20, 0, 20}}); std::shared_ptr hEventSelection = rEventSelection.get(HIST("hEventCut")); - rTrackSelection.add("hTrackCut", "No. of tracks after cuts", kTH1D, {{20, 0, 20}}); - std::shared_ptr hTrackSelection = rTrackSelection.get(HIST("hTrackCut")); + rEventSelection.add("hTrackCut", "No. of tracks after cuts", kTH1D, {{20, 0, 20}}); + std::shared_ptr hTrackSelection = rEventSelection.get(HIST("hTrackCut")); auto check = [](bool enabled) { return enabled ? "" : " #otimes"; }; // check if a cut is enabled and put #otimes beside that label if not enabled @@ -249,8 +249,8 @@ struct Kstarqa { std::vector trackCutLabels = { "All Tracks", std::string("GlobalTracks") + check(configGp.isGlobalTracks.value), - std::string("pT > ") + std::to_string(configGp.cfgCutPT.value), - std::string("|#eta| < ") + std::to_string(configGp.cfgCutEtaMax.value), + std::string("pT > ") + std::to_string(configGp.cfgCutPT), + std::string("|#eta| < ") + std::to_string(configGp.cfgCutEtaMax), std::string("DCAxy < ") + std::to_string(configGp.cfgCutDCAxyMax.value) + check(!configGp.isApplyPtDepDCACut.value), std::string("DCAz < ") + std::to_string(configGp.cfgCutDCAz.value), std::string("ITS clusters > ") + std::to_string(configGp.cfgITScluster.value), @@ -506,27 +506,24 @@ struct Kstarqa { template bool selectionTrack(const T& candidate) { - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 0); + + rEventSelection.fill(HIST("hTrackCut"), 0); if (configGp.isGlobalTracks) { if (!candidate.isGlobalTrack()) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 1); + rEventSelection.fill(HIST("hTrackCut"), 1); if (std::abs(candidate.pt()) < configGp.cfgCutPT) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 2); + rEventSelection.fill(HIST("hTrackCut"), 2); if (std::abs(candidate.eta()) > configGp.cfgCutEtaMax) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 3); + rEventSelection.fill(HIST("hTrackCut"), 3); if (!configGp.isApplyPtDepDCACut) { if (std::abs(candidate.dcaXY()) > configGp.cfgCutDCAxyMax) @@ -537,8 +534,7 @@ struct Kstarqa { return false; } - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 4); + rEventSelection.fill(HIST("hTrackCut"), 4); if (!configGp.isApplyPtDepDCACut) { if (std::abs(candidate.dcaZ()) > configGp.cfgCutDCAz) @@ -548,26 +544,22 @@ struct Kstarqa { return false; } - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 5); + rEventSelection.fill(HIST("hTrackCut"), 5); if (candidate.itsNCls() < configGp.cfgITScluster) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 6); + rEventSelection.fill(HIST("hTrackCut"), 6); if (candidate.tpcNClsFound() < configGp.cfgTPCcluster) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 7); + rEventSelection.fill(HIST("hTrackCut"), 7); if (configGp.hasITS && !candidate.hasITS()) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 8); + rEventSelection.fill(HIST("hTrackCut"), 8); if (configGp.isITSTPCRefit) { if (!(candidate.flags() & o2::aod::track::ITSrefit) || @@ -576,37 +568,31 @@ struct Kstarqa { } } - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 9); + rEventSelection.fill(HIST("hTrackCut"), 9); if (configGp.cfgPVContributor && !candidate.isPVContributor()) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 10); + rEventSelection.fill(HIST("hTrackCut"), 10); } else if (!configGp.isGlobalTracks) { - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 0); + rEventSelection.fill(HIST("hTrackCut"), 0); if (!candidate.isGlobalTrackWoDCA()) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 1); + rEventSelection.fill(HIST("hTrackCut"), 1); if (std::abs(candidate.pt()) < configGp.cfgCutPT) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 2); + rEventSelection.fill(HIST("hTrackCut"), 2); if (std::abs(candidate.eta()) > configGp.cfgCutEtaMax) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 3); + rEventSelection.fill(HIST("hTrackCut"), 3); if (!configGp.isApplyPtDepDCACut) { if (std::abs(candidate.dcaXY()) > configGp.cfgCutDCAxyMax) @@ -617,8 +603,7 @@ struct Kstarqa { return false; } - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 4); + rEventSelection.fill(HIST("hTrackCut"), 4); if (!configGp.isApplyPtDepDCACut) { if (std::abs(candidate.dcaZ()) > configGp.cfgCutDCAz) @@ -628,14 +613,12 @@ struct Kstarqa { return false; } - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 5); + rEventSelection.fill(HIST("hTrackCut"), 5); if (candidate.itsNCls() < configGp.cfgITScluster) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 6); + rEventSelection.fill(HIST("hTrackCut"), 6); if (candidate.tpcNClsFound() < configGp.cfgTPCcluster) return false; @@ -643,8 +626,7 @@ struct Kstarqa { if (configGp.hasITS && !candidate.hasITS()) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 7); + rEventSelection.fill(HIST("hTrackCut"), 7); if (configGp.isITSTPCRefit) { if (!(candidate.flags() & o2::aod::track::ITSrefit) || @@ -653,14 +635,12 @@ struct Kstarqa { } } - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 8); + rEventSelection.fill(HIST("hTrackCut"), 8); if (configGp.cfgPVContributor && !candidate.isPVContributor()) return false; - if (fillHist) - rEventSelection.fill(HIST("hTrackCut"), 9); + rEventSelection.fill(HIST("hTrackCut"), 9); } return true; From c262063fda0d5edea679611d49cfc9aa5e5b0d40 Mon Sep 17 00:00:00 2001 From: sawan Date: Tue, 22 Sep 2026 22:22:38 +0530 Subject: [PATCH 3/4] Added process function for cross checks and mass shift --- PWGLF/Tasks/Resonances/doublephimeson.cxx | 825 ++++++++++++++++++++++ 1 file changed, 825 insertions(+) diff --git a/PWGLF/Tasks/Resonances/doublephimeson.cxx b/PWGLF/Tasks/Resonances/doublephimeson.cxx index b027bfdd2d3..beb670129ab 100644 --- a/PWGLF/Tasks/Resonances/doublephimeson.cxx +++ b/PWGLF/Tasks/Resonances/doublephimeson.cxx @@ -27,6 +27,9 @@ #include #include +#include "TLorentzVector.h" +#include "TMath.h" +#include "TMatrixD.h" #include #include #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) @@ -34,6 +37,7 @@ #include #include #include +#include #include #include @@ -75,6 +79,8 @@ struct doublephimeson { Configurable cfgCrossPhiHigh{"cfgCrossPhiHigh", 1.03, "Upper edge of phi mass window for cross-pairing (ghost) veto"}; Configurable useParametrized{"useParametrized", false, "Use pT dependent mass peak and width"}; Configurable useCrossPairRejection{"useCrossPairRejection", true, "Use cross pair phi signal compatibilaty"}; + Configurable cRotations{"cRotations", 10, "Number of rotations for rotational background"}; + Configurable applyMomentumShift{"applyMomentumShift", false, "Apply momentum shift to kaons to check effect on phi mass peak"}; Configurable cfgFillDataDrivenPhiResolution{ "cfgFillDataDrivenPhiResolution", true, "Fill the single-phi daughter-kinematics sparse used for data-driven X resolution and inter-dataset momentum-scale calibration"}; @@ -224,6 +230,7 @@ struct doublephimeson { ConfigurableAxis configThnAxisPhiPtVertex{"configThnAxisPhiPtVertex", {100, 0.0, 100.0}, "phi pT (GeV/c)"}; ConfigurableAxis configThnAxisDecayLength{"configThnAxisDecayLength", {200, 0.0, 1.0}, "3D decay length (cm)"}; ConfigurableAxis configThnAxisFitChi2Ndf{"configThnAxisFitChi2Ndf", {200, 0.0, 100.0}, "four-kaon fit chi2/NDF"}; + ConfigurableAxis configThnAxisFitProbability{"configThnAxisFitProbability", {100, 0.0, 1.0}, "four-kaon fit probability"}; ConfigurableAxis configThnAxisRmsDcaSig{"configThnAxisRmsDcaSig", {300, 0.0, 15.0}, "RMS DCA significance"}; // Data-driven mass-resolution inputs. @@ -289,7 +296,9 @@ struct doublephimeson { histos.add("hnsigmaTPCKaonMinus", "hnsigmaTPCKaonMinus", kTH2F, {{1000, -3.0, 3.0f}, {100, 0.0f, 10.0f}}); histos.add("hnsigmaTPCTOFKaon", "hnsigmaTPCTOFKaon", kTH3F, {{500, -3.0, 3.0f}, {500, -3.0, 3.0f}, {100, 0.0f, 10.0f}}); histos.add("hPhiMassVsPt", "hPhiMassVsPt", kTH2F, {{40, 1.0, 1.04f}, {1000, 0.0f, 100.0f}}); + histos.add("hPhiMassVsPtShifted", "hPhiMassVsPtShifted", kTH2F, {{40, 1.0, 1.04f}, {1000, 0.0f, 100.0f}}); histos.add("hPhiMass", "hPhiMass", kTH3F, {{40, 1.0, 1.04f}, {40, 1.0, 1.04f}, {250, 0.0f, 100.0f}}); + histos.add("hPhiMassShifted", "hPhiMassShifted", kTH3F, {{40, 1.0, 1.04f}, {40, 1.0, 1.04f}, {250, 0.0f, 100.0f}}); histos.add("hPhiMassNormalized", "hPhiMassNormalized", kTH3F, {{100, -10.0, 10.0f}, {100, -10.0, 10.0f}, {250, 0.0f, 100.0f}}); histos.add("hPhiMass2", "hPhiMass2", kTH2F, {{40, 1.0, 1.04f}, {40, 1.0f, 1.04f}}); histos.add("hkPlusDeltaetaDeltaPhi", "hkPlusDeltaetaDeltaPhi", kTH2F, {{400, -2.0, 2.0}, {640, -2.0 * TMath::Pi(), 2.0 * TMath::Pi()}}); @@ -316,6 +325,7 @@ struct doublephimeson { AxisSpec axisDoublePhiPID{50, 0.0, 5.0, "max daughter n_{#sigma}^{comb}"}; const AxisSpec thnAxisDecayLength{configThnAxisDecayLength, "#it{L}_{3D} (cm)"}; const AxisSpec thnAxisFitChi2Ndf{configThnAxisFitChi2Ndf, "#chi^{2}/NDF"}; + const AxisSpec thnAxisFitProbability{configThnAxisFitProbability, "fit probability"}; const AxisSpec thnAxisRmsDcaSig{configThnAxisRmsDcaSig, "RMS DCA significance"}; const AxisSpec ddPhiMassAxis{cfgDDPhiMassAxis, "m_{K^{+}K^{-}} (GeV/c^{2})"}; @@ -366,6 +376,45 @@ struct doublephimeson { thnAxisNumPhi, axisDoublePhiPID}); + histos.add("SEMassPhiPhi", "SEMassPhiPhi", HistType::kTHnSparseF, + { + thnAxisInvMass, // M(phi-phi) + thnAxisPt, // pT(phi-phi) + thnAxisInvMassDeltaPhi // DeltaM_phi + }); + + histos.add("SEMassPhiPhiRefitted", "SEMassPhiPhiRefitted", HistType::kTHnSparseF, + { + thnAxisInvMass, // M(phi-phi) + thnAxisPt, // pT(phi-phi) + thnAxisInvMassDeltaPhi, // DeltaM_phi + thnAxisFitChi2Ndf, // chi2/NDF of the 4-kaon kinematic fit + thnAxisFitProbability, // fit probability of the 4-kaon kinematic fit + thnAxisInvMassPhi, // m(phi1) + thnAxisInvMassPhi // m(phi2) + }); + + histos.add("SEMassPhiPhiShifted", "SEMassPhiPhiShifted", HistType::kTHnSparseF, + { + thnAxisInvMass, // M(phi-phi) + thnAxisPt, // pT(phi-phi) + thnAxisInvMassDeltaPhi, // DeltaM_phi + thnAxisFitChi2Ndf, // chi2/NDF of the 4-kaon kinematic fit + thnAxisFitProbability, // fit probability of the 4-kaon kinematic fit + thnAxisInvMassPhi, // m(phi1) + thnAxisInvMassPhi // m(phi2) + }); + + histos.add("SEMassPhiPhiRotational", "SEMassPhiPhiRotational", HistType::kTHnSparseF, + { + thnAxisInvMass, // M(phi-phi) + thnAxisPt // pT(phi-phi) + }); + + histos.add("SEMassUnlike_VertexVars", "SEMassUnlike_VertexVars", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisInvMassDeltaPhi, thnAxisInvMassPhi, thnAxisInvMassPhi, thnAxisDecayLength, thnAxisFitChi2Ndf, thnAxisRmsDcaSig}); + histos.add("NPhiPerEvent", "NPhiPerEvent", HistType::kTH1F, {{20, 0, 20}}); + histos.add("NEvents", "NEvents", HistType::kTH1F, {{2, 0, 2}}); + histos.add("SEMassUnlike_VertexVars", "SEMassUnlike_VertexVars", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisInvMassDeltaPhi, thnAxisInvMassPhi, thnAxisInvMassPhi, thnAxisDecayLength, thnAxisFitChi2Ndf, thnAxisRmsDcaSig}); // Single-phi calibration input. Axis order: @@ -402,6 +451,7 @@ struct doublephimeson { ddCalibKaonPtAxis, ddCalibKaonPtAxis, ddCalibKaonPtAxis, ddCalibKaonPtAxis}); } + TRandom* rn = new TRandom(); // get kstar TLorentzVector trackSum, PartOneCMS, PartTwoCMS, trackRelK; @@ -994,6 +1044,252 @@ struct doublephimeson { return false; } + struct FourKFitResult { + bool converged = false; + + double chi2 = -1.0; + int ndf = 2; + double probability = -1.0; + + double refittedMass = -1.0; + + TLorentzVector k11; + TLorentzVector k12; + TLorentzVector k21; + TLorentzVector k22; + + double phi1Mass = -1.0; + double phi2Mass = -1.0; + }; + + FourKFitResult fitFourKaons( + const TLorentzVector& k11Init, + const TLorentzVector& k12Init, + const TLorentzVector& k21Init, + const TLorentzVector& k22Init, + double mK, + double mPhi, + double sigmaPFrac = 0.01, + int maxIterations = 20, + double tolerance = 1e-6) + { + FourKFitResult result; + + constexpr int N = 12; + constexpr int NC = 2; + + // ============================================================ + // Initial state vector + // + // x = (px1,py1,pz1, + // px2,py2,pz2, + // px3,py3,pz3, + // px4,py4,pz4) + // ============================================================ + + TMatrixD x(N, 1); + TMatrixD x0(N, 1); + + const TLorentzVector* pInit[4] = {&k11Init, &k12Init, &k21Init, &k22Init}; + for (int i = 0; i < 4; ++i) { + + x(3 * i + 0, 0) = pInit[i]->Px(); + x(3 * i + 1, 0) = pInit[i]->Py(); + x(3 * i + 2, 0) = pInit[i]->Pz(); + + x0(3 * i + 0, 0) = pInit[i]->Px(); + x0(3 * i + 1, 0) = pInit[i]->Py(); + x0(3 * i + 2, 0) = pInit[i]->Pz(); + } + + // ============================================================ + // Approximate covariance matrix + // + // This is ONLY for Stage 1/2. + // Later replace this with the actual O2 track covariance. + // ============================================================ + + TMatrixD V(N, N); + V.Zero(); + + for (int i = 0; i < 4; ++i) { + + const double p = pInit[i]->P(); + + const double sigmaP = sigmaPFrac * std::max(p, 1e-6); + + V(3 * i + 0, 3 * i + 0) = sigmaP * sigmaP; + V(3 * i + 1, 3 * i + 1) = sigmaP * sigmaP; + V(3 * i + 2, 3 * i + 2) = sigmaP * sigmaP; + } + + TMatrixD Vinv = V; + Vinv.Invert(); + + // ============================================================ + // Helper to construct a kaon 4-vector from x + // ============================================================ + + auto makeLV = [&](const TMatrixD& xx, int particle) { + const double px = xx(3 * particle + 0, 0); + const double py = xx(3 * particle + 1, 0); + const double pz = xx(3 * particle + 2, 0); + + TLorentzVector p; + p.SetXYZM(px, py, pz, mK); + + return p; + }; + + // ============================================================ + // Iterative constrained fit + // ============================================================ + + for (int iteration = 0; iteration < maxIterations; ++iteration) { + + TLorentzVector k11 = makeLV(x, 0); + TLorentzVector k12 = makeLV(x, 1); + TLorentzVector k21 = makeLV(x, 2); + TLorentzVector k22 = makeLV(x, 3); + TLorentzVector phi1 = k11 + k12; + TLorentzVector phi2 = k21 + k22; + // ---------------------------------------------------------- + // Constraints + // + // f1 = M(phi1)^2 - Mphi^2 + // f2 = M(phi2)^2 - Mphi^2 + // ---------------------------------------------------------- + + TMatrixD f(NC, 1); + + f(0, 0) = phi1.M2() - mPhi * mPhi; + + f(1, 0) = phi2.M2() - mPhi * mPhi; + + // ---------------------------------------------------------- + // Check whether already sufficiently close + // ---------------------------------------------------------- + + const double scale = mPhi * mPhi; + + if (std::abs(f(0, 0)) < tolerance * scale && + std::abs(f(1, 0)) < tolerance * scale) { + + result.converged = true; + break; + } + + // ---------------------------------------------------------- + // Numerical Jacobian + // + // A(i,j) = df_i / dx_j + // ---------------------------------------------------------- + + TMatrixD A(NC, N); + A.Zero(); + + for (int j = 0; j < N; ++j) { + + TMatrixD xp = x; + TMatrixD xm = x; + + const double step = 1e-5 * std::max(std::abs(x(j, 0)), 1.0); + + xp(j, 0) += step; + xm(j, 0) -= step; + + TLorentzVector p1p = makeLV(xp, 0); + TLorentzVector p2p = makeLV(xp, 1); + TLorentzVector p3p = makeLV(xp, 2); + TLorentzVector p4p = makeLV(xp, 3); + TLorentzVector p1m = makeLV(xm, 0); + TLorentzVector p2m = makeLV(xm, 1); + TLorentzVector p3m = makeLV(xm, 2); + TLorentzVector p4m = makeLV(xm, 3); + + const double f1p = (p1p + p2p).M2() - mPhi * mPhi; + const double f1m = (p1m + p2m).M2() - mPhi * mPhi; + const double f2p = (p3p + p4p).M2() - mPhi * mPhi; + const double f2m = (p3m + p4m).M2() - mPhi * mPhi; + A(0, j) = (f1p - f1m) / (2.0 * step); + A(1, j) = (f2p - f2m) / (2.0 * step); + } + + // ---------------------------------------------------------- + // C = A V A^T + // ---------------------------------------------------------- + + TMatrixD AT(TMatrixD::kTransposed, A); + TMatrixD C = A * V * AT; + TMatrixD Cinv = C; + Cinv.Invert(); + + // ---------------------------------------------------------- + // Delta x + // + // dx = -V A^T (A V A^T)^-1 f + // ---------------------------------------------------------- + + TMatrixD dx = -1.0 * V * AT * Cinv * f; + + x += dx; + + // ---------------------------------------------------------- + // Check convergence + // ---------------------------------------------------------- + + double maxCorrection = 0.0; + + for (int j = 0; j < N; ++j) { + + maxCorrection = std::max(maxCorrection, std::abs(dx(j, 0))); + } + + if (maxCorrection < tolerance) { + + result.converged = true; + break; + } + } + + // ============================================================ + // Construct final particles + // ============================================================ + + result.k11 = makeLV(x, 0); + result.k12 = makeLV(x, 1); + result.k21 = makeLV(x, 2); + result.k22 = makeLV(x, 3); + + TLorentzVector phi1Fit = result.k11 + result.k12; + TLorentzVector phi2Fit = result.k21 + result.k22; + TLorentzVector pairFit = result.k11 + result.k12 + result.k21 + result.k22; + + result.phi1Mass = phi1Fit.M(); + result.phi2Mass = phi2Fit.M(); + result.refittedMass = pairFit.M(); + // ============================================================ + // chi2 + // + // chi2 = (x-x0)^T V^-1 (x-x0) + // ============================================================ + + TMatrixD deltaX = x; + deltaX -= x0; + TMatrixD deltaXT(TMatrixD::kTransposed, deltaX); + TMatrixD chi2Matrix = deltaXT * Vinv * deltaX; + + result.chi2 = chi2Matrix(0, 0); + result.ndf = NC; + + if (result.converged) { + + result.probability = TMath::Prob(result.chi2, result.ndf); + } + + return result; + } + TLorentzVector exotic, Phid1, Phid2; TLorentzVector Phi1kaonplus, Phi1kaonminus, Phi2kaonplus, Phi2kaonminus; // TLorentzVector exoticRot, Phid1Rot; @@ -2322,6 +2618,535 @@ struct doublephimeson { } PROCESS_SWITCH(doublephimeson, processPairOpti6, "Process fitted phi-phi pairs with vertex variables", false); + void processOpti7(aod::RedPhiEvents::iterator const& collision, aod::PhiTracks const& phitracks) + { + if (additionalEvsel && (collision.numPos() < 2 || collision.numNeg() < 2)) { + return; + } + constexpr double mPhiPDG = o2::constants::physics::MassPhi; + constexpr double mKPDG = o2::constants::physics::MassKPlus; + int phimult = 0; + + for (const auto& Phitrack : phitracks) { + const double kpluspt = std::hypot(Phitrack.phid1Px(), Phitrack.phid1Py()); + const double kminuspt = std::hypot(Phitrack.phid2Px(), Phitrack.phid2Py()); + + // pT cut on kaon tracks + if (kpluspt > maxKaonPt || kminuspt > maxKaonPt) { + continue; + } + + // Mass window selection for phi candidates + if (Phitrack.phiMass() < minPhiMass1 || Phitrack.phiMass() > maxPhiMass1) { + continue; + } + + TLorentzVector phi; + phi.SetXYZM(Phitrack.phiPx(), Phitrack.phiPy(), Phitrack.phiPz(), Phitrack.phiMass()); + + // pT cut on phi daughters + if (phi.Pt() < minPhiPt || phi.Pt() > maxPhiPt) { + continue; + } + + // PID selection for kaon tracks + if (!selectionPID(Phitrack.phid1TPC(), Phitrack.phid1TOF(), Phitrack.phid1TOFHit(), strategyPID1, kpluspt)) { + continue; + } + if (!selectionPID(Phitrack.phid2TPC(), Phitrack.phid2TOF(), Phitrack.phid2TOFHit(), strategyPID2, kminuspt)) { + continue; + } + phimult++; + } + + if (phimult < 2) { + return; + } + + for (auto const& Phitrack1 : phitracks) { + for (auto const& Phitrack2 : phitracks) { + + // Avoid double counting + if (Phitrack2.index() <= Phitrack1.index()) { + continue; + } + + // pT cut kaon tracks + const double kplus1pt = std::hypot(Phitrack1.phid1Px(), Phitrack1.phid1Py()); + const double kminus1pt = std::hypot(Phitrack1.phid2Px(), Phitrack1.phid2Py()); + const double kplus2pt = std::hypot(Phitrack2.phid1Px(), Phitrack2.phid1Py()); + const double kminus2pt = std::hypot(Phitrack2.phid2Px(), Phitrack2.phid2Py()); + + if (kplus1pt > maxKaonPt || kminus1pt > maxKaonPt || kplus2pt > maxKaonPt || kminus2pt > maxKaonPt) { + continue; + } + + // Mass window selection for phi candidates + if (Phitrack1.phiMass() < minPhiMass1 || Phitrack1.phiMass() > maxPhiMass1 || Phitrack2.phiMass() < minPhiMass1 || Phitrack2.phiMass() > maxPhiMass1) { + continue; + } + + // pT cut on phi daughters + TLorentzVector phi1, phi2; + phi1.SetXYZM(Phitrack1.phiPx(), Phitrack1.phiPy(), Phitrack1.phiPz(), Phitrack1.phiMass()); + phi2.SetXYZM(Phitrack2.phiPx(), Phitrack2.phiPy(), Phitrack2.phiPz(), Phitrack2.phiMass()); + if (phi1.Pt() < minPhiPt || phi1.Pt() > maxPhiPt || phi2.Pt() < minPhiPt || phi2.Pt() > maxPhiPt) { + continue; + } + + // PID selection for kaon tracks + if (!selectionPID(Phitrack1.phid1TPC(), Phitrack1.phid1TOF(), Phitrack1.phid1TOFHit(), strategyPID1, kplus1pt) || + !selectionPID(Phitrack1.phid2TPC(), Phitrack1.phid2TOF(), Phitrack1.phid2TOFHit(), strategyPID2, kminus1pt) || + !selectionPID(Phitrack2.phid1TPC(), Phitrack2.phid1TOF(), Phitrack2.phid1TOFHit(), strategyPID1, kplus2pt) || + !selectionPID(Phitrack2.phid2TPC(), Phitrack2.phid2TOF(), Phitrack2.phid2TOFHit(), strategyPID2, kminus2pt)) { + continue; + } + + // Check for shared daughters + if (Phitrack1.phid1Index() == Phitrack2.phid1Index() || + Phitrack1.phid1Index() == Phitrack2.phid2Index() || + Phitrack1.phid2Index() == Phitrack2.phid1Index() || + Phitrack1.phid2Index() == Phitrack2.phid2Index()) { + continue; + } + + TLorentzVector pair = phi1 + phi2; + // Mass window range for the phi-phi pair + if (pair.Pt() < minExoticPt || pair.M() < minExoticMass || pair.M() > maxExoticMass) { + continue; + } + + double deltaM = std::hypot(Phitrack1.phiMass() - mPhiPDG, Phitrack2.phiMass() - mPhiPDG); + + histos.fill(HIST("SEMassPhiPhi"), + pair.M(), + pair.Pt(), + deltaM); + } + } + } + PROCESS_SWITCH(doublephimeson, processOpti7, "Process optimised save-event for cross-checks", false); + + void processOpti8(aod::RedPhiEvents::iterator const& collision, aod::PhiTracks const& phitracks) + { + histos.fill(HIST("NEvents"), 0.5); + if (additionalEvsel && (collision.numPos() < 2 || collision.numNeg() < 2)) { + return; + } + histos.fill(HIST("NEvents"), 1.5); + constexpr double mPhiPDG = o2::constants::physics::MassPhi; + constexpr double mKPDG = o2::constants::physics::MassKPlus; + int phimult = 0; + + for (const auto& Phitrack : phitracks) { + const double kpluspt = std::hypot(Phitrack.phid1Px(), Phitrack.phid1Py()); + const double kminuspt = std::hypot(Phitrack.phid2Px(), Phitrack.phid2Py()); + + histos.fill(HIST("hnsigmaTPCTOFKaonBefore"), Phitrack.phid1TPC(), Phitrack.phid1TOF(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonPlusBefore"), Phitrack.phid1TPC(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonMinusBefore"), Phitrack.phid2TPC(), kminuspt); + + // pT cut on kaon tracks + if (kpluspt > maxKaonPt || kminuspt > maxKaonPt) { + continue; + } + + // Mass window selection for phi candidates + if (Phitrack.phiMass() < minPhiMass1 || Phitrack.phiMass() > maxPhiMass1) { + continue; + } + + TLorentzVector phi; + phi.SetXYZM(Phitrack.phiPx(), Phitrack.phiPy(), Phitrack.phiPz(), Phitrack.phiMass()); + + // pT cut on phi daughters + if (phi.Pt() < minPhiPt || phi.Pt() > maxPhiPt) { + continue; + } + + // PID selection for kaon tracks + if (!selectionPID(Phitrack.phid1TPC(), Phitrack.phid1TOF(), Phitrack.phid1TOFHit(), strategyPID1, kpluspt)) { + continue; + } + if (!selectionPID(Phitrack.phid2TPC(), Phitrack.phid2TOF(), Phitrack.phid2TOFHit(), strategyPID2, kminuspt)) { + continue; + } + + histos.fill(HIST("hnsigmaTPCTOFKaon"), Phitrack.phid1TPC(), Phitrack.phid1TOF(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonPlus"), Phitrack.phid1TPC(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonMinus"), Phitrack.phid2TPC(), kminuspt); + + phimult++; + } + + histos.fill(HIST("NPhiPerEvent"), phimult); + + if (phimult < 2) { + return; + } + + for (auto const& Phitrack1 : phitracks) { + + // pT cut kaon tracks + const double kplus1pt = std::hypot(Phitrack1.phid1Px(), Phitrack1.phid1Py()); + const double kminus1pt = std::hypot(Phitrack1.phid2Px(), Phitrack1.phid2Py()); + if (kplus1pt > maxKaonPt || kminus1pt > maxKaonPt) { + continue; + } + + // PID selection for kaon tracks + if (!selectionPID(Phitrack1.phid1TPC(), Phitrack1.phid1TOF(), Phitrack1.phid1TOFHit(), strategyPID1, kplus1pt) || + !selectionPID(Phitrack1.phid2TPC(), Phitrack1.phid2TOF(), Phitrack1.phid2TOFHit(), strategyPID2, kminus1pt)) { + continue; + } + + TLorentzVector phi1; + phi1.SetXYZM(Phitrack1.phiPx(), Phitrack1.phiPy(), Phitrack1.phiPz(), Phitrack1.phiMass()); + if (phi1.Pt() < minPhiPt || phi1.Pt() > maxPhiPt) { + continue; + } + + // Mass window selection for phi candidates + if (Phitrack1.phiMass() < minPhiMass1 || Phitrack1.phiMass() > maxPhiMass1) { + continue; + } + + histos.fill(HIST("hPhiMassVsPt"), phi1.M(), phi1.Pt()); + + for (auto const& Phitrack2 : phitracks) { + + // Avoid double counting + if (Phitrack2.index() <= Phitrack1.index()) { + continue; + } + + // pT cut kaon tracks + const double kplus2pt = std::hypot(Phitrack2.phid1Px(), Phitrack2.phid1Py()); + const double kminus2pt = std::hypot(Phitrack2.phid2Px(), Phitrack2.phid2Py()); + + if (kplus2pt > maxKaonPt || kminus2pt > maxKaonPt) { + continue; + } + + // Mass window selection for phi candidates + if (Phitrack2.phiMass() < minPhiMass1 || Phitrack2.phiMass() > maxPhiMass1) { + continue; + } + + // pT cut on phi + TLorentzVector phi2; + phi2.SetXYZM(Phitrack2.phiPx(), Phitrack2.phiPy(), Phitrack2.phiPz(), Phitrack2.phiMass()); + if (phi2.Pt() < minPhiPt || phi2.Pt() > maxPhiPt) { + continue; + } + + // PID selection for kaon tracks + if (!selectionPID(Phitrack2.phid1TPC(), Phitrack2.phid1TOF(), Phitrack2.phid1TOFHit(), strategyPID1, kplus2pt) || + !selectionPID(Phitrack2.phid2TPC(), Phitrack2.phid2TOF(), Phitrack2.phid2TOFHit(), strategyPID2, kminus2pt)) { + continue; + } + + // Check for shared daughters + if (Phitrack1.phid1Index() == Phitrack2.phid1Index() || + Phitrack1.phid1Index() == Phitrack2.phid2Index() || + Phitrack1.phid2Index() == Phitrack2.phid1Index() || + Phitrack1.phid2Index() == Phitrack2.phid2Index()) { + continue; + } + + for (int i = 0; i < cRotations; i++) { + double thetaRot = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / 10, o2::constants::math::PI + o2::constants::math::PI / 10); + + TLorentzVector daughterRot; + daughterRot.SetXYZM(phi1.Px() * std::cos(thetaRot) - phi1.Py() * std::sin(thetaRot), phi1.Px() * std::sin(thetaRot) + phi1.Py() * std::cos(thetaRot), phi1.Pz(), phi1.M()); + + TLorentzVector pairRot = daughterRot + phi2; + + if (pairRot.Pt() < minExoticPt || pairRot.M() < minExoticMass || pairRot.M() > maxExoticMass) { + continue; + } + + histos.fill(HIST("SEMassPhiPhiRotational"), + pairRot.M(), + pairRot.Pt()); + } + + TLorentzVector pair = phi1 + phi2; + // Mass window range for the phi-phi pair + if (pair.Pt() < minExoticPt || pair.M() < minExoticMass || pair.M() > maxExoticMass) { + continue; + } + histos.fill(HIST("hPhiMass"), phi1.M(), phi2.M(), pair.Pt()); + + // ===================================================================== + // 4-KAON KINEMATIC FIT + // 2 x M(KK) = M(phi) constraints + // ===================================================================== + + TLorentzVector k11, k12, k21, k22; + k11.SetXYZM(Phitrack1.phid1Px(), Phitrack1.phid1Py(), Phitrack1.phid1Pz(), mKPDG); + k12.SetXYZM(Phitrack1.phid2Px(), Phitrack1.phid2Py(), Phitrack1.phid2Pz(), mKPDG); + k21.SetXYZM(Phitrack2.phid1Px(), Phitrack2.phid1Py(), Phitrack2.phid1Pz(), mKPDG); + k22.SetXYZM(Phitrack2.phid2Px(), Phitrack2.phid2Py(), Phitrack2.phid2Pz(), mKPDG); + + FourKFitResult fitResult = fitFourKaons(k11, k12, k21, k22, mKPDG, mPhiPDG, 0.01, 20, 1e-6); + + double refittedMass = pair.M(); + double fitChi2 = -1.0; + double fitProb = -1.0; + + if (fitResult.converged) { + refittedMass = fitResult.refittedMass; + fitChi2 = fitResult.chi2; + fitProb = fitResult.probability; + } + + double deltaM = std::hypot(Phitrack1.phiMass() - mPhiPDG, Phitrack2.phiMass() - mPhiPDG); + + if (fitResult.converged) { + histos.fill(HIST("SEMassPhiPhiRefitted"), + pair.M(), + pair.Pt(), + deltaM, + fitChi2, + fitProb, + phi1.M(), + phi2.M()); + } + } + } + } + PROCESS_SWITCH(doublephimeson, processOpti8, "Process optimised save-event with phi-phi pairs with 4K fit", false); + + double getMomentumCorrection(double phiPt) + { + constexpr int kNBins = 14; + constexpr double kPtEdges[kNBins + 1] = {0.5, 0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 15.0, 20.0}; + constexpr double kEpsilon[kNBins] = {0.0024146539379, 0.00762912304405, 0.00817386742107, 0.00982766792836, 0.0106788436324, 0.0112954252224, 0.0115583787246, 0.0117416086153, 0.012422150094, 0.011968584565, 0.0125561781092, 0.0139372556405, 0.014995530318, 0.0165687496159}; + + if (phiPt < kPtEdges[0] || phiPt >= kPtEdges[kNBins]) + return 0.0; + + auto it = std::upper_bound(std::begin(kPtEdges), std::end(kPtEdges), phiPt); + int binIndex = std::distance(std::begin(kPtEdges), it) - 1; + + return kEpsilon[binIndex]; + } + + TLorentzVector CorrectKaonMomentum(const TLorentzVector& kaon, double epsilon) + { + const double scale = 1.0 + epsilon; + constexpr double mKPDG = o2::constants::physics::MassKPlus; + TLorentzVector corrected; + + // Scale the 3-momentum + TVector3 pCorr = scale * kaon.Vect(); + + // Recalculate energy using fixed kaon mass + corrected.SetVectM(pCorr, mKPDG); + + return corrected; + } + + void processOpti9(aod::RedPhiEvents::iterator const& collision, aod::PhiTracks const& phitracks) + { + if (additionalEvsel && (collision.numPos() < 2 || collision.numNeg() < 2)) { + return; + } + constexpr double mPhiPDG = o2::constants::physics::MassPhi; + constexpr double mKPDG = o2::constants::physics::MassKPlus; + int phimult = 0; + + for (const auto& Phitrack : phitracks) { + const double kpluspt = std::hypot(Phitrack.phid1Px(), Phitrack.phid1Py()); + const double kminuspt = std::hypot(Phitrack.phid2Px(), Phitrack.phid2Py()); + + histos.fill(HIST("hnsigmaTPCTOFKaonBefore"), Phitrack.phid1TPC(), Phitrack.phid1TOF(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonPlusBefore"), Phitrack.phid1TPC(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonMinusBefore"), Phitrack.phid2TPC(), kminuspt); + + // pT cut on kaon tracks + if (kpluspt > maxKaonPt || kminuspt > maxKaonPt) { + continue; + } + + // Mass window selection for phi candidates + if (Phitrack.phiMass() < minPhiMass1 || Phitrack.phiMass() > maxPhiMass1) { + continue; + } + + TLorentzVector phi; + phi.SetXYZM(Phitrack.phiPx(), Phitrack.phiPy(), Phitrack.phiPz(), Phitrack.phiMass()); + + if (phi.Pt() < minPhiPt || phi.Pt() > maxPhiPt) { + continue; + } + + // PID selection for kaon tracks + if (!selectionPID(Phitrack.phid1TPC(), Phitrack.phid1TOF(), Phitrack.phid1TOFHit(), strategyPID1, kpluspt)) { + continue; + } + if (!selectionPID(Phitrack.phid2TPC(), Phitrack.phid2TOF(), Phitrack.phid2TOFHit(), strategyPID2, kminuspt)) { + continue; + } + + histos.fill(HIST("hnsigmaTPCTOFKaon"), Phitrack.phid1TPC(), Phitrack.phid1TOF(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonPlus"), Phitrack.phid1TPC(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonMinus"), Phitrack.phid2TPC(), kminuspt); + + phimult++; + } + + if (phimult < 2) { + return; + } + + for (auto const& Phitrack1 : phitracks) { + // pT cut kaon tracks + const double kplus1pt = std::hypot(Phitrack1.phid1Px(), Phitrack1.phid1Py()); + const double kminus1pt = std::hypot(Phitrack1.phid2Px(), Phitrack1.phid2Py()); + if (kplus1pt > maxKaonPt || kminus1pt > maxKaonPt) { + continue; + } + + // PID selection for kaon tracks + if (!selectionPID(Phitrack1.phid1TPC(), Phitrack1.phid1TOF(), Phitrack1.phid1TOFHit(), strategyPID1, kplus1pt) || + !selectionPID(Phitrack1.phid2TPC(), Phitrack1.phid2TOF(), Phitrack1.phid2TOFHit(), strategyPID2, kminus1pt)) { + continue; + } + + // pT cut on Phi resonance + TLorentzVector phi1; + phi1.SetXYZM(Phitrack1.phiPx(), Phitrack1.phiPy(), Phitrack1.phiPz(), Phitrack1.phiMass()); + if (phi1.Pt() < minPhiPt || phi1.Pt() > maxPhiPt) { + continue; + } + + // Mass window selection for phi candidates + if (Phitrack1.phiMass() < minPhiMass1 || Phitrack1.phiMass() > maxPhiMass1) { + continue; + } + histos.fill(HIST("hPhiMassVsPt"), phi1.M(), phi1.Pt()); + + // Uncorrected phi pT for scale factor lookup + const double phi1UncorrPt = std::hypot(Phitrack1.phiPx(), Phitrack1.phiPy()); + + // Correction based on ORIGINAL phi pT + const double epsilon1 = getMomentumCorrection(phi1UncorrPt); + + // Original kaons + TLorentzVector kplus1, kminus1; + kplus1.SetXYZM(Phitrack1.phid1Px(), Phitrack1.phid1Py(), Phitrack1.phid1Pz(), mKPDG); + kminus1.SetXYZM(Phitrack1.phid2Px(), Phitrack1.phid2Py(), Phitrack1.phid2Pz(), mKPDG); + + // Corrected kaon momenta + TLorentzVector kplus1Corr = CorrectKaonMomentum(kplus1, epsilon1); + TLorentzVector kminus1Corr = CorrectKaonMomentum(kminus1, epsilon1); + TLorentzVector phi1Corr = kplus1Corr + kminus1Corr; + + histos.fill(HIST("hPhiMassVsPtShifted"), phi1Corr.M(), phi1Corr.Pt()); + + for (auto const& Phitrack2 : phitracks) { + + if (Phitrack2.index() <= Phitrack1.index()) { + continue; + } + + // pT cut kaon tracks + const double kplus2pt = std::hypot(Phitrack2.phid1Px(), Phitrack2.phid1Py()); + const double kminus2pt = std::hypot(Phitrack2.phid2Px(), Phitrack2.phid2Py()); + + if (kplus2pt > maxKaonPt || kminus2pt > maxKaonPt) { + continue; + } + + // Mass window selection for phi candidates + if (Phitrack2.phiMass() < minPhiMass1 || Phitrack2.phiMass() > maxPhiMass1) { + continue; + } + + // Uncorrected phi pT for scale factor lookup + const double phi2UncorrPt = std::hypot(Phitrack2.phiPx(), Phitrack2.phiPy()); + + TLorentzVector phi2; + phi2.SetXYZM(Phitrack2.phiPx(), Phitrack2.phiPy(), Phitrack2.phiPz(), Phitrack2.phiMass()); + if (phi2.Pt() < minPhiPt || phi2.Pt() > maxPhiPt) { + continue; + } + + // PID selection for kaon original tracks + if (!selectionPID(Phitrack2.phid1TPC(), Phitrack2.phid1TOF(), Phitrack2.phid1TOFHit(), strategyPID1, kplus2pt) || !selectionPID(Phitrack2.phid2TPC(), Phitrack2.phid2TOF(), Phitrack2.phid2TOFHit(), strategyPID2, kminus2pt)) { + continue; + } + + // Check for shared daughters + if (Phitrack1.phid1Index() == Phitrack2.phid1Index() || + Phitrack1.phid1Index() == Phitrack2.phid2Index() || + Phitrack1.phid2Index() == Phitrack2.phid1Index() || + Phitrack1.phid2Index() == Phitrack2.phid2Index()) { + continue; + } + + // Correction based on ORIGINAL phi pT + const double epsilon2 = getMomentumCorrection(phi2UncorrPt); + + // Original kaons + TLorentzVector kplus2, kminus2; + + kplus2.SetXYZM(Phitrack2.phid1Px(), Phitrack2.phid1Py(), Phitrack2.phid1Pz(), mKPDG); + kminus2.SetXYZM(Phitrack2.phid2Px(), Phitrack2.phid2Py(), Phitrack2.phid2Pz(), mKPDG); + + // Correct kaon momenta + TLorentzVector kplus2Corr = CorrectKaonMomentum(kplus2, epsilon2); + TLorentzVector kminus2Corr = CorrectKaonMomentum(kminus2, epsilon2); + + double kplus2CorrPt = kplus2Corr.Pt(); + double kminus2CorrPt = kminus2Corr.Pt(); + + // Corrected phi + TLorentzVector phi2Corr = kplus2Corr + kminus2Corr; + + // Reconstruct double-phi pair from SHIFTED phi candidates + TLorentzVector pair = phi1 + phi2; + TLorentzVector pairShifted = phi1Corr + phi2Corr; + + if (pair.Pt() < minExoticPt || pair.M() < minExoticMass || pair.M() > maxExoticMass) { + continue; + } + + histos.fill(HIST("hPhiMass"), phi1.M(), phi2.M(), pair.Pt()); + histos.fill(HIST("hPhiMassShifted"), phi1Corr.M(), phi2Corr.M(), pairShifted.Pt()); + + double deltaMShifted = std::hypot(phi1Corr.M() - mPhiPDG, phi2Corr.M() - mPhiPDG); + double deltaM = std::hypot(phi1.M() - mPhiPDG, phi2.M() - mPhiPDG); + + // 4-Kaon Kinematic Fit using SHIFTED kaon vectors + FourKFitResult fitResult = fitFourKaons(kplus1Corr, kminus1Corr, kplus2Corr, kminus2Corr, mKPDG, mPhiPDG, 0.01, 20, 1e-6); + + double refittedMass = pairShifted.M(); + double fitChi2 = -1.0; + double fitProb = -1.0; + + if (fitResult.converged) { + refittedMass = fitResult.refittedMass; + fitChi2 = fitResult.chi2; + fitProb = fitResult.probability; + } + + // Fill shifted double-phi THnSparse + histos.fill(HIST("SEMassPhiPhiShifted"), + pairShifted.M(), + pairShifted.Pt(), + deltaMShifted, + fitChi2, + fitProb, + phi1Corr.M(), + phi2Corr.M()); + } + } + } + PROCESS_SWITCH(doublephimeson, processOpti9, "Process optimised save-event for phi-phi pairs after kaon momentum shift", false); + SliceCache cache; using BinningTypeVertexContributor = ColumnBinningPolicy; From 4c9cb71dafa694f0cef07d6d175f8af9268497e7 Mon Sep 17 00:00:00 2001 From: sawan Date: Tue, 22 Sep 2026 22:44:25 +0530 Subject: [PATCH 4/4] solved clang format errors --- PWGLF/Tasks/Resonances/doublephimeson.cxx | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/PWGLF/Tasks/Resonances/doublephimeson.cxx b/PWGLF/Tasks/Resonances/doublephimeson.cxx index beb670129ab..3850a5cf2f9 100644 --- a/PWGLF/Tasks/Resonances/doublephimeson.cxx +++ b/PWGLF/Tasks/Resonances/doublephimeson.cxx @@ -27,9 +27,6 @@ #include #include -#include "TLorentzVector.h" -#include "TMath.h" -#include "TMatrixD.h" #include #include #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) @@ -37,6 +34,7 @@ #include #include #include +#include #include #include