Skip to content

Commit 685cf1f

Browse files
IgorPtakIgorPtak
authored
[PWGCF] Add Bertsch-Pratt LCMS decomposition for femto pairs (#16380)
1 parent d5c009f commit 685cf1f

1 file changed

Lines changed: 80 additions & 2 deletions

File tree

‎PWGCF/Femto/Core/pairHistManager.h‎

Lines changed: 80 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,7 @@
3939
#include <map>
4040
#include <string>
4141
#include <string_view>
42+
#include <tuple>
4243
#include <unordered_set>
4344
#include <vector>
4445

@@ -111,7 +112,11 @@ enum PairHist {
111112

112113
// angular
113114
kDeltaEtaDeltaPhi,
114-
115+
// Bertsch-Pratt 3D decomposition in LCMS
116+
kQout,
117+
kQside,
118+
kQlong,
119+
kQoutQsideQlong,
115120
kPairHistogramLast
116121
};
117122

@@ -174,6 +179,10 @@ struct ConfPairBinning : o2::framework::ConfigurableGroup {
174179
o2::framework::Configurable<int> transverseMassType{"transverseMassType", static_cast<int>(modes::TransverseMassType::kAveragePdgMass), "Type of transverse mass (0-> Average Pdg Mass, 1-> Reduced Pdg Mass, 2-> Mt from combined 4 vector)"};
175180
o2::framework::ConfigurableAxis binningDeltaEta{"binningDeltaEta", {{35, -1.6, 1.6}}, "Delta eta"};
176181
o2::framework::ConfigurableAxis binningDeltaPhi{"binningDeltaPhi", {{35, -o2::constants::math::PIHalf, 3 * o2::constants::math::PIHalf}}, "Delta phi"};
182+
o2::framework::Configurable<bool> plotBertschPratt{"plotBertschPratt", false, "Enable 1D projections and 3D (q_out, q_side, q_long) Bertsch-Pratt histograms in LCMS"};
183+
o2::framework::ConfigurableAxis qout{"qout", {{300, -1.5f, 1.5f}}, "q_{out} (GeV/c) in LCMS"};
184+
o2::framework::ConfigurableAxis qside{"qside", {{300, -1.5f, 1.5f}}, "q_{side} (GeV/c) in LCMS"};
185+
o2::framework::ConfigurableAxis qlong{"qlong", {{300, -1.5f, 1.5f}}, "q_{long} (GeV/c) in LCMS"};
177186
};
178187

179188
struct ConfPairCuts : o2::framework::ConfigurableGroup {
@@ -250,6 +259,11 @@ constexpr std::array<histmanager::HistInfo<PairHist>, kPairHistogramLast>
250259
{kMeVtz1VsMult1VsCent1VsVtz2VsMult2VsCent2, o2::framework::HistType::kTHnSparseF, "hVtz1VsMult1VsCent1VsVtz2VsMult2VsCent2", "Mixing bins; V_{z,1} (cm); multiplicity_{1}; centrality_{1} (%); V_{z,2} (cm); multiplicity_{2}; centrality_{2} (%)"},
251260
// angular
252261
{kDeltaEtaDeltaPhi, o2::framework::HistType::kTH2F, "hDeltaEtaDeltaPhi", "#Delta#phi vs #Delta#eta; #Delta#phi; #Delta#eta"},
262+
// Bertsch-Pratt 3D decomposition in LCMS
263+
{kQout, o2::framework::HistType::kTH1F, "hQout", "q_{out} in LCMS; q_{out} (GeV/#it{c}); Entries"},
264+
{kQside, o2::framework::HistType::kTH1F, "hQside", "q_{side} in LCMS; q_{side} (GeV/#it{c}); Entries"},
265+
{kQlong, o2::framework::HistType::kTH1F, "hQlong", "q_{long} in LCMS; q_{long} (GeV/#it{c}); Entries"},
266+
{kQoutQsideQlong, o2::framework::HistType::kTH3F, "hQoutQsideQlong", "Bertsch-Pratt 3D; q_{out} (GeV/#it{c}); q_{side} (GeV/#it{c}); q_{long} (GeV/#it{c})"},
253267
}};
254268

