Skip to content

Commit 1d65360

Browse files
njacazionjacazio
authored
[PWGLF] h-d corr. task add deltaY plots, extend axis range (#18001)
1 parent 320e288 commit 1d65360

1 file changed

Lines changed: 58 additions & 28 deletions

File tree

‎PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx‎

Lines changed: 58 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@
3636
#include <Framework/InitContext.h>
3737
#include <Framework/O2DatabasePDGPlugin.h>
3838
#include <Framework/runDataProcessing.h>
39+
#include <ReconstructionDataFormats/PID.h>
3940

4041
#include <TH1.h>
4142
#include <TH2.h>
@@ -47,6 +48,7 @@
4748

4849
#include <chrono>
4950
#include <cmath>
51+
#include <cstddef>
5052
#include <cstdint>
5153
#include <memory>
5254
#include <string>
@@ -146,6 +148,12 @@ struct HadronNucleiCorrelation {
146148
ConfigurableAxis axisNSigma{"axisNSigma", {35, -7.f, 7.f}, "n#sigma"};
147149
ConfigurableAxis deltaPhiAxis = {"deltaPhiAxis", {46, -1 * o2::constants::math::PIHalf, 3 * o2::constants::math::PIHalf}, "#Delta#phi (rad)"};
148150

151+
struct : ConfigurableGroup {
152+
std::string prefix = "axes"; // JSON group name
153+
ConfigurableAxis axisDeltaEta{"axisDeltaEta", {300, -1.5, 1.5}, "#Delta#eta"};
154+
ConfigurableAxis axisDeltaRap{"axisDeltaRap", {300, -1.5, 1.5}, "#Delta y"};
155+
} settingsAxes;
156+
149157
using FilteredCollisions = soa::Filtered<aod::SingleCollSels>;
150158
using FilteredCollisionsExtra = soa::Filtered<soa::Join<aod::SingleCollSels, aod::SingleCollExtras>>;
151159
using SimCollisions = soa::Filtered<aod::McCollisions>;
@@ -165,10 +173,10 @@ struct HadronNucleiCorrelation {
165173
// std::unique_ptr<o2::aod::singletrackselector::FemtoPair<TrkTypeMC>> PairMC = std::make_unique<o2::aod::singletrackselector::FemtoPair<TrkTypeMC>>();
166174

167175
// Data histograms
168-
std::vector<std::shared_ptr<TH3>> hEtaPhiSameEv;
169-
std::vector<std::shared_ptr<TH3>> hEtaPhiMixdEv;
170-
std::vector<std::shared_ptr<TH3>> hCorrEtaPhiSameEv;
171-
std::vector<std::shared_ptr<TH3>> hCorrEtaPhiMixdEv;
176+
std::vector<std::shared_ptr<TH3>> hDeltaPhiSameEv;
177+
std::vector<std::shared_ptr<TH3>> hDeltaPhiMixdEv;
178+
std::vector<std::shared_ptr<TH3>> hCorrDeltaPhiSameEv;
179+
std::vector<std::shared_ptr<TH3>> hCorrDeltaPhiMixdEv;
172180

173181
int nBinspT = 0;
174182
TH2* hEffPtEtaProton = nullptr;
@@ -185,9 +193,32 @@ struct HadronNucleiCorrelation {
185193
// compact candidate vectors instead of re-reading the full MC particle table for every pair
186194
struct GenCandidate {
187195
float ptVal, etaVal, phiVal;
188-
float pt() const { return ptVal; }
189-
float eta() const { return etaVal; }
190-
float phi() const { return phiVal; }
196+
[[nodiscard]] float pt() const { return ptVal; }
197+
[[nodiscard]] float eta() const { return etaVal; }
198+
[[nodiscard]] float phi() const { return phiVal; }
199+
[[nodiscard]] float pz() const { return ptVal * std::sinh(etaVal); }
200+
[[nodiscard]] float p() const { return std::hypot(ptVal, pz()); }
201+
[[nodiscard]] float energyForMass(const float mass) const { return std::hypot(p(), mass); }
202+
[[nodiscard]] float rapidityForMass(const float mass) const
203+
{
204+
const float e = energyForMass(mass);
205+
return 0.5f * std::log((e + pz()) / (e - pz()));
206+
};
207+
[[nodiscard]] float mass(const int pdg) const
208+
{
209+
switch (pdg) {
210+
case o2::constants::physics::Pdg::kDeuteron:
211+
case -o2::constants::physics::Pdg::kDeuteron:
212+
return o2::track::PID::getMass(o2::track::PID::Deuteron);
213+
case PDG_t::kProton:
214+
case -PDG_t::kProton:
215+
return o2::track::PID::getMass(o2::track::PID::Proton);
216+
default:
217+
LOG(fatal) << "Unhandled pdg " << pdg;
218+
return 0.f;
219+
}
220+
}
221+
[[nodiscard]] float rapidityForPdg(const int pdg) const { return rapidityForMass(mass(pdg)); }
191222
};
192223

193224
// Per-collision information needed for generated-level same- and mixed-event pairing
@@ -223,8 +254,8 @@ struct HadronNucleiCorrelation {
223254
const AxisSpec ptAxis = {200, -10.f, 10.f, "#it{p}_{T} GeV/#it{c}"};
224255
const AxisSpec ptAxisSmall = {100, -5.f, 5.f, "#it{p}_{T} GeV/#it{c}"};
225256

226-
const AxisSpec deltaEtaAxis = {300, -1.5, 1.5, "#Delta#eta"};
227-
const AxisSpec deltaRapAxis = {300, -1.5, 1.5, "#Delta y"};
257+
const AxisSpec deltaEtaAxis = {settingsAxes.axisDeltaEta, "#Delta#eta"};
258+
const AxisSpec deltaRapAxis = {settingsAxes.axisDeltaRap, "#Delta y"};
228259

229260
if (doprocessSameEvent || doprocessSameEventEvSel) {
230261
registry.add("hNEvents", "hNEvents", {HistType::kTH1D, {{7, 0.f, 7.f}}});
@@ -288,17 +319,17 @@ struct HadronNucleiCorrelation {
288319
const TString ptTag = Form("pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10);
289320
const TString ptInterval = Form("(%.1f<p_{T}^{assoc} <%.1f GeV/c)", pTBins.value.at(i), pTBins.value.at(i + 1));
290321
if (doRapidity) {
291-
hEtaPhiSameEv.push_back(registry.add<TH3>(Form("hEtaPhi_%s_SE_%s", name.Data(), ptTag.Data()), "Raw #Delta y #Delta#phi " + ptInterval, {HistType::kTH3F, {deltaRapAxis, deltaPhiAxis, ptBinnedAxis}}));
292-
hEtaPhiMixdEv.push_back(registry.add<TH3>(Form("hEtaPhi_%s_ME_%s", name.Data(), ptTag.Data()), "Raw #Delta y #Delta#phi " + ptInterval, {HistType::kTH3F, {deltaRapAxis, deltaPhiAxis, ptBinnedAxis}}));
322+
hDeltaPhiSameEv.push_back(registry.add<TH3>(Form("hEtaPhi_%s_SE_%s", name.Data(), ptTag.Data()), "Raw #Delta y #Delta#phi " + ptInterval, {HistType::kTH3F, {deltaRapAxis, deltaPhiAxis, ptBinnedAxis}}));
323+
hDeltaPhiMixdEv.push_back(registry.add<TH3>(Form("hEtaPhi_%s_ME_%s", name.Data(), ptTag.Data()), "Raw #Delta y #Delta#phi " + ptInterval, {HistType::kTH3F, {deltaRapAxis, deltaPhiAxis, ptBinnedAxis}}));
293324

294-
hCorrEtaPhiSameEv.push_back(registry.add<TH3>(Form("hCorrEtaPhi_%s_SE_%s", name.Data(), ptTag.Data()), "#Delta y #Delta#phi " + ptInterval, {HistType::kTH3F, {deltaRapAxis, deltaPhiAxis, ptBinnedAxis}}));
295-
hCorrEtaPhiMixdEv.push_back(registry.add<TH3>(Form("hCorrEtaPhi_%s_ME_%s", name.Data(), ptTag.Data()), "#Delta y #Delta#phi " + ptInterval, {HistType::kTH3F, {deltaRapAxis, deltaPhiAxis, ptBinnedAxis}}));
325+
hCorrDeltaPhiSameEv.push_back(registry.add<TH3>(Form("hCorrEtaPhi_%s_SE_%s", name.Data(), ptTag.Data()), "#Delta y #Delta#phi " + ptInterval, {HistType::kTH3F, {deltaRapAxis, deltaPhiAxis, ptBinnedAxis}}));
326+
hCorrDeltaPhiMixdEv.push_back(registry.add<TH3>(Form("hCorrEtaPhi_%s_ME_%s", name.Data(), ptTag.Data()), "#Delta y #Delta#phi " + ptInterval, {HistType::kTH3F, {deltaRapAxis, deltaPhiAxis, ptBinnedAxis}}));
296327
} else {
297-
hEtaPhiSameEv.push_back(registry.add<TH3>(Form("hEtaPhi_%s_SE_%s", name.Data(), ptTag.Data()), "Raw #Delta#eta#Delta#phi " + ptInterval, {HistType::kTH3F, {deltaEtaAxis, deltaPhiAxis, ptBinnedAxis}}));
298-
hEtaPhiMixdEv.push_back(registry.add<TH3>(Form("hEtaPhi_%s_ME_%s", name.Data(), ptTag.Data()), "Raw #Delta#eta#Delta#phi " + ptInterval, {HistType::kTH3F, {deltaEtaAxis, deltaPhiAxis, ptBinnedAxis}}));
328+
hDeltaPhiSameEv.push_back(registry.add<TH3>(Form("hEtaPhi_%s_SE_%s", name.Data(), ptTag.Data()), "Raw #Delta#eta#Delta#phi " + ptInterval, {HistType::kTH3F, {deltaEtaAxis, deltaPhiAxis, ptBinnedAxis}}));
329+
hDeltaPhiMixdEv.push_back(registry.add<TH3>(Form("hEtaPhi_%s_ME_%s", name.Data(), ptTag.Data()), "Raw #Delta#eta#Delta#phi " + ptInterval, {HistType::kTH3F, {deltaEtaAxis, deltaPhiAxis, ptBinnedAxis}}));
299330

300-
hCorrEtaPhiSameEv.push_back(registry.add<TH3>(Form("hCorrEtaPhi_%s_SE_%s", name.Data(), ptTag.Data()), "#Delta#eta#Delta#phi " + ptInterval, {HistType::kTH3F, {deltaEtaAxis, deltaPhiAxis, ptBinnedAxis}}));
301-
hCorrEtaPhiMixdEv.push_back(registry.add<TH3>(Form("hCorrEtaPhi_%s_ME_%s", name.Data(), ptTag.Data()), "#Delta#eta#Delta#phi " + ptInterval, {HistType::kTH3F, {deltaEtaAxis, deltaPhiAxis, ptBinnedAxis}}));
331+
hCorrDeltaPhiSameEv.push_back(registry.add<TH3>(Form("hCorrEtaPhi_%s_SE_%s", name.Data(), ptTag.Data()), "#Delta#eta#Delta#phi " + ptInterval, {HistType::kTH3F, {deltaEtaAxis, deltaPhiAxis, ptBinnedAxis}}));
332+
hCorrDeltaPhiMixdEv.push_back(registry.add<TH3>(Form("hCorrEtaPhi_%s_ME_%s", name.Data(), ptTag.Data()), "#Delta#eta#Delta#phi " + ptInterval, {HistType::kTH3F, {deltaEtaAxis, deltaPhiAxis, ptBinnedAxis}}));
302333
}
303334
}
304335
}
@@ -640,14 +671,14 @@ struct HadronNucleiCorrelation {
640671
}
641672

642673
if (ME) {
643-
hEtaPhiMixdEv[k]->Fill(deltaEta, deltaPhi, part1.pt());
674+
hDeltaPhiMixdEv[k]->Fill(deltaEta, deltaPhi, part1.pt());
644675
if (corr0 != 0 && corr1 != 0) {
645-
hCorrEtaPhiMixdEv[k]->Fill(deltaEta, deltaPhi, part1.pt(), 1. / (corr0 * corr1));
676+
hCorrDeltaPhiMixdEv[k]->Fill(deltaEta, deltaPhi, part1.pt(), 1. / (corr0 * corr1));
646677
}
647678
} else {
648-
hEtaPhiSameEv[k]->Fill(deltaEta, deltaPhi, part1.pt());
679+
hDeltaPhiSameEv[k]->Fill(deltaEta, deltaPhi, part1.pt());
649680
if (corr0 != 0 && corr1 != 0) {
650-
hCorrEtaPhiSameEv[k]->Fill(deltaEta, deltaPhi, part1.pt(), 1. / (corr0 * corr1));
681+
hCorrDeltaPhiSameEv[k]->Fill(deltaEta, deltaPhi, part1.pt(), 1. / (corr0 * corr1));
651682
}
652683
} // SE
653684
} // pT condition
@@ -663,17 +694,16 @@ struct HadronNucleiCorrelation {
663694
float deltaEta = part0.eta() - part1.eta();
664695
float deltaPhi = part0.phi() - part1.phi();
665696
deltaPhi = RecoDecay::constrainAngle(deltaPhi, -1 * o2::constants::math::PIHalf);
666-
697+
// Here we have to use doRapidity
698+
const float deltaRapidity = part0.rapidityForPdg(pdgPart0) - part1.rapidityForPdg(pdgPart1);
667699
for (int k = 0; k < nBinspT; k++) {
668-
669700
if (part0.pt() >= pTBins.value.at(k) && part0.pt() < pTBins.value.at(k + 1)) {
670-
671701
if (ME) {
672-
hEtaPhiMixdEv[k]->Fill(deltaEta, deltaPhi, part1.pt());
673-
hCorrEtaPhiMixdEv[k]->Fill(deltaEta, deltaPhi, part1.pt());
702+
hDeltaPhiMixdEv[k]->Fill(doRapidity ? deltaRapidity : deltaEta, deltaPhi, part1.pt());
703+
hCorrDeltaPhiMixdEv[k]->Fill(doRapidity ? deltaRapidity : deltaEta, deltaPhi, part1.pt());
674704
} else {
675-
hEtaPhiSameEv[k]->Fill(deltaEta, deltaPhi, part1.pt());
676-
hCorrEtaPhiSameEv[k]->Fill(deltaEta, deltaPhi, part1.pt());
705+
hDeltaPhiSameEv[k]->Fill(doRapidity ? deltaRapidity : deltaEta, deltaPhi, part1.pt());
706+
hCorrDeltaPhiSameEv[k]->Fill(doRapidity ? deltaRapidity : deltaEta, deltaPhi, part1.pt());
677707
} // SE
678708
} // pT condition
679709
} // nBinspT loop

0 commit comments

Comments
 (0)