2828#include < Math/Vector4D.h> // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h)
2929#include < Math/Vector4Dfwd.h>
3030
31+ #include < algorithm>
3132#include < array>
3233#include < cmath>
3334#include < cstdint>
@@ -44,6 +45,10 @@ enum TripletHist {
4445 // standard 1D
4546 kQ3 ,
4647 kMt , // averate mt of all pairs in the triplet
48+ // kstar betweeen single particles
49+ kKstar12 ,
50+ kKstar13 ,
51+ kKstar23 ,
4752 // standard 2D
4853 kPt1VsQ3 ,
4954 kPt2VsQ3 ,
@@ -106,13 +111,15 @@ struct ConfTripletBinning : o2::framework::ConfigurableGroup {
106111 std::string prefix = std::string(" TripletBinning" );
107112 o2::framework::Configurable<bool > plot1D{" plot1D" , true , " Enable 1D histograms" };
108113 o2::framework::Configurable<bool > plot2D{" plot2D" , true , " Enable 2D histograms" };
114+ o2::framework::Configurable<bool > plotKstar{" plotKstar" , false , " Enable kstar histograms" };
109115 o2::framework::Configurable<bool > plotPt1VsPt2VsPt3{" plotPt1VsPt2VsPt3" , false , " Enable 3D histogram (Pt1 vs Pt2 vs Pt3)" };
110116 o2::framework::Configurable<bool > plotQ3VsPt1VsPt2VsPt3{" plotQ3VsPt1VsPt2VsPt3" , false , " Enable 4D histogram (Q3 vs Pt1 vs Pt2 vs Pt3)" };
111117 o2::framework::Configurable<bool > plotQ3VsMtVsMult{" plotQ3VsMtVsMult" , false , " Enable 3D histogram (Q3 vs Mt vs Mult)" };
112118 o2::framework::Configurable<bool > plotQ3VsMtVsMultVsCent{" plotQ3VsMtVsMultVsCent" , false , " Enable 4D histogram (Q3 vs Mt vs Mult Vs Cent)" };
113119 o2::framework::Configurable<bool > plotQ3VsMtVsPt1VsPt2VsPt3VsMult{" plotQ3VsMtVsPt1VsPt2VsPt3VsMult" , false , " Enable 6D histogram (Q3 vs Mt Vs Pt1 vs Pt2 vs Pt3 vs Mult)" };
114120 o2::framework::Configurable<bool > plotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent{" plotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent" , false , " Enable 7D histogram (Q3 vs Mt Vs Pt1 vs Pt2 vs Pt3 vs Mult vs Cent)" };
115121 o2::framework::ConfigurableAxis q3{" q3" , {{600 , 0 , 6 }}, " q3" };
122+ o2::framework::ConfigurableAxis kstar{" kstar" , {{600 , 0 , 6 }}, " kstar (between pairs of particles)" };
116123 o2::framework::ConfigurableAxis mt{" mt" , {{500 , 0.8 , 5.8 }}, " mt" };
117124 o2::framework::ConfigurableAxis multiplicity{" multiplicity" , {{50 , 0 , 200 }}, " multiplicity" };
118125 o2::framework::ConfigurableAxis centrality{" centrality" , {{10 , 0 , 100 }}, " centrality (mult. percentile)" };
@@ -142,6 +149,9 @@ constexpr std::array<histmanager::HistInfo<TripletHist>, kTripletHistogramLast>
142149 // 1D
143150 {kQ3 , o2::framework::HistType::kTH1F , " hQ3" , " Q_{3}; Q_{3} (GeV/#it{c}); Entries" },
144151 {kMt , o2::framework::HistType::kTH1F , " hMt" , " transverse mass; m_{T} (GeV/#it{c}^{2}); Entries" },
152+ {kKstar12 , o2::framework::HistType::kTH1F , " hKstar12" , " k* between particle 1 and particle 2; k* (GeV/#it{c}); Entries" },
153+ {kKstar13 , o2::framework::HistType::kTH1F , " hKstar13" , " k* between particle 1 and particle 3; k* (GeV/#it{c}); Entries" },
154+ {kKstar23 , o2::framework::HistType::kTH1F , " hKstar23" , " k* between particle 2 and particle 3; k* (GeV/#it{c}); Entries" },
145155 // 2D
146156 {kPt1VsQ3 , o2::framework::HistType::kTH2F , " hPt1VsQ3" , " p_{T,1} vs Q_{3}; p_{T,1} (GeV/#it{c}); Q_{3} (GeV/#it{c})" },
147157 {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<histmanager::HistInfo<TripletHist>, kTripletHistogramLast>
173183#define TRIPLET_HIST_ANALYSIS_MAP (conf, confMixing ) \
174184 {kQ3 , {(conf).q3 }}, \
175185 {kMt , {(conf).mt }}, \
186+ {kKstar12 , {(conf).kstar }}, \
187+ {kKstar13 , {(conf).kstar }}, \
188+ {kKstar23 , {(conf).kstar }}, \
176189 {kPt1VsQ3 , {(conf).pt1 , (conf).q3 }}, \
177190 {kPt2VsQ3 , {(conf).pt2 , (conf).q3 }}, \
178191 {kPt3VsQ3 , {(conf).pt3 , (conf).q3 }}, \
@@ -261,6 +274,7 @@ class TripletHistManager
261274 mMtMin = ConfTripletCuts.mtMin .value ;
262275 mMtMax = ConfTripletCuts.mtMax .value ;
263276
277+ mPlotKstar = ConfTripletBinning.plotKstar .value ;
264278 mPlotPt1VsPt2VsPt3 = ConfTripletBinning.plotPt1VsPt2VsPt3 .value ;
265279 mPlotQ3VsPt1VsPt2VsPt3 = ConfTripletBinning.plotQ3VsPt1VsPt2VsPt3 .value ;
266280 mPlotQ3VsMtVsMult = ConfTripletBinning.plotQ3VsMtVsMult .value ;
@@ -314,6 +328,10 @@ class TripletHistManager
314328 // set Q3
315329 mQ3 = getQ3 (mParticle1 , mParticle2 , mParticle3 );
316330
331+ mKstar12 = getKstar (mParticle1 , mParticle2 );
332+ mKstar13 = getKstar (mParticle1 , mParticle3 );
333+ mKstar23 = getKstar (mParticle2 , mParticle3 );
334+
317335 // if one of the particles has a mass getter, we cache the value for the filling later
318336 if constexpr (utils::HasMass<T1 >) {
319337 mMass1 = particle1.mass ();
@@ -495,6 +513,24 @@ class TripletHistManager
495513 return static_cast <float >(std::sqrt (-q));
496514 }
497515
516+ float getKstar (ROOT ::Math::PtEtaPhiMVector const & part1, ROOT ::Math::PtEtaPhiMVector const & part2)
517+ {
518+ // Use Cartesian 4-vectors: addition/M2() become pure arithmetic
519+ const ROOT ::Math::PxPyPzEVector p1 (part1);
520+ const ROOT ::Math::PxPyPzEVector p2 (part2);
521+
522+ // Mandelstam s = (p1 + p2)^2
523+ const double s = (p1 + p2).M2 ();
524+ const double m1sq = p1.M2 ();
525+ const double m2sq = p2.M2 ();
526+
527+ // Källen function λ(s, m1^2, m2^2) = (s - m1^2 - m2^2)² - 4*m1^2*m2^2
528+ const double kallen = (s - m1sq - m2sq) * (s - m1sq - m2sq) - 4.0 * m1sq * m2sq;
529+
530+ // k* = 0.5 * sqrt(λ/s)
531+ return static_cast <float >(0.5 * std::sqrt (std::max (0.0 , kallen) / s));
532+ }
533+
498534 void initAnalysis (std::map<TripletHist, std::vector<o2::framework::AxisSpec>> const & Specs)
499535 {
500536 std::string analysisDir = std::string (prefix) + std::string (AnalysisDir);
@@ -512,6 +548,12 @@ class TripletHistManager
512548 // TODO: implement histograms differential im mass of the particles
513549 }
514550
551+ if (mPlotKstar ) {
552+ mHistogramRegistry ->add (analysisDir + getHistNameV2 (kKstar12 , HistTable), getHistDesc (kKstar12 , HistTable), getHistType (kKstar12 , HistTable), {Specs.at (kKstar12 )});
553+ mHistogramRegistry ->add (analysisDir + getHistNameV2 (kKstar13 , HistTable), getHistDesc (kKstar13 , HistTable), getHistType (kKstar13 , HistTable), {Specs.at (kKstar13 )});
554+ mHistogramRegistry ->add (analysisDir + getHistNameV2 (kKstar23 , HistTable), getHistDesc (kKstar23 , HistTable), getHistType (kKstar23 , HistTable), {Specs.at (kKstar23 )});
555+ }
556+
515557 if (mPlotPt1VsPt2VsPt3 ) {
516558 mHistogramRegistry ->add (analysisDir + getHistNameV2 (kPt1VsPt2VsPt3 , HistTable), getHistDesc (kPt1VsPt2VsPt3 , HistTable), getHistType (kPt1VsPt2VsPt3 , HistTable), {Specs.at (kPt1VsPt2VsPt3 )});
517559 }
@@ -577,6 +619,12 @@ class TripletHistManager
577619 mHistogramRegistry ->fill (HIST (prefix) + HIST (AnalysisDir) + HIST (getHistName (kQ3VsCent , HistTable)), mQ3 , mCent );
578620 }
579621
622+ if (mPlotKstar ) {
623+ mHistogramRegistry ->fill (HIST (prefix) + HIST (AnalysisDir) + HIST (getHistName (kKstar12 , HistTable)), mKstar12 );
624+ mHistogramRegistry ->fill (HIST (prefix) + HIST (AnalysisDir) + HIST (getHistName (kKstar13 , HistTable)), mKstar13 );
625+ mHistogramRegistry ->fill (HIST (prefix) + HIST (AnalysisDir) + HIST (getHistName (kKstar23 , HistTable)), mKstar23 );
626+ }
627+
580628 if (mPlotPt1VsPt2VsPt3 ) {
581629 mHistogramRegistry ->fill (HIST (prefix) + HIST (AnalysisDir) + HIST (getHistName (kPt1VsPt2VsPt3 , HistTable)), mParticle1 .Pt (), mParticle2 .Pt (), mParticle3 .Pt ());
582630 }
@@ -658,6 +706,9 @@ class TripletHistManager
658706 float mMt = 0 .f;
659707 float mMult = 0 .f;
660708 float mCent = 0 .f;
709+ float mKstar12 = 0 .f;
710+ float mKstar13 = 0 .f;
711+ float mKstar23 = 0 .f;
661712
662713 // mc
663714 ROOT ::Math::PtEtaPhiMVector mTrueParticle1 ;
@@ -680,6 +731,7 @@ class TripletHistManager
680731 bool mPlot1d = true ;
681732 bool mPlot2d = true ;
682733
734+ bool mPlotKstar = false ;
683735 bool mPlotPt1VsPt2VsPt3 = false ;
684736 bool mPlotQ3VsPt1VsPt2VsPt3 = false ;
685737 bool mPlotQ3VsMtVsMult = false ;
0 commit comments