255269
#define PAIR_HIST_ANALYSIS_MAP(confAnalysis, confMixing) \
@@ -294,7 +308,11 @@ constexpr std::array<histmanager::HistInfo<PairHist>, kPairHistogramLast>
294308
{kMeMixingWindowRaw, {confMixing.particleBinning}}, \
295309
{kMeMixingWindowEffective, {confMixing.particleBinning}}, \
296310
{kMeNpart1VsNpart2, {confMixing.particleBinning, confMixing.particleBinning}}, \
297-
{kMeVtz1VsMult1VsCent1VsVtz2VsMult2VsCent2, {confMixing.vtxBins, confMixing.multBins, confMixing.centBins, confMixing.vtxBins, confMixing.multBins, confMixing.centBins}},
311+
{kMeVtz1VsMult1VsCent1VsVtz2VsMult2VsCent2, {confMixing.vtxBins, confMixing.multBins, confMixing.centBins, confMixing.vtxBins, confMixing.multBins, confMixing.centBins}}, \
312+
{kQout, {confAnalysis.qout}}, \
313+
{kQside, {confAnalysis.qside}}, \
314+
{kQlong, {confAnalysis.qlong}}, \
315+
{kQoutQsideQlong, {confAnalysis.qout, confAnalysis.qside, confAnalysis.qlong}},
298316

299317
#define PAIR_HIST_MC_MAP(conf) \
300318
{kTrueKstarVsKstar, {conf.kstar, conf.kstar}}, \
@@ -388,6 +406,7 @@ class PairHistManager
388406

389407
mPlotDalitz = ConfPairBinning.plotDalitz.value;
390408
mPlotDeltaEtaDeltaPhi = ConfPairBinning.plotDeltaEtaDeltaPhi.value;
409+
mPlotBertschPratt = ConfPairBinning.plotBertschPratt.value;
391410

392411
// transverse mass type
393412
mMtType = static_cast<modes::TransverseMassType>(ConfPairBinning.transverseMassType.value);
@@ -488,6 +507,10 @@ class PairHistManager
488507
// set kstar
489508
mKstar = getKstar(mParticle1, mParticle2);
490509

510+
if (mPlotBertschPratt) {
511+
std::tie(mQout, mQside, mQlong) = computeBertschPrattLCMS(mParticle1, mParticle2);
512+
}
513+
491514
if (mPlotDeltaEtaDeltaPhi) {
492515
mDeltaEta = particle1.eta() - particle2.eta();
493516
mDeltaPhi = RecoDecay::constrainAngle(particle1.phi() - particle2.phi(), -o2::constants::math::PIHalf);
@@ -735,6 +758,12 @@ class PairHistManager
735758
if (mPlotDeltaEtaDeltaPhi) {
736759
mHistogramRegistry->add(analysisDir + getHistNameV2(kDeltaEtaDeltaPhi, HistTable), getHistDesc(kDeltaEtaDeltaPhi, HistTable), getHistType(kDeltaEtaDeltaPhi, HistTable), {Specs.at(kDeltaEtaDeltaPhi)});
737760
}
761+
if (mPlotBertschPratt) {
762+
mHistogramRegistry->add(analysisDir + getHistNameV2(kQout, HistTable), getHistDesc(kQout, HistTable), getHistType(kQout, HistTable), {Specs.at(kQout)});
763+
mHistogramRegistry->add(analysisDir + getHistNameV2(kQside, HistTable), getHistDesc(kQside, HistTable), getHistType(kQside, HistTable), {Specs.at(kQside)});
764+
mHistogramRegistry->add(analysisDir + getHistNameV2(kQlong, HistTable), getHistDesc(kQlong, HistTable), getHistType(kQlong, HistTable), {Specs.at(kQlong)});
765+
mHistogramRegistry->add(analysisDir + getHistNameV2(kQoutQsideQlong, HistTable), getHistDesc(kQoutQsideQlong, HistTable), getHistType(kQoutQsideQlong, HistTable), {Specs.at(kQoutQsideQlong)});
766+
}
738767
}
739768

