|
39 | 39 | #include <map> |
40 | 40 | #include <string> |
41 | 41 | #include <string_view> |
| 42 | +#include <tuple> |
42 | 43 | #include <unordered_set> |
43 | 44 | #include <vector> |
44 | 45 |
|
@@ -111,7 +112,11 @@ enum PairHist { |
111 | 112 |
|
112 | 113 | // angular |
113 | 114 | kDeltaEtaDeltaPhi, |
114 | | - |
| 115 | + // Bertsch-Pratt 3D decomposition in LCMS |
| 116 | + kQout, |
| 117 | + kQside, |
| 118 | + kQlong, |
| 119 | + kQoutQsideQlong, |
115 | 120 | kPairHistogramLast |
116 | 121 | }; |
117 | 122 |
|
@@ -174,6 +179,10 @@ struct ConfPairBinning : o2::framework::ConfigurableGroup { |
174 | 179 | 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)"}; |
175 | 180 | o2::framework::ConfigurableAxis binningDeltaEta{"binningDeltaEta", {{35, -1.6, 1.6}}, "Delta eta"}; |
176 | 181 | 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"}; |
177 | 186 | }; |
178 | 187 |
|
179 | 188 | struct ConfPairCuts : o2::framework::ConfigurableGroup { |
@@ -250,6 +259,11 @@ constexpr std::array<histmanager::HistInfo<PairHist>, kPairHistogramLast> |
250 | 259 | {kMeVtz1VsMult1VsCent1VsVtz2VsMult2VsCent2, o2::framework::HistType::kTHnSparseF, "hVtz1VsMult1VsCent1VsVtz2VsMult2VsCent2", "Mixing bins; V_{z,1} (cm); multiplicity_{1}; centrality_{1} (%); V_{z,2} (cm); multiplicity_{2}; centrality_{2} (%)"}, |
251 | 260 | // angular |
252 | 261 | {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})"}, |
253 | 267 | }}; |
254 | 268 |
|
255 | 269 | #define PAIR_HIST_ANALYSIS_MAP(confAnalysis, confMixing) \ |
@@ -294,7 +308,11 @@ constexpr std::array<histmanager::HistInfo<PairHist>, kPairHistogramLast> |
294 | 308 | {kMeMixingWindowRaw, {confMixing.particleBinning}}, \ |
295 | 309 | {kMeMixingWindowEffective, {confMixing.particleBinning}}, \ |
296 | 310 | {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}}, |
298 | 316 |
|
299 | 317 | #define PAIR_HIST_MC_MAP(conf) \ |
300 | 318 | {kTrueKstarVsKstar, {conf.kstar, conf.kstar}}, \ |
@@ -388,6 +406,7 @@ class PairHistManager |
388 | 406 |
|
389 | 407 | mPlotDalitz = ConfPairBinning.plotDalitz.value; |
390 | 408 | mPlotDeltaEtaDeltaPhi = ConfPairBinning.plotDeltaEtaDeltaPhi.value; |
| 409 | + mPlotBertschPratt = ConfPairBinning.plotBertschPratt.value; |
391 | 410 |
|
392 | 411 | // transverse mass type |
393 | 412 | mMtType = static_cast<modes::TransverseMassType>(ConfPairBinning.transverseMassType.value); |
@@ -488,6 +507,10 @@ class PairHistManager |
488 | 507 | // set kstar |
489 | 508 | mKstar = getKstar(mParticle1, mParticle2); |
490 | 509 |
|
| 510 | + if (mPlotBertschPratt) { |
| 511 | + std::tie(mQout, mQside, mQlong) = computeBertschPrattLCMS(mParticle1, mParticle2); |
| 512 | + } |
| 513 | + |
491 | 514 | if (mPlotDeltaEtaDeltaPhi) { |
492 | 515 | mDeltaEta = particle1.eta() - particle2.eta(); |
493 | 516 | mDeltaPhi = RecoDecay::constrainAngle(particle1.phi() - particle2.phi(), -o2::constants::math::PIHalf); |
@@ -735,6 +758,12 @@ class PairHistManager |
735 | 758 | if (mPlotDeltaEtaDeltaPhi) { |
736 | 759 | mHistogramRegistry->add(analysisDir + getHistNameV2(kDeltaEtaDeltaPhi, HistTable), getHistDesc(kDeltaEtaDeltaPhi, HistTable), getHistType(kDeltaEtaDeltaPhi, HistTable), {Specs.at(kDeltaEtaDeltaPhi)}); |
737 | 760 | } |
| 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 | + } |
738 | 767 | } |
739 | 768 |
|
740 | 769 | void initMc(std::map<PairHist, std::vector<o2::framework::AxisSpec>> const& Specs) |
@@ -850,6 +879,12 @@ class PairHistManager |
850 | 879 | if (mPlotDeltaEtaDeltaPhi) { |
851 | 880 | mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kDeltaEtaDeltaPhi, HistTable)), mDeltaPhi, mDeltaEta); |
852 | 881 | } |
| 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 | + } |
853 | 888 | } |
854 | 889 |
|
855 | 890 | void fillMc() |
@@ -920,6 +955,43 @@ class PairHistManager |
920 | 955 | return static_cast<float>(0.5 * std::sqrt(std::max(0.0, kallen) / s)); |
921 | 956 | } |
922 | 957 |
|
| 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 | + |
923 | 995 | o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; |
924 | 996 | bool mUsePdgMass = true; |
925 | 997 | double mPdgMass1 = 0.; |
@@ -995,6 +1067,12 @@ class PairHistManager |
995 | 1067 | bool mPlotDalitz = false; |
996 | 1068 | bool mPlotDeltaEtaDeltaPhi = false; |
997 | 1069 |
|
| 1070 | + bool mPlotBertschPratt = false; |
| 1071 | + |
| 1072 | + float mQout = 0.f; |
| 1073 | + float mQside = 0.f; |
| 1074 | + float mQlong = 0.f; |
| 1075 | + |
998 | 1076 | float mDeltaEta = 0.f; |
999 | 1077 | float mDeltaPhi = 0.f; |
1000 | 1078 |
|
|
0 commit comments