Skip to content

Commit 2e27a7b

Browse files
njacazionjacazio
authored
[PWGLF] Add on-the-fly coal to h-d corr (#18091)
1 parent 7de6b9e commit 2e27a7b

1 file changed

Lines changed: 259 additions & 4 deletions

File tree

‎PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx‎

Lines changed: 259 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,7 @@
4646
#include <TParticlePDG.h>
4747
#include <TString.h>
4848

49+
#include <algorithm>
4950
#include <chrono>
5051
#include <cmath>
5152
#include <cstddef>
@@ -91,11 +92,20 @@ struct HadronNucleiCorrelation {
9192
Configurable<bool> doQA{"doQA", true, "save QA histograms"};
9293
Configurable<bool> doMCQA{"doMCQA", false, "save MC QA histograms"};
9394
Configurable<bool> isMC{"isMC", false, "is MC"};
94-
Configurable<bool> isMCGen{"isMCGen", false, "is isMCGen"};
9595
Configurable<bool> isPrim{"isPrim", true, "is isPrim"};
9696
Configurable<bool> doCorrection{"doCorrection", false, "do efficiency correction"};
9797
Configurable<bool> doQuadraticPID{"doQuadraticPID", false, "do PID with sum in quadrature of TOF and TPC"};
9898

99+
struct : ConfigurableGroup {
100+
std::string prefix = "Coalescence"; // JSON group name
101+
Configurable<bool> doMCGenCoalescence{"doMCGenCoalescence", false, "do a simple coalescence on the generated level"};
102+
// Coalescence parameters used in Eur. Phys. J. A 59 (2023) 72:
103+
// p_max = 0.148 GeV/c
104+
// r_max = 2 fm
105+
Configurable<float> pMax{"pMax", 0.148, "maximum momentum for coalescence"};
106+
Configurable<float> rMax{"rMax", 2.0, "maximum radius for coalescence"};
107+
} settingsCoalescence;
108+
99109
Configurable<std::string> fCorrectionPath{"fCorrectionPath", "", "Correction path to file"};
100110
Configurable<std::string> fCorrectionHisto{"fCorrectionHisto", "", "Correction histogram"};
101111
Configurable<std::string> cfgUrl{"cfgUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
@@ -516,11 +526,11 @@ struct HadronNucleiCorrelation {
516526
registry.add("hReco_Pt_Proton_TPCEl", "Reco (anti)protons in reco collisions", {HistType::kTH1F, {ptAxisSmall}});
517527
registry.add("hReco_Pt_Proton_TPCEl_or_TOF", "Reco (anti)protons in reco collisions", {HistType::kTH1F, {ptAxisSmall}});
518528
registry.add("hReco_Pt_Deuteron_TPCEl", "Reco (anti)deuterons in reco collisions", {HistType::kTH1F, {ptAxisSmall}});
519-
registry.add("hReco_Pt_Deuteron_TPCEl_or_TOF", "Reco (anti)protons in reco collisions", {HistType::kTH1F, {ptAxisSmall}});
529+
registry.add("hReco_Pt_Deuteron_TPCEl_or_TOF", "Reco (anti)deuterons in reco collisions", {HistType::kTH1F, {ptAxisSmall}});
520530
}
521531
}
522532

523-
if (isMCGen) {
533+
if (doprocessMixedEventGen || doprocessSameEventGen) {
524534
registry.add("Generated/hNEventsMC", "hNEventsMC", {HistType::kTH1D, {{1, 0.f, 1.f}}});
525535
registry.get<TH1>(HIST("Generated/hNEventsMC"))->GetXaxis()->SetBinLabel(1, "All");
526536

@@ -549,6 +559,8 @@ struct HadronNucleiCorrelation {
549559
registry.add("Generated/hAntiDeuteronsVsPt", "hAntiDeuteronsVsPt", {HistType::kTH1D, {ptAxisGen}});
550560
registry.add("Generated/hProtonsVsPt", "hProtonsVsPt", {HistType::kTH1D, {ptAxisGen}});
551561
registry.add("Generated/hAntiProtonsVsPt", "hAntiProtonsVsPt", {HistType::kTH1D, {ptAxisGen}});
562+
} else if (settingsCoalescence.doMCGenCoalescence.value) {
563+
LOG(fatal) << "This is not a Gen Run but the Gen coalescence is required. Turn off doMCGenCoalescence";
552564
}
553565
}
554566

@@ -1695,10 +1707,253 @@ struct HadronNucleiCorrelation {
16951707

16961708
registry.fill(HIST("Generated/hNEventsMC"), 0.5);
16971709

1710+
// Local representation of a generated particle.
1711+
struct GenCoalescenceCandidate {
1712+
GenCoalescenceCandidate(float pt,
1713+
float eta,
1714+
float phi,
1715+
int pdg,
1716+
float vx,
1717+
float vy,
1718+
float vz,
1719+
float vt,
1720+
bool physicalPrimary)
1721+
: mKinematics{pt, eta, phi},
1722+
mPdg{pdg},
1723+
mVx{vx},
1724+
mVy{vy},
1725+
mVz{vz},
1726+
mVt{vt},
1727+
mPhysicalPrimary{physicalPrimary}
1728+
{
1729+
}
1730+
[[nodiscard]] float pt() const { return mKinematics.pt(); }
1731+
[[nodiscard]] float eta() const { return mKinematics.eta(); }
1732+
[[nodiscard]] float phi() const { return mKinematics.phi(); }
1733+
1734+
[[nodiscard]] float px() const { return pt() * std::cos(phi()); }
1735+
[[nodiscard]] float py() const { return pt() * std::sin(phi()); }
1736+
[[nodiscard]] float pz() const { return mKinematics.pz(); }
1737+
1738+
[[nodiscard]] float vx() const { return mVx; }
1739+
[[nodiscard]] float vy() const { return mVy; }
1740+
[[nodiscard]] float vz() const { return mVz; }
1741+
[[nodiscard]] float vt() const { return mVt; }
1742+
1743+
[[nodiscard]] float mass() const
1744+
{
1745+
switch (std::abs(mPdg)) {
1746+
case PDG_t::kProton:
1747+
return o2::track::PID::getMass(o2::track::PID::Proton);
1748+
case PDG_t::kNeutron:
1749+
return o2::constants::physics::MassNeutron;
1750+
case o2::constants::physics::Pdg::kDeuteron:
1751+
return o2::track::PID::getMass(o2::track::PID::Deuteron);
1752+
default:
1753+
LOG(fatal) << "Unhandled pdg " << mPdg;
1754+
return 0.f;
1755+
}
1756+
}
1757+
[[nodiscard]] float energy() const { return std::hypot(mKinematics.p(), mass()); }
1758+
[[nodiscard]] float rapidity() const { return mKinematics.rapidityForMass(mass()); }
1759+
[[nodiscard]] float y() const { return rapidity(); }
1760+
[[nodiscard]] bool isPhysicalPrimary() const { return mPhysicalPrimary; }
1761+
[[nodiscard]] int pdgCode() const { return mPdg; }
1762+
1763+
private:
1764+
GenCandidate mKinematics;
1765+
1766+
int mPdg = 0;
1767+
float mVx = 0.f;
1768+
float mVy = 0.f;
1769+
float mVz = 0.f;
1770+
float mVt = 0.f;
1771+
1772+
bool mPhysicalPrimary = false;
1773+
};
1774+
1775+
std::vector<GenCoalescenceCandidate> particlesToProcess;
1776+
particlesToProcess.reserve(mcParticles.size());
1777+
1778+
for (const auto& particle : mcParticles) {
1779+
switch (particle.pdgCode()) {
1780+
case PDG_t::kProton:
1781+
case PDG_t::kProtonBar:
1782+
case PDG_t::kNeutron:
1783+
case PDG_t::kNeutronBar:
1784+
case o2::constants::physics::Pdg::kDeuteron:
1785+
case -o2::constants::physics::Pdg::kDeuteron:
1786+
particlesToProcess.emplace_back(particle.pt(),
1787+
particle.eta(),
1788+
particle.phi(),
1789+
particle.pdgCode(),
1790+
particle.vx(),
1791+
particle.vy(),
1792+
particle.vz(),
1793+
particle.vt(),
1794+
particle.isPhysicalPrimary());
1795+
break;
1796+
default:
1797+
break;
1798+
}
1799+
}
1800+
1801+
if (settingsCoalescence.doMCGenCoalescence.value) {
1802+
std::vector<bool> consumed(particlesToProcess.size(), false);
1803+
std::vector<GenCoalescenceCandidate> deuterons;
1804+
1805+
// Try p+n -> d and pbar+nbar -> dbar independently.
1806+
for (const int sign : {+1, -1}) {
1807+
const int protonPDG = sign * PDG_t::kProton;
1808+
const int neutronPDG = sign * PDG_t::kNeutron;
1809+
const int deuteronPDG = sign * o2::constants::physics::Pdg::kDeuteron;
1810+
1811+
for (size_t ip = 0; ip < particlesToProcess.size(); ++ip) {
1812+
if (consumed[ip] || particlesToProcess[ip].pdgCode() != protonPDG) {
1813+
continue;
1814+
}
1815+
1816+
const auto& proton = particlesToProcess[ip];
1817+
1818+
for (size_t in = 0; in < particlesToProcess.size(); ++in) {
1819+
if (consumed[in] || particlesToProcess[in].pdgCode() != neutronPDG) {
1820+
continue;
1821+
}
1822+
1823+
const auto& neutron = particlesToProcess[in];
1824+
1825+
const float ep = proton.energy();
1826+
const float en = neutron.energy();
1827+
1828+
// Propagate the particle freezing out first to the freeze-out time of the other particle
1829+
float xp = proton.vx();
1830+
float yp = proton.vy();
1831+
float zp = proton.vz();
1832+
const float tp = proton.vt();
1833+
1834+
float xn = neutron.vx();
1835+
float yn = neutron.vy();
1836+
float zn = neutron.vz();
1837+
const float tn = neutron.vt();
1838+
1839+
const float commonTime = std::max(static_cast<float>(tp), static_cast<float>(tn));
1840+
1841+
if (tp < commonTime) {
1842+
const float dt = commonTime - tp;
1843+
xp += proton.px() / ep * dt;
1844+
yp += proton.py() / ep * dt;
1845+
zp += proton.pz() / ep * dt;
1846+
}
1847+
1848+
if (tn < commonTime) {
1849+
const float dt = commonTime - tn;
1850+
xn += neutron.px() / en * dt;
1851+
yn += neutron.py() / en * dt;
1852+
zn += neutron.pz() / en * dt;
1853+
}
1854+
1855+
// Velocity of the p-n centre-of-mass frame.
1856+
const float totalE = ep + en;
1857+
const float totalPx = proton.px() + neutron.px();
1858+
const float totalPy = proton.py() + neutron.py();
1859+
const float totalPz = proton.pz() + neutron.pz();
1860+
1861+
const float bx = totalPx / totalE;
1862+
const float by = totalPy / totalE;
1863+
const float bz = totalPz / totalE;
1864+
1865+
const float beta2 = bx * bx + by * by + bz * bz;
1866+
if (beta2 >= 1.) {
1867+
continue;
1868+
}
1869+
1870+
const float gamma = 1. / std::sqrt(1. - beta2);
1871+
1872+
// Boost the proton momentum into the p-n rest frame.
1873+
// In that frame p_p* = -p_n*, therefore |p_p*| is the relative momentum entering the coalescence cut.
1874+
float pxStar = proton.px();
1875+
float pyStar = proton.py();
1876+
float pzStar = proton.pz();
1877+
1878+
if (beta2 > 0.) {
1879+
const float betaDotP = bx * proton.px() + by * proton.py() + bz * proton.pz();
1880+
1881+
const float factor = ((gamma - 1.) * betaDotP / beta2) - gamma * ep;
1882+
1883+
pxStar += factor * bx;
1884+
pyStar += factor * by;
1885+
pzStar += factor * bz;
1886+
}
1887+
1888+
const float pRelative = std::sqrt(pxStar * pxStar + pyStar * pyStar + pzStar * pzStar);
1889+
1890+
if (pRelative >= settingsCoalescence.pMax.value) {
1891+
continue;
1892+
}
1893+
1894+
// The two particles have already been propagated to equal time in the lab. Transform their spatial separation to the pair CM frame.
1895+
float dx = xp - xn;
1896+
float dy = yp - yn;
1897+
float dz = zp - zn;
1898+
1899+
if (beta2 > 0.) {
1900+
const float betaDotR = bx * dx + by * dy + bz * dz;
1901+
const float factor = (gamma - 1.) * betaDotR / beta2;
1902+
1903+
dx += factor * bx;
1904+
dy += factor * by;
1905+
dz += factor * bz;
1906+
}
1907+
1908+
const float rRelative = std::sqrt(dx * dx + dy * dy + dz * dz);
1909+
1910+
if (rRelative >= settingsCoalescence.rMax.value) {
1911+
continue;
1912+
}
1913+
1914+
// Successful coalescence: remove the constituent proton and neutron and replace them by a deuteron carrying their total three-momentum
1915+
consumed[ip] = true;
1916+
consumed[in] = true;
1917+
1918+
const float deuteronPt = std::hypot(totalPx, totalPy);
1919+
const float deuteronPhi = std::atan2(totalPy, totalPx);
1920+
const float deuteronEta = std::asinh(totalPz / std::max(deuteronPt, 1.e-12f));
1921+
1922+
deuterons.emplace_back(deuteronPt,
1923+
deuteronEta,
1924+
deuteronPhi,
1925+
deuteronPDG,
1926+
0.5f * (xp + xn),
1927+
0.5f * (yp + yn),
1928+
0.5f * (zp + zn),
1929+
commonTime,
1930+
proton.isPhysicalPrimary() && neutron.isPhysicalPrimary());
1931+
break;
1932+
}
1933+
}
1934+
}
1935+
1936+
// Construct the post-coalescence particle list.
1937+
std::vector<GenCoalescenceCandidate> coalescedParticles;
1938+
coalescedParticles.reserve(particlesToProcess.size() + deuterons.size());
1939+
1940+
for (size_t i = 0; i < particlesToProcess.size(); ++i) {
1941+
if (!consumed[i]) {
1942+
coalescedParticles.push_back(std::move(particlesToProcess[i]));
1943+
}
1944+
}
1945+
1946+
for (auto& deuteron : deuterons) {
1947+
coalescedParticles.push_back(std::move(deuteron));
1948+
}
1949+
1950+
particlesToProcess = std::move(coalescedParticles);
1951+
}
1952+
16981953
// Pairing candidates are collected during the QA loop below (which already visits every particle)
16991954
GenCollisionCache genCache;
17001955

1701-
for (const auto& particle : mcParticles) {
1956+
for (const auto& particle : particlesToProcess) {
17021957
auto fillGeneratedQa = [this, &particle](const float binPosition) {
17031958
switch (particle.pdgCode()) {
17041959
case PDG_t::kProton:

0 commit comments

Comments
 (0)