Skip to content

Commit f84eacc

Browse files
njacazionjacazio
authored
Update hadronnucleicorrelation.cxx
1 parent 23b2982 commit f84eacc

1 file changed

Lines changed: 98 additions & 67 deletions

File tree

‎PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx‎

Lines changed: 98 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -1708,48 +1708,72 @@ struct HadronNucleiCorrelation {
17081708

17091709
// Local representation of a generated particle.
17101710
struct GenCoalescenceCandidate {
1711-
GenCandidate kinematics;
1712-
1713-
int pdg = 0;
1714-
float x = 0.f;
1715-
float y = 0.f;
1716-
float z = 0.f;
1717-
float t = 0.f;
1718-
1719-
bool physicalPrimary = false;
1720-
1721-
[[nodiscard]] float pt() const { return kinematics.pt(); }
1722-
[[nodiscard]] float eta() const { return kinematics.eta(); }
1723-
[[nodiscard]] float phi() const { return kinematics.phi(); }
1724-
[[nodiscard]] float pz() const { return kinematics.pz(); }
1711+
GenCoalescenceCandidate(float pt,
1712+
float eta,
1713+
float phi,
1714+
int pdg,
1715+
float vx,
1716+
float vy,
1717+
float vz,
1718+
float vt,
1719+
bool physicalPrimary)
1720+
: mKinematics{pt, eta, phi},
1721+
mPdg{pdg},
1722+
mVx{vx},
1723+
mVy{vy},
1724+
mVz{vz},
1725+
mVt{vt},
1726+
mPhysicalPrimary{physicalPrimary}
1727+
{
1728+
}
1729+
[[nodiscard]] float pt() const { return mKinematics.pt(); }
1730+
[[nodiscard]] float eta() const { return mKinematics.eta(); }
1731+
[[nodiscard]] float phi() const { return mKinematics.phi(); }
17251732

17261733
[[nodiscard]] float px() const { return pt() * std::cos(phi()); }
17271734
[[nodiscard]] float py() const { return pt() * std::sin(phi()); }
1735+
[[nodiscard]] float pz() const { return mKinematics.pz(); }
1736+
1737+
[[nodiscard]] float vx() const { return mVx; }
1738+
[[nodiscard]] float vy() const { return mVy; }
1739+
[[nodiscard]] float vz() const { return mVz; }
1740+
[[nodiscard]] float vt() const { return mVt; }
17281741

17291742
[[nodiscard]] float mass() const
17301743
{
1731-
switch (std::abs(pdg)) {
1744+
switch (std::abs(mPdg)) {
17321745
case PDG_t::kProton:
17331746
return o2::track::PID::getMass(o2::track::PID::Proton);
17341747
case PDG_t::kNeutron:
17351748
return o2::constants::physics::MassNeutron;
17361749
case o2::constants::physics::Pdg::kDeuteron:
17371750
return o2::track::PID::getMass(o2::track::PID::Deuteron);
17381751
default:
1739-
LOG(fatal) << "Unhandled pdg " << pdg;
1752+
LOG(fatal) << "Unhandled pdg " << mPdg;
17401753
return 0.f;
17411754
}
17421755
}
1743-
[[nodiscard]] float energy() const { return std::hypot(kinematics.p(), mass()); }
1744-
[[nodiscard]] float rapidity() const { return kinematics.rapidityForMass(mass()); }
1745-
[[nodiscard]] bool isPhysicalPrimary() const { return physicalPrimary; }
1746-
[[nodiscard]] int pdgCode() const { return pdg; }
1756+
[[nodiscard]] float energy() const { return std::hypot(mKinematics.p(), mass()); }
1757+
[[nodiscard]] float rapidity() const { return mKinematics.rapidityForMass(mass()); }
1758+
[[nodiscard]] float y() const { return rapidity(); }
1759+
[[nodiscard]] bool isPhysicalPrimary() const { return mPhysicalPrimary; }
1760+
[[nodiscard]] int pdgCode() const { return mPdg; }
1761+
1762+
private:
1763+
GenCandidate mKinematics;
1764+
1765+
int mPdg = 0;
1766+
float mVx = 0.f;
1767+
float mVy = 0.f;
1768+
float mVz = 0.f;
1769+
float mVt = 0.f;
1770+
1771+
bool mPhysicalPrimary = false;
17471772
};
17481773

17491774
std::vector<GenCoalescenceCandidate> particlesToProcess;
17501775
particlesToProcess.reserve(mcParticles.size());
17511776

1752-
int64_t localIndex = 0;
17531777
for (const auto& particle : mcParticles) {
17541778
switch (particle.pdgCode()) {
17551779
case PDG_t::kProton:
@@ -1758,13 +1782,15 @@ struct HadronNucleiCorrelation {
17581782
case PDG_t::kNeutronBar:
17591783
case o2::constants::physics::Pdg::kDeuteron:
17601784
case -o2::constants::physics::Pdg::kDeuteron:
1761-
particlesToProcess.push_back({.kinematics = {particle.pt(), particle.eta(), particle.phi()},
1762-
.pdg = particle.pdgCode(),
1763-
.x = particle.vx(),
1764-
.y = particle.vy(),
1765-
.z = particle.vz(),
1766-
.t = particle.vt(),
1767-
.physicalPrimary = particle.isPhysicalPrimary()});
1785+
particlesToProcess.emplace_back(particle.pt(),
1786+
particle.eta(),
1787+
particle.phi(),
1788+
particle.pdgCode(),
1789+
particle.vx(),
1790+
particle.vy(),
1791+
particle.vz(),
1792+
particle.vt(),
1793+
particle.isPhysicalPrimary());
17681794
break;
17691795
default:
17701796
break;
@@ -1773,7 +1799,7 @@ struct HadronNucleiCorrelation {
17731799

17741800
if (settingsCoalescence.doMCGenCoalescence.value) {
17751801
std::vector<bool> consumed(particlesToProcess.size(), false);
1776-
std::vector<GenParticle> deuterons;
1802+
std::vector<GenCoalescenceCandidate> deuterons;
17771803

17781804
// Try p+n -> d and pbar+nbar -> dbar independently.
17791805
for (const int sign : {+1, -1}) {
@@ -1799,35 +1825,37 @@ struct HadronNucleiCorrelation {
17991825
const float en = neutron.energy();
18001826

18011827
// Propagate the particle freezing out first to the freeze-out time of the other particle
1802-
float xp = proton.vx;
1803-
float yp = proton.vy;
1804-
float zp = proton.vz;
1805-
1806-
float xn = neutron.vx;
1807-
float yn = neutron.vy;
1808-
float zn = neutron.vz;
1809-
1810-
const float commonTime = std::max(static_cast<float>(proton.vt), static_cast<float>(neutron.vt));
1811-
1812-
if (proton.vt < commonTime) {
1813-
const float dt = commonTime - proton.vt;
1814-
xp += proton.px / ep * dt;
1815-
yp += proton.py / ep * dt;
1816-
zp += proton.pz / ep * dt;
1828+
float xp = proton.vx();
1829+
float yp = proton.vy();
1830+
float zp = proton.vz();
1831+
const float tp = proton.vt();
1832+
1833+
float xn = neutron.vx();
1834+
float yn = neutron.vy();
1835+
float zn = neutron.vz();
1836+
const float tn = neutron.vt();
1837+
1838+
const float commonTime = std::max(static_cast<float>(tp), static_cast<float>(tn));
1839+
1840+
if (tp < commonTime) {
1841+
const float dt = commonTime - tp;
1842+
xp += proton.px() / ep * dt;
1843+
yp += proton.py() / ep * dt;
1844+
zp += proton.pz() / ep * dt;
18171845
}
18181846

1819-
if (neutron.vt < commonTime) {
1820-
const float dt = commonTime - neutron.vt;
1821-
xn += neutron.px / en * dt;
1822-
yn += neutron.py / en * dt;
1823-
zn += neutron.pz / en * dt;
1847+
if (tn < commonTime) {
1848+
const float dt = commonTime - tn;
1849+
xn += neutron.px() / en * dt;
1850+
yn += neutron.py() / en * dt;
1851+
zn += neutron.pz() / en * dt;
18241852
}
18251853

18261854
// Velocity of the p-n centre-of-mass frame.
18271855
const float totalE = ep + en;
1828-
const float totalPx = proton.px + neutron.px;
1829-
const float totalPy = proton.py + neutron.py;
1830-
const float totalPz = proton.pz + neutron.pz;
1856+
const float totalPx = proton.px() + neutron.px();
1857+
const float totalPy = proton.py() + neutron.py();
1858+
const float totalPz = proton.pz() + neutron.pz();
18311859

18321860
const float bx = totalPx / totalE;
18331861
const float by = totalPy / totalE;
@@ -1842,12 +1870,12 @@ struct HadronNucleiCorrelation {
18421870

18431871
// Boost the proton momentum into the p-n rest frame.
18441872
// In that frame p_p* = -p_n*, therefore |p_p*| is the relative momentum entering the coalescence cut.
1845-
float pxStar = proton.px;
1846-
float pyStar = proton.py;
1847-
float pzStar = proton.pz;
1873+
float pxStar = proton.px();
1874+
float pyStar = proton.py();
1875+
float pzStar = proton.pz();
18481876

18491877
if (beta2 > 0.) {
1850-
const float betaDotP = bx * proton.px + by * proton.py + bz * proton.pz;
1878+
const float betaDotP = bx * proton.px() + by * proton.py() + bz * proton.pz();
18511879

18521880
const float factor = ((gamma - 1.) * betaDotP / beta2) - gamma * ep;
18531881

@@ -1886,23 +1914,26 @@ struct HadronNucleiCorrelation {
18861914
consumed[ip] = true;
18871915
consumed[in] = true;
18881916

1889-
deuterons.push_back({deuteronPDG,
1890-
static_cast<float>(totalPx),
1891-
static_cast<float>(totalPy),
1892-
static_cast<float>(totalPz),
1893-
static_cast<float>(0.5 * (xp + xn)),
1894-
static_cast<float>(0.5 * (yp + yn)),
1895-
static_cast<float>(0.5 * (zp + zn)),
1896-
static_cast<float>(commonTime),
1897-
proton.isPhysicalPrimary() && neutron.isPhysicalPrimary(),
1898-
localIndex++});
1917+
const float deuteronPt = std::hypot(totalPx, totalPy);
1918+
const float deuteronPhi = std::atan2(totalPy, totalPx);
1919+
const float deuteronEta = std::asinh(totalPz / std::max(deuteronPt, 1.e-12f));
1920+
1921+
deuterons.emplace_back(deuteronPt,
1922+
deuteronEta,
1923+
deuteronPhi,
1924+
deuteronPDG,
1925+
0.5f * (xp + xn),
1926+
0.5f * (yp + yn),
1927+
0.5f * (zp + zn),
1928+
commonTime,
1929+
proton.isPhysicalPrimary() && neutron.isPhysicalPrimary());
18991930
break;
19001931
}
19011932
}
19021933
}
19031934

19041935
// Construct the post-coalescence particle list.
1905-
std::vector<GenParticle> coalescedParticles;
1936+
std::vector<GenCoalescenceCandidate> coalescedParticles;
19061937
coalescedParticles.reserve(particlesToProcess.size() + deuterons.size());
19071938

19081939
for (size_t i = 0; i < particlesToProcess.size(); ++i) {

0 commit comments

Comments
 (0)