diff --git a/PWGCF/Femto/Core/tripletHistManager.h b/PWGCF/Femto/Core/tripletHistManager.h index ca1e0779369..8d5600c3331 100644 --- a/PWGCF/Femto/Core/tripletHistManager.h +++ b/PWGCF/Femto/Core/tripletHistManager.h @@ -28,6 +28,7 @@ #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include +#include #include #include #include @@ -44,6 +45,10 @@ enum TripletHist { // standard 1D kQ3, kMt, // averate mt of all pairs in the triplet + // kstar betweeen single particles + kKstar12, + kKstar13, + kKstar23, // standard 2D kPt1VsQ3, kPt2VsQ3, @@ -106,6 +111,7 @@ struct ConfTripletBinning : o2::framework::ConfigurableGroup { std::string prefix = std::string("TripletBinning"); o2::framework::Configurable plot1D{"plot1D", true, "Enable 1D histograms"}; o2::framework::Configurable plot2D{"plot2D", true, "Enable 2D histograms"}; + o2::framework::Configurable plotKstar{"plotKstar", false, "Enable kstar histograms"}; o2::framework::Configurable plotPt1VsPt2VsPt3{"plotPt1VsPt2VsPt3", false, "Enable 3D histogram (Pt1 vs Pt2 vs Pt3)"}; o2::framework::Configurable plotQ3VsPt1VsPt2VsPt3{"plotQ3VsPt1VsPt2VsPt3", false, "Enable 4D histogram (Q3 vs Pt1 vs Pt2 vs Pt3)"}; o2::framework::Configurable plotQ3VsMtVsMult{"plotQ3VsMtVsMult", false, "Enable 3D histogram (Q3 vs Mt vs Mult)"}; @@ -113,6 +119,7 @@ struct ConfTripletBinning : o2::framework::ConfigurableGroup { o2::framework::Configurable plotQ3VsMtVsPt1VsPt2VsPt3VsMult{"plotQ3VsMtVsPt1VsPt2VsPt3VsMult", false, "Enable 6D histogram (Q3 vs Mt Vs Pt1 vs Pt2 vs Pt3 vs Mult)"}; o2::framework::Configurable plotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent{"plotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent", false, "Enable 7D histogram (Q3 vs Mt Vs Pt1 vs Pt2 vs Pt3 vs Mult vs Cent)"}; o2::framework::ConfigurableAxis q3{"q3", {{600, 0, 6}}, "q3"}; + o2::framework::ConfigurableAxis kstar{"kstar", {{600, 0, 6}}, "kstar (between pairs of particles)"}; o2::framework::ConfigurableAxis mt{"mt", {{500, 0.8, 5.8}}, "mt"}; o2::framework::ConfigurableAxis multiplicity{"multiplicity", {{50, 0, 200}}, "multiplicity"}; o2::framework::ConfigurableAxis centrality{"centrality", {{10, 0, 100}}, "centrality (mult. percentile)"}; @@ -142,6 +149,9 @@ constexpr std::array, kTripletHistogramLast> // 1D {kQ3, o2::framework::HistType::kTH1F, "hQ3", "Q_{3}; Q_{3} (GeV/#it{c}); Entries"}, {kMt, o2::framework::HistType::kTH1F, "hMt", "transverse mass; m_{T} (GeV/#it{c}^{2}); Entries"}, + {kKstar12, o2::framework::HistType::kTH1F, "hKstar12", "k* between particle 1 and particle 2; k* (GeV/#it{c}); Entries"}, + {kKstar13, o2::framework::HistType::kTH1F, "hKstar13", "k* between particle 1 and particle 3; k* (GeV/#it{c}); Entries"}, + {kKstar23, o2::framework::HistType::kTH1F, "hKstar23", "k* between particle 2 and particle 3; k* (GeV/#it{c}); Entries"}, // 2D {kPt1VsQ3, o2::framework::HistType::kTH2F, "hPt1VsQ3", "p_{T,1} vs Q_{3}; p_{T,1} (GeV/#it{c}); Q_{3} (GeV/#it{c})"}, {kPt2VsQ3, o2::framework::HistType::kTH2F, "hPt2VsQ3", "p_{T,2} vs Q_{3}; p_{T,2} (GeV/#it{c}); Q_{3} (GeV/#it{c})"}, @@ -173,6 +183,9 @@ constexpr std::array, kTripletHistogramLast> #define TRIPLET_HIST_ANALYSIS_MAP(conf, confMixing) \ {kQ3, {(conf).q3}}, \ {kMt, {(conf).mt}}, \ + {kKstar12, {(conf).kstar}}, \ + {kKstar13, {(conf).kstar}}, \ + {kKstar23, {(conf).kstar}}, \ {kPt1VsQ3, {(conf).pt1, (conf).q3}}, \ {kPt2VsQ3, {(conf).pt2, (conf).q3}}, \ {kPt3VsQ3, {(conf).pt3, (conf).q3}}, \ @@ -261,6 +274,7 @@ class TripletHistManager mMtMin = ConfTripletCuts.mtMin.value; mMtMax = ConfTripletCuts.mtMax.value; + mPlotKstar = ConfTripletBinning.plotKstar.value; mPlotPt1VsPt2VsPt3 = ConfTripletBinning.plotPt1VsPt2VsPt3.value; mPlotQ3VsPt1VsPt2VsPt3 = ConfTripletBinning.plotQ3VsPt1VsPt2VsPt3.value; mPlotQ3VsMtVsMult = ConfTripletBinning.plotQ3VsMtVsMult.value; @@ -314,6 +328,10 @@ class TripletHistManager // set Q3 mQ3 = getQ3(mParticle1, mParticle2, mParticle3); + mKstar12 = getKstar(mParticle1, mParticle2); + mKstar13 = getKstar(mParticle1, mParticle3); + mKstar23 = getKstar(mParticle2, mParticle3); + // if one of the particles has a mass getter, we cache the value for the filling later if constexpr (utils::HasMass) { mMass1 = particle1.mass(); @@ -495,6 +513,24 @@ class TripletHistManager return static_cast(std::sqrt(-q)); } + float getKstar(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) + { + // Use Cartesian 4-vectors: addition/M2() become pure arithmetic + const ROOT::Math::PxPyPzEVector p1(part1); + const ROOT::Math::PxPyPzEVector p2(part2); + + // Mandelstam s = (p1 + p2)^2 + const double s = (p1 + p2).M2(); + const double m1sq = p1.M2(); + const double m2sq = p2.M2(); + + // Källen function λ(s, m1^2, m2^2) = (s - m1^2 - m2^2)² - 4*m1^2*m2^2 + const double kallen = (s - m1sq - m2sq) * (s - m1sq - m2sq) - 4.0 * m1sq * m2sq; + + // k* = 0.5 * sqrt(λ/s) + return static_cast(0.5 * std::sqrt(std::max(0.0, kallen) / s)); + } + void initAnalysis(std::map> const& Specs) { std::string analysisDir = std::string(prefix) + std::string(AnalysisDir); @@ -512,6 +548,12 @@ class TripletHistManager // TODO: implement histograms differential im mass of the particles } + if (mPlotKstar) { + mHistogramRegistry->add(analysisDir + getHistNameV2(kKstar12, HistTable), getHistDesc(kKstar12, HistTable), getHistType(kKstar12, HistTable), {Specs.at(kKstar12)}); + mHistogramRegistry->add(analysisDir + getHistNameV2(kKstar13, HistTable), getHistDesc(kKstar13, HistTable), getHistType(kKstar13, HistTable), {Specs.at(kKstar13)}); + mHistogramRegistry->add(analysisDir + getHistNameV2(kKstar23, HistTable), getHistDesc(kKstar23, HistTable), getHistType(kKstar23, HistTable), {Specs.at(kKstar23)}); + } + if (mPlotPt1VsPt2VsPt3) { mHistogramRegistry->add(analysisDir + getHistNameV2(kPt1VsPt2VsPt3, HistTable), getHistDesc(kPt1VsPt2VsPt3, HistTable), getHistType(kPt1VsPt2VsPt3, HistTable), {Specs.at(kPt1VsPt2VsPt3)}); } @@ -577,6 +619,12 @@ class TripletHistManager mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3VsCent, HistTable)), mQ3, mCent); } + if (mPlotKstar) { + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kKstar12, HistTable)), mKstar12); + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kKstar13, HistTable)), mKstar13); + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kKstar23, HistTable)), mKstar23); + } + if (mPlotPt1VsPt2VsPt3) { mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kPt1VsPt2VsPt3, HistTable)), mParticle1.Pt(), mParticle2.Pt(), mParticle3.Pt()); } @@ -658,6 +706,9 @@ class TripletHistManager float mMt = 0.f; float mMult = 0.f; float mCent = 0.f; + float mKstar12 = 0.f; + float mKstar13 = 0.f; + float mKstar23 = 0.f; // mc ROOT::Math::PtEtaPhiMVector mTrueParticle1; @@ -680,6 +731,7 @@ class TripletHistManager bool mPlot1d = true; bool mPlot2d = true; + bool mPlotKstar = false; bool mPlotPt1VsPt2VsPt3 = false; bool mPlotQ3VsPt1VsPt2VsPt3 = false; bool mPlotQ3VsMtVsMult = false;