Skip to content

Fixes for modules at negative eta #14716

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

Merged
njacazio merged 3 commits into from
Oct 3, 2025
Merged

Fixes for modules at negative eta #14716

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
33f4a4c
Fixes for modules at negative eta
ekryshen Oct 2, 2025
5189d08
Add braces in if statements
ekryshen Oct 2, 2025
1ab8c4e
retrigger checks
ekryshen Oct 3, 2025
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
1 change: 1 addition & 0 deletions Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/Geometry.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@ class Geometry
double getSamplingAlpha() { return mSamplingAlpha; }
double getCrystalDeltaPhi() { return 2 * std::atan(mCrystalModW / 2 / mRMin); }
double getSamplingDeltaPhi() { return 2 * std::atan(mSamplingModW / 2 / mRMin); }
double getFrontFaceMaxEta(int i);
double getCrystalPhiMin();
double getSamplingPhiMin();
int getNModulesZ() { return mNModulesZ; }
Expand Down
46 changes: 30 additions & 16 deletions Detectors/Upgrades/ALICE3/ECal/base/src/Geometry.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,11 +73,18 @@ double Geometry::getSamplingPhiMin()
return (superModuleDeltaPhi - samplingDeltaPhi * mNSamplingModulesPhi) / 2.;
}

double Geometry::getFrontFaceMaxEta(int i)
{
double theta = std::atan(mRMin / getFrontFaceZatMinR(i));
return -std::log(std::tan(theta / 2.));
}

//==============================================================================
void Geometry::fillFrontFaceCenterCoordinates()
{
if (mFrontFaceCenterR.size() > 0)
if (mFrontFaceCenterR.size() > 0) {
return;
}
mFrontFaceCenterTheta.resize(mNCrystalModulesZ + mNSamplingModulesZ);
mFrontFaceZatMinR.resize(mNCrystalModulesZ + mNSamplingModulesZ);
mFrontFaceCenterR.resize(mNCrystalModulesZ + mNSamplingModulesZ);
Expand Down Expand Up @@ -153,7 +160,7 @@ int Geometry::getCellID(int moduleId, int sectorId, bool isCrystal)
if (sectorId % 2 == 0) { // sampling at positive eta
cellID = sectorId / 2 * mNModulesZ + moduleId + mNSamplingModulesZ + mNCrystalModulesZ * 2;
} else { // sampling at negative eta
cellID = sectorId / 2 * mNModulesZ - moduleId + mNSamplingModulesZ;
cellID = sectorId / 2 * mNModulesZ - moduleId + mNSamplingModulesZ - 1;
}
}
return cellID;
Expand Down Expand Up @@ -206,13 +213,15 @@ void Geometry::detIdToGlobalPosition(int detId, double& x, double& y, double& z)
{
int chamber, sector, iphi, iz;
detIdToRelIndex(detId, chamber, sector, iphi, iz);
double r = 0;
if (iz < mNSamplingModulesZ + mNCrystalModulesZ) {
z = -mFrontFaceCenterZ[mNSamplingModulesZ + mNCrystalModulesZ - iz - 1];
r = mFrontFaceCenterR[mNSamplingModulesZ + mNCrystalModulesZ - iz - 1];
} else {
z = +mFrontFaceCenterZ[iz % (mNSamplingModulesZ + mNCrystalModulesZ)];
z = mFrontFaceCenterZ[iz % (mNSamplingModulesZ + mNCrystalModulesZ)];
r = mFrontFaceCenterR[iz % (mNSamplingModulesZ + mNCrystalModulesZ)];
}
double phi = chamber == 1 ? mFrontFaceCenterCrystalPhi[iphi] : mFrontFaceCenterSamplingPhi[iphi];
double r = mFrontFaceCenterR[iz % (mNSamplingModulesZ + mNCrystalModulesZ)];
x = r * std::cos(phi);
y = r * std::sin(phi);
}
Expand All @@ -224,40 +233,45 @@ int Geometry::areNeighboursVertex(int detId1, int detId2) const
int ch2, sector2, iphi2, iz2;
detIdToRelIndex(detId1, ch1, sector1, iphi1, iz1);
detIdToRelIndex(detId2, ch2, sector2, iphi2, iz2);
if (sector1 != sector2 || ch1 != ch2)
if (sector1 != sector2 || ch1 != ch2) {
return 0;
if (std::abs(iphi1 - iphi2) <= 1 && std::abs(iz1 - iz2) <= 1)
}
if (std::abs(iphi1 - iphi2) <= 1 && std::abs(iz1 - iz2) <= 1) {
return 1;
}
return 0;
}