740769
void initMc(std::map<PairHist, std::vector<o2::framework::AxisSpec>> const& Specs)
@@ -850,6 +879,12 @@ class PairHistManager
850879
if (mPlotDeltaEtaDeltaPhi) {
851880
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kDeltaEtaDeltaPhi, HistTable)), mDeltaPhi, mDeltaEta);
852881
}
882+
if (mPlotBertschPratt) {
883+
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQout, HistTable)), mQout);
884+
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQside, HistTable)), mQside);
885+
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQlong, HistTable)), mQlong);
886+
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQoutQsideQlong, HistTable)), mQout, mQside, mQlong);
887+
}
853888
}
854889

855890
void fillMc()
@@ -920,6 +955,43 @@ class PairHistManager
920955
return static_cast<float>(0.5 * std::sqrt(std::max(0.0, kallen) / s));
921956
}
922957

958+
std::tuple<float, float, float> computeBertschPrattLCMS(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2)
959+
{
960+
const ROOT::Math::PxPyPzEVector p1(part1);
961+
const ROOT::Math::PxPyPzEVector p2(part2);
962+
const ROOT::Math::PxPyPzEVector pSum = p1 + p2;
963+
964+
const double tPx = pSum.Px();
965+
const double tPy = pSum.Py();
966+
const double tPz = pSum.Pz();
967+
const double tE = pSum.E();
968+
969+
const double tPt = std::sqrt(tPx * tPx + tPy * tPy);
970+
const double tMt = std::sqrt(tE * tE - tPz * tPz);
971+
972+
static constexpr double kMinTransverseMomentum = 1e-9;
973+
if (tPt < kMinTransverseMomentum || tMt < kMinTransverseMomentum) {
974+
return {0.0, 0.0, 0.0};
975+
}
976+
977+
const double betaL = tPz / tE;
978+
const double gammaL = tE / tMt;
979+
980+
const double kout1 = (p1.Px() * tPx + p1.Py() * tPy) / tPt;
981+
const double kside1 = (-p1.Px() * tPy + p1.Py() * tPx) / tPt;
982+
const double klong1 = gammaL * (p1.Pz() - betaL * p1.E());
983+
984+
const double kout2 = (p2.Px() * tPx + p2.Py() * tPy) / tPt;
985+
const double kside2 = (p2.Py() * tPx - p2.Px() * tPy) / tPt;
986+
const double klong2 = gammaL * (p2.Pz() - betaL * p2.E());
987+
988+
float qOut = static_cast<float>(kout1 - kout2);
989+
float qSide = static_cast<float>(kside1 - kside2);
990+
float qLong = static_cast<float>(klong1 - klong2);
991+
992+
return {qOut, qSide, qLong};
993+
}
994+
923995
o2::framework::HistogramRegistry* mHistogramRegistry = nullptr;
924996
bool mUsePdgMass = true;
925997
double mPdgMass1 = 0.;
@@ -995,6 +1067,12 @@ class PairHistManager
9951067
bool mPlotDalitz = false;
9961068
bool mPlotDeltaEtaDeltaPhi = false;
9971069

1070+
bool mPlotBertschPratt = false;
1071+
1072+
float mQout = 0.f;
1073+
float mQside = 0.f;
1074+
float mQlong = 0.f;
1075+
9981076
float mDeltaEta = 0.f;
9991077
float mDeltaPhi = 0.f;
10001078

0 commit comments

Comments
 (0)