Skip to content

[DPG] New method to select primaries #17248

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rliotino99 wants to merge 2 commits into
base: master
Choose a base branch
from
+130 −56
Open

[DPG] New method to select primaries #17248

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6e9b0e7
New method to select primaries
rliotino99 Jul 28, 2026
199d2ed
Fixing code-check
rliotino99 Jul 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
186 changes: 130 additions & 56 deletions DPG/Tasks/AOTTrack/PID/HMPID/hmpidTableProducer.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,8 @@
#include <Framework/runDataProcessing.h>

#include <TGeoManager.h>
#include <TMath.h>
#include <TPDGCode.h>

#include <HMPIDBase/Param.h>

Expand Down Expand Up @@ -73,6 +75,12 @@ struct HmpidTableProducer {
Configurable<bool> requireTPC{"requireTPC", true, "Require TPC track"};
Configurable<bool> requireTOF{"requireTOF", true, "Require TOF track"};

Configurable<bool> useInAbsorberGeomMethod{"useInAbsorberGeomMethod", false, "Use geometrical method to check if daughters are born in absorber"};

// (reference) 473 cm - was the legacy value in run2 simulation
Configurable<float> survivalThresholdRich2{"survivalThresholdRich2", 437.5f, "survivalThresholdRich2"};
Configurable<float> survivalThresholdRich4{"survivalThresholdRich4", 439.0f, "survivalThresholdRich4"};

using CollisionCandidates = o2::soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFV0As>;

using TrackCandidates = soa::Join<aod::Tracks, aod::TracksExtra,
Expand All @@ -94,30 +102,9 @@ struct HmpidTableProducer {

static constexpr int Rich2 = 2, Rich4 = 4;

// -----------------------------------------------------------------------
// HMPID absorber geometry (hardcoded from HMPIDSimulation/Detector.cxx,
// Detector::ConstructGeometry / Detector::createAbsorber).
// The experiment is finalised and this geometry will not change, so the
// values are copied here instead of being re-derived from TGeoManager/CCDB
// at runtime. If the detector geometry code is ever revisited, these
// constants must be updated accordingly.
//
// Each absorber is a box (TGeoBBox) whose LOCAL->GLOBAL transform is built as:
// pMatrix->SetTranslation(T);
// pMatrix->RotateZ(theta);
// which yields, for a local point p: p_glob = Rz(theta) * p_loc + T.
// In particular the box CENTER in global coordinates is exactly T (the
// rotation does not affect T, since it is applied to p_loc only, not to
// the already-set translation). Only the box AXES are rotated by theta
// with respect to the global x,y axes.
//
// To test whether a global point lies inside the box we invert the
// transform: p_loc = Rz(-theta) * (p_glob - T), then compare component-wise
// against the box half-widths.
// -----------------------------------------------------------------------
// (reference) HMPID Detector class in O2
static constexpr double AbsThetaDeg = 33.5;
const double mAbsCosT = std::cos(AbsThetaDeg * TMath::DegToRad());
const double mAbsSinT = std::sin(AbsThetaDeg * TMath::DegToRad());
double mAbsCosT = 0., mAbsSinT = 0.;

// Rich2 absorber: trans2 = {435.5, 0., -155.}, thickness 40mm -> halfX = 2cm
static constexpr double AbsRich2CenterX = 435.5, AbsRich2CenterZ = -155.;
Expand All @@ -132,6 +119,9 @@ struct HmpidTableProducer {

void init(o2::framework::InitContext&)
{
mAbsCosT = std::cos(AbsThetaDeg * TMath::DegToRad());
mAbsSinT = std::sin(AbsThetaDeg * TMath::DegToRad());

ccdb->setURL(ccdbConfig.ccdbUrl);
ccdb->setCaching(true);
ccdb->setLocalObjectValidityChecking();
Expand Down Expand Up @@ -159,6 +149,10 @@ struct HmpidTableProducer {
kTH1F, {{4, -0.5, 3.5, ""}});

histos.add("hProdVertex", ";X (cm);Y (cm);Z (cm)", HistType::kTH3F, {{500, -500., 500.}, {500, -500., 500.}, {500, -500., 500.}});
histos.add("hDaughterRCyl_Rich2", "hDaughterRCyl_Rich2", kTH1F, {{600, 0., 600.}});
histos.add("hDaughterRCyl_Rich4", "hDaughterRCyl_Rich4", kTH1F, {{600, 0., 600.}});
histos.add("hDaughterRSph_Rich2", "hDaughterRSph_Rich2", kTH1F, {{600, 0., 600.}});
histos.add("hDaughterRSph_Rich4", "hDaughterRSph_Rich4", kTH1F, {{600, 0., 600.}});
}

// -----------------------------------------------------------------------
Expand Down Expand Up @@ -274,8 +268,9 @@ struct HmpidTableProducer {
(x[1] - planePoint[1]) * planeNormal[1] +
(x[2] - planePoint[2]) * planeNormal[2];

if (std::abs(dist) >= std::abs(distPrev))
if (std::abs(dist) >= std::abs(distPrev)) {
return false;
}

distPrev = dist;
s -= dist;
Expand Down Expand Up @@ -318,14 +313,16 @@ struct HmpidTableProducer {
// Intersection track - radiator plane
std::array<double, 3> xRad{}, pAtRad{};

if (!intersectHelixPlane(bz, charge, x, p, pRad, nRad, xRad, pAtRad))
if (!intersectHelixPlane(bz, charge, x, p, pRad, nRad, xRad, pAtRad)) {
continue;
}

// Intersection track - PC plane
std::array<double, 3> xPc{}, pAtPc{};

if (!intersectHelixPlane(bz, charge, xRad, pAtRad, pPc, nPc, xPc, pAtPc))
if (!intersectHelixPlane(bz, charge, xRad, pAtRad, pPc, nPc, xPc, pAtPc)) {
continue;
}

double theta = 0., phi = 0.;
param->mars2LorsVec(ch, pAtRad.data(), theta, phi);
Expand All @@ -334,8 +331,9 @@ struct HmpidTableProducer {
param->mars2Lors(ch, xPc.data(), xL, yL);

// Use isInside to check Chamber intersected
if (param->isInside(xL, yL, param->distCut()))
if (param->isInside(xL, yL, param->distCut())) {
return ch;
}
}

// No chamber intersected
Expand All @@ -361,17 +359,32 @@ struct HmpidTableProducer {
return false;
}

// subtract the box center (translation is not rotated, see geometry block above)
const double rx = vx * mAbsCosT + vy * mAbsSinT;
const double ry = -vx * mAbsSinT + vy * mAbsCosT;
const double rz = vz;
// translate to box center
const double lx = vx - centerX;
const double ly = vy; // centerY = 0
const double lz = vz - centerZ;

// rotate by -theta into the box local frame
const double lx = rx - centerX;
const double ly = ry; // centerY = 0
const double lz = rz - centerZ;
const double rx = lx * mAbsCosT + ly * mAbsSinT;
const double ry = -lx * mAbsSinT + ly * mAbsCosT;
const double rz = lz;

return std::abs(rx) <= halfX && std::abs(ry) <= AbsHalfY && std::abs(rz) <= AbsHalfZ;
}

bool survivedAbsorber(double vx, double vy, int chamber)
{
float thresholdR = 0.;
if (chamber == Rich2) {
thresholdR = survivalThresholdRich2;
} else if (chamber == Rich4) {
thresholdR = survivalThresholdRich4;
} else {
return false;
}

return std::abs(lx) <= halfX && std::abs(ly) <= AbsHalfY && std::abs(lz) <= AbsHalfZ;
const float r = std::hypot(vx, vy);
return r > thresholdR;
}

void processEvent(CollisionCandidates::iterator const& col,
Expand All @@ -397,26 +410,30 @@ struct HmpidTableProducer {

const auto& globalTrack = t.template track_as<TTrackTable>();

if (!globalTrack.has_collision())
if (!globalTrack.has_collision()) {
continue;
}

const auto& col = globalTrack.template collision_as<CollisionCandidates>();
initCCDB(col.template bc_as<aod::BCsWithTimestamps>());
uint32_t collId = col.globalIndex();

if ((requireITS && !globalTrack.hasITS()) ||
(requireTPC && !globalTrack.hasTPC()) ||
(requireTOF && !globalTrack.hasTOF()))
(requireTOF && !globalTrack.hasTOF())) {
continue;
}

if (mCollisionsWithHmpid.insert(collId).second)
if (mCollisionsWithHmpid.insert(collId).second) {
histos.fill(HIST("eventsHmpid"), 0.5);
}

// clusSize diagnostics
histos.fill(HIST("hClusSize"), t.hmpidClusSize());
bool isCorrupt = (t.hmpidClusSize() <= 0);
if (isCorrupt)
if (isCorrupt) {
histos.fill(HIST("hClusSizeCorrupt"), t.hmpidClusSize());
}

// --- M2: clusSize encoding ---
int chamberM2 = t.hmpidClusSize() / 1000000;
Expand Down Expand Up @@ -464,20 +481,22 @@ struct HmpidTableProducer {
// bin 2 = clusSize <= 0, M1 recovery (corrupt, M1 ok)
// bin 3 = clusSize <= 0, M1 fails (corrupt, skipped)

if (!isCorrupt && chamberM3 >= 0)
if (!isCorrupt && chamberM3 >= 0) {
histos.fill(HIST("hChamberAssignment"), 0.);
else if (!isCorrupt && chamberM3 < 0)
} else if (!isCorrupt && chamberM3 < 0) {
histos.fill(HIST("hChamberAssignment"), 1.);
else if (isCorrupt && chamberM3 >= 0)
} else if (isCorrupt && chamberM3 >= 0) {
histos.fill(HIST("hChamberAssignment"), 2.);
else
} else {
histos.fill(HIST("hChamberAssignment"), 3.);
}

histos.fill(HIST("hChamberM3"), chamberM3);
histos.fill(HIST("hChamberM3vsM2"), chamberM2, chamberM3);

if (chamberM3 < 0)
if (chamberM3 < 0) {
continue;
}

std::vector<float> hmpidPhotsCharge2(o2::aod::kDimPhotonsCharge, 0.f);

Expand Down Expand Up @@ -514,23 +533,78 @@ struct HmpidTableProducer {

if ((chamberM3 == Rich2 || chamberM3 == Rich4) && mc.has_daughters()) {
auto dIds = mc.daughtersIds();

for (int32_t idx = dIds.front(); idx <= dIds.back(); ++idx) {
auto daughter = mcParticles.rawIteratorAt(idx);

histos.fill(HIST("hProdVertex"), daughter.vx(), daughter.vy(), daughter.vz());

if (isInAbsorber(daughter.vx(), daughter.vy(), daughter.vz(), chamberM3)) {
interactionInAbsorber = true;
break;
}
} // end loop daughters
bool foundRelevantDaughter = false; // true if at least one non-delta/photon daughter was examined

if (useInAbsorberGeomMethod) {
for (int32_t idx = dIds.front(); idx <= dIds.back(); ++idx) {
auto daughter = mcParticles.rawIteratorAt(idx);

int absPdg = std::abs(daughter.pdgCode());
if (absPdg == kElectron || absPdg == kGamma) {
continue;
}

foundRelevantDaughter = true;

// diagnostics on daughters distribution
histos.fill(HIST("hProdVertex"), daughter.vx(), daughter.vy(), daughter.vz());

double rCyl = std::hypot(daughter.vx(), daughter.vy());
double rSph = std::hypot(daughter.vx(), daughter.vy(), daughter.vz());
if (chamberM3 == Rich2) {
histos.fill(HIST("hDaughterRCyl_Rich2"), rCyl);
histos.fill(HIST("hDaughterRSph_Rich2"), rSph);
} else {
histos.fill(HIST("hDaughterRCyl_Rich4"), rCyl);
histos.fill(HIST("hDaughterRSph_Rich4"), rSph);
}

if (isInAbsorber(daughter.vx(), daughter.vy(), daughter.vz(), chamberM3)) {
interactionInAbsorber = true;
}
} // end loop daughters
} else {
bool survived = false;
for (int32_t idx = dIds.front(); idx <= dIds.back(); ++idx) {
auto daughter = mcParticles.rawIteratorAt(idx);

// skip delta rays (e-/e+), photons, and HMPID Cherenkov/feedback
int absPdg = std::abs(daughter.pdgCode());
if (absPdg == kElectron || absPdg == kGamma) {
continue;
}

foundRelevantDaughter = true;

histos.fill(HIST("hProdVertex"), daughter.vx(), daughter.vy(), daughter.vz());

double rCyl = std::hypot(daughter.vx(), daughter.vy());
double rSph = std::hypot(daughter.vx(), daughter.vy(), daughter.vz());
if (chamberM3 == Rich2) {
histos.fill(HIST("hDaughterRCyl_Rich2"), rCyl);
histos.fill(HIST("hDaughterRSph_Rich2"), rSph);
} else {
histos.fill(HIST("hDaughterRCyl_Rich4"), rCyl);
histos.fill(HIST("hDaughterRSph_Rich4"), rSph);
}

if (survivedAbsorber(daughter.vx(), daughter.vy(), chamberM3)) {
survived = true;
}
} // end loop daughters

// No relevant daughter found (only delta rays/photons, or no
// daughters at all): no evidence of a genuine interaction -> treat
// as primary/survived, consistent with the "no daughters" case.
interactionInAbsorber = foundRelevantDaughter ? !survived : false;
}
} // end if has_daughters

hmpidAnalysisMC(mc.pdgCode(), mc.vx(), mc.vy(), mc.vz(),
mc.isPhysicalPrimary(), mc.getProcess(), interactionInAbsorber);
} else {
hmpidAnalysisMC(-1, 0.f, 0.f, 0.f, false, -100, false);
// No MC truth associated to this track
hmpidAnalysisMC(-999, -999.f, -999.f, -999.f, false, -100, false);
}
} // end if constexpr (isMC)

Expand Down
Loading