//==============================================================================
bool Geometry::isAtTheEdge(int cellId)
{
auto [row, col] = globalRowColFromIndex(cellId);
if (col == 0)
if (col == 0) {
return 1;
if (col == mNSamplingModulesZ)
} else if (col == mNSamplingModulesZ) {
return 1;
if (col == mNSamplingModulesZ - 1)
} else if (col == mNSamplingModulesZ - 1) {
return 1;
if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ)
} else if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ) {
return 1;
if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ - 1)
} else if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ - 1) {
return 1;
if (col == mNModulesZ - 1)
} else if (col == mNModulesZ - 1) {
return 1;
}
for (int m = 0; m <= mNSuperModules; m++) {
if (isCrystal(cellId)) {
if (row == m * mNCrystalModulesPhi)
if (row == m * mNCrystalModulesPhi) {
return 1;
if (row == m * mNCrystalModulesPhi - 1)
} else if (row == m * mNCrystalModulesPhi - 1) {
return 1;
}
} else {
if (row == m * mNSamplingModulesPhi)
if (row == m * mNSamplingModulesPhi) {
return 1;
if (row == m * mNSamplingModulesPhi - 1)
} else if (row == m * mNSamplingModulesPhi - 1) {
return 1;
}
}
}
return 0;
Expand Down
42 changes: 25 additions & 17 deletions Detectors/Upgrades/ALICE3/ECal/reconstruction/include/ECalReconstruction/Clusterizer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,25 +48,33 @@ class Clusterizer
void setClusteringThreshold(double threshold) { mClusteringThreshold = threshold; }
void setCrystalDigitThreshold(double threshold) { mCrystalDigitThreshold = threshold; }
void setSamplingDigitThreshold(double threshold) { mSamplingDigitThreshold = threshold; }
void setCrystalEnergyCorrectionPars(std::vector<double> pars) { mCrystalEnergyCorrectionPars = pars; }
void setSamplingEnergyCorrectionPars(std::vector<double> pars) { mSamplingEnergyCorrectionPars = pars; }
void setCrystalZCorrectionPars(std::vector<double> pars) { mCrystalZCorrectionPars = pars; }
void setSamplingZCorrectionPars(std::vector<double> pars) { mSamplingZCorrectionPars = pars; }

private:
std::vector<std::vector<int>> mDigitIndices; // 2D map of digit indices used for recursive cluster finding
bool mUnfoldClusters = true; // to perform cluster unfolding
double mCrystalDigitThreshold = 0.040; // minimal energy of crystal digit
double mSamplingDigitThreshold = 0.100; // minimal energy of sampling digit
double mClusteringThreshold = 0.050; // minimal energy of digit to start clustering (GeV)
double mClusteringTimeGate = 1e9; // maximal time difference between digits to be accepted to clusters (in ns)
int mNLMMax = 30; // maximal number of local maxima in unfolding
double mLogWeight = 4.; // cutoff used in log. weight calculation
double mUnfogingEAccuracy = 1.e-4; // accuracy of energy calculation in unfoding prosedure (GeV)
double mUnfogingXZAccuracy = 1.e-2; // accuracy of position calculation in unfolding procedure (cm)
int mNMaxIterations = 100; // maximal number of iterations in unfolding procedure
double mLocalMaximumCut = 0.015; // minimal height of local maximum over neighbours
bool mApplyCorrectionZ = 1; // z-correction
bool mApplyCorrectionE = 1; // energy-correction
TF1* fCrystalShowerShape; //! Crystal shower shape
TF1* fSamplingShowerShape; //! Sampling shower shape
TF1* fCrystalRMS; //! Crystal RMS
std::vector<std::vector<int>> mDigitIndices; // 2D map of digit indices used for recursive cluster finding
bool mUnfoldClusters = true; // to perform cluster unfolding
double mCrystalDigitThreshold = 0.040; // minimal energy of crystal digit
double mSamplingDigitThreshold = 0.100; // minimal energy of sampling digit
double mClusteringThreshold = 0.050; // minimal energy of digit to start clustering (GeV)
double mClusteringTimeGate = 1e9; // maximal time difference between digits to be accepted to clusters (in ns)
int mNLMMax = 30; // maximal number of local maxima in unfolding
double mLogWeight = 4.; // cutoff used in log. weight calculation
double mUnfogingEAccuracy = 1.e-4; // accuracy of energy calculation in unfoding prosedure (GeV)
double mUnfogingXZAccuracy = 1.e-2; // accuracy of position calculation in unfolding procedure (cm)
int mNMaxIterations = 100; // maximal number of iterations in unfolding procedure
double mLocalMaximumCut = 0.015; // minimal height of local maximum over neighbours
bool mApplyCorrectionZ = 1; // apply z-correction
bool mApplyCorrectionE = 1; // apply energy-correction
TF1* fCrystalShowerShape; //! Crystal shower shape
TF1* fSamplingShowerShape; //! Sampling shower shape
TF1* fCrystalRMS; //! Crystal RMS
std::vector<double> mCrystalEnergyCorrectionPars; // crystal energy-correction parameters
std::vector<double> mSamplingEnergyCorrectionPars; // sampling energy-correction parameters
std::vector<double> mCrystalZCorrectionPars; // crystal z-correction parameters
std::vector<double> mSamplingZCorrectionPars; // sampling z-correction parameters
};

} // namespace ecal
Expand Down
103 changes: 71 additions & 32 deletions Detectors/Upgrades/ALICE3/ECal/reconstruction/src/Clusterizer.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,45 @@ Clusterizer::Clusterizer(bool applyCorrectionZ, bool applyCorrectionE)
mDigitIndices.resize(geo.getNrows(), std::vector<int>(geo.getNcols(), -1));
mApplyCorrectionZ = applyCorrectionZ;
mApplyCorrectionE = applyCorrectionE;

mCrystalEnergyCorrectionPars.reserve(6);
mCrystalEnergyCorrectionPars[0] = 0.00444;
mCrystalEnergyCorrectionPars[1] = -1.322;
mCrystalEnergyCorrectionPars[2] = 1.021;
mCrystalEnergyCorrectionPars[3] = 0.0018;
mCrystalEnergyCorrectionPars[4] = 0.;
mCrystalEnergyCorrectionPars[5] = 0.;

mSamplingEnergyCorrectionPars.reserve(6);
mSamplingEnergyCorrectionPars[0] = 0.0033;
mSamplingEnergyCorrectionPars[1] = -2.09;
mSamplingEnergyCorrectionPars[2] = 1.007;
mSamplingEnergyCorrectionPars[3] = 0.0667;
mSamplingEnergyCorrectionPars[4] = -0.108;
mSamplingEnergyCorrectionPars[5] = 0.0566;

mCrystalZCorrectionPars.reserve(9);
mCrystalZCorrectionPars[0] = -0.005187;
mCrystalZCorrectionPars[1] = 0.7301;
mCrystalZCorrectionPars[2] = -0.7382;
mCrystalZCorrectionPars[3] = 0.;
mCrystalZCorrectionPars[4] = 0.;
mCrystalZCorrectionPars[5] = 0.;
mCrystalZCorrectionPars[6] = 0.;
mCrystalZCorrectionPars[7] = 0.;
mCrystalZCorrectionPars[8] = 0.;

mSamplingZCorrectionPars.reserve(9);
mSamplingZCorrectionPars[0] = -2.137;
mSamplingZCorrectionPars[1] = 6.400;
mSamplingZCorrectionPars[2] = -3.342;
mSamplingZCorrectionPars[3] = -0.1364;
mSamplingZCorrectionPars[4] = 0.4019;
mSamplingZCorrectionPars[5] = -0.1969;
mSamplingZCorrectionPars[6] = 0.008223;
mSamplingZCorrectionPars[7] = -0.02425;
mSamplingZCorrectionPars[8] = 0.01190;

fCrystalShowerShape = new TF1("fCrystal", "x<[1] ? [0]*exp([3]*x+[4]*x*x+[5]*x*x*x) : (x<[2] ? [0]*[6]*exp([7]*x+[8]*x*x) : [0]*[9]*exp([10]*x+[11]*x*x))", 0, 15);
double pc[12];
pc[0] = 1. / 13.15;
Expand Down Expand Up @@ -94,22 +133,24 @@ void Clusterizer::findClusters(const gsl::span<const Digit>& digits, std::vector
void Clusterizer::addDigitToCluster(Cluster& cluster, int row, int col, const gsl::span<const Digit>& digits)
{
auto& geo = Geometry::instance();
if (row < 0 || row >= geo.getNrows() || col < 0 || col >= geo.getNcols())
if (row < 0 || row >= geo.getNrows() || col < 0 || col >= geo.getNcols()) {
return;
}
int digitIndex = mDigitIndices[row][col];
LOGP(debug, " checking row={} and col={} digitIndex={}", row, col, digitIndex);
if (digitIndex < 0)
if (digitIndex < 0) {
return;

}
const Digit& digit = digits[digitIndex];
if (cluster.getMultiplicity() > 0) {
// check if new digit is in the same chamber and sector
const Digit& digit2 = digits[cluster.getDigitIndex(0)];
auto [sector1, ch1] = geo.getSectorChamber(digit.getTower());
auto [sector2, ch2] = geo.getSectorChamber(digit2.getTower());
LOGP(debug, " checking if sector and chamber are the same: ({},{}) ({},{})", sector1, ch1, sector2, ch2);
if (sector1 != sector2 || ch1 != ch2)
if (sector1 != sector2 || ch1 != ch2) {
return;
}
}

mDigitIndices[row][col] = -1;
Expand Down Expand Up @@ -140,11 +181,13 @@ void Clusterizer::makeClusters(const gsl::span<const Digit>& digits, std::vector
auto [row, col] = geo.globalRowColFromIndex(digit.getTower());
bool isCrystal = geo.isCrystal(digit.getTower());
if (isCrystal) {
if (digit.getEnergy() < mCrystalDigitThreshold)
if (digit.getEnergy() < mCrystalDigitThreshold) {
continue;
}
} else {
if (digit.getEnergy() < mSamplingDigitThreshold)
if (digit.getEnergy() < mSamplingDigitThreshold) {
continue;
}
}
mDigitIndices[row][col] = i;
}
Expand All @@ -153,10 +196,12 @@ void Clusterizer::makeClusters(const gsl::span<const Digit>& digits, std::vector
for (int i = 0; i < nDigits; i++) {
const Digit& digitSeed = digits[i];
auto [row, col] = geo.globalRowColFromIndex(digitSeed.getTower());
if (mDigitIndices[row][col] < 0)
if (mDigitIndices[row][col] < 0) {
continue; // digit was already added in one of the clusters
if (digitSeed.getEnergy() < mClusteringThreshold)
}
if (digitSeed.getEnergy() < mClusteringThreshold) {
continue;
}
LOGP(debug, " starting new cluster at row={} and col={}", row, col);
auto& cluster = clusters.emplace_back();
addDigitToCluster(cluster, row, col, digits);
Expand Down Expand Up @@ -343,8 +388,9 @@ void Clusterizer::evalClusters(std::vector<Cluster>& clusters)
double xi, yi, zi;
geo.detIdToGlobalPosition(towerId, xi, yi, zi);
double r = std::sqrt((x - xi) * (x - xi) + (y - yi) * (y - yi) + (z - zi) * (z - zi));
if (r > 2.2)
if (r > 2.2) {
continue;
}
double frac = fCrystalShowerShape->Eval(r);
double rms = fCrystalRMS->Eval(r);
chi2 += std::pow((energy / ee - frac) / rms, 2.);
Expand All @@ -354,38 +400,30 @@ void Clusterizer::evalClusters(std::vector<Cluster>& clusters)

// correct cluster energy and z position
float eta = std::abs(cluster.getEta());
float eCor = 1;
float zCor = 0;
bool isCrystal = geo.isCrystal(cluster.getDigitTowerId(0));
if (isCrystal) {
eCor = 0.00444 * std::pow(ee, -1.322) + (1.021 + 0.0018 * eta);
if (mApplyCorrectionE)
ee *= eCor;
if (mApplyCorrectionZ)
zCor = (-0.00518682 + 0.730052 * eta - 0.73817 * eta * eta);
} else {
eCor = 0.0033 * std::pow(ee, -2.09) + (1.007 + 0.0667 * eta - 0.108 * eta * eta + 0.0566 * eta * eta * eta);
if (mApplyCorrectionE)
ee *= eCor;
if (mApplyCorrectionZ)
zCor = (-2.13679 + 6.40009 * eta - 3.34233 * eta * eta) + (-0.136425 + 0.401887 * eta - 0.196851 * eta * eta) * ee + (0.00822276 - 0.0242512 * eta + 0.0118986 * eta * eta) * ee * ee;
if (mApplyCorrectionE) {
std::vector<double>& pe = isCrystal ? mCrystalEnergyCorrectionPars : mSamplingEnergyCorrectionPars;
ee *= pe[0] * std::pow(ee, pe[1]) + pe[2] + pe[3] * eta + pe[4] * eta * eta + pe[5] * eta * eta * eta;
cluster.setE(ee);
}
if (mApplyCorrectionZ) {
std::vector<double>& pz = isCrystal ? mCrystalZCorrectionPars : mSamplingZCorrectionPars;
float zCor = (pz[0] + pz[1] * eta + pz[2] * eta * eta) + (pz[3] + pz[4] * eta + pz[5] * eta * eta) * ee + (pz[6] + pz[7] * eta + pz[8] * eta * eta) * ee * ee;
cluster.setZ(z > 0 ? z - zCor : z + zCor);
}

cluster.setE(ee);
cluster.setZ(cluster.getZ() - zCor);

// check if cluster is at the edge of detector module
bool isEdge = 0;
for (size_t i = 0; i < cluster.getMultiplicity(); i++) {
int towerId = cluster.getDigitTowerId(i);
if (!geo.isAtTheEdge(towerId))
continue;
isEdge = 1;
break;
if (geo.isAtTheEdge(towerId)) {
isEdge = 1;
break;
}
}
cluster.setEdgeFlag(isEdge);

LOGF(debug, "Cluster coordinates: (%6.2f,%6.2f,%6.2f), eCor=%6.2f zCor=%6.2f", cluster.getX(), cluster.getY(), cluster.getZ(), eCor, zCor);
LOGF(debug, "Cluster coordinates: (%6.2f,%6.2f,%6.2f)", cluster.getX(), cluster.getY(), cluster.getZ());
}
}

Expand All @@ -403,8 +441,9 @@ int Clusterizer::getNumberOfLocalMax(Cluster& clu, int* maxAt, float* maxAtEnerg
for (int i = 0; i < n; i++) {
isLocalMax[i] = false;
float en1 = clu.getDigitEnergy(i);
if (en1 > mClusteringThreshold)
if (en1 > mClusteringThreshold) {
isLocalMax[i] = true;
}
}

for (int i = 0; i < n; i++) {
Expand Down
4 changes: 2 additions & 2 deletions Detectors/Upgrades/ALICE3/ECal/simulation/src/Detector.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -358,9 +358,9 @@ bool Detector::ProcessHits(FairVolume* vol)
return false;
}

if (isCrystal)
if (isCrystal) {
LOGP(debug, "Processing crystal {}", volName.Data());
else {
} else {
eloss *= mSamplingFactorTransportModel;
LOGP(debug, "Processing scintillator {}", volName.Data());
}
Expand Down
8 changes: 5 additions & 3 deletions Detectors/Upgrades/ALICE3/ECal/simulation/src/Digitizer.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,16 +57,18 @@ void Digitizer::processHits(const std::vector<Hit>* hits, std::vector<Digit>& di
bool isCrystal = geo.isCrystal(cellID);
if (isCrystal) { // crystal
double elossSmearedNpe = gRandom->Poisson(eloss * mCrystalPePerGeV) / mCrystalPePerGeV;
if (mSmearCrystal)
if (mSmearCrystal) {
elossSmeared = elossSmearedNpe * gRandom->Gaus(1, 0.007); // light attenuation in crystals
} else { // sampling
}
} else { // sampling
elossSmeared *= mSamplingFraction;
}

Digit& digit = mArrayD[cellID];
digit.setAmplitude(digit.getAmplitude() + elossSmeared);
if (t < digit.getTimeStamp())
if (t < digit.getTimeStamp()) {
digit.setTimeStamp(t); // setting earliest time, TODO: add time smearing
}
LOGF(debug, " crystal: %d cellID = %5d, eloss = %8.5f elossSmeared = %8.5f time = %8.5f", isCrystal, cellID, eloss, elossSmeared, t);

// Adding MC info
Expand Down