Skip to content

Commit b1904ff

Browse files
ekryshenekryshen
authored
Fixes for modules at negative eta (#14716)
* Fixes for modules at negative eta * Add braces in if statements * retrigger checks
1 parent a430fb1 commit b1904ff

File tree

6 files changed

+134
-70
lines changed

6 files changed

+134
-70
lines changed

‎Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/Geometry.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,7 @@ class Geometry
6161
double getSamplingAlpha() { return mSamplingAlpha; }
6262
double getCrystalDeltaPhi() { return 2 * std::atan(mCrystalModW / 2 / mRMin); }
6363
double getSamplingDeltaPhi() { return 2 * std::atan(mSamplingModW / 2 / mRMin); }
64+
double getFrontFaceMaxEta(int i);
6465
double getCrystalPhiMin();
6566
double getSamplingPhiMin();
6667
int getNModulesZ() { return mNModulesZ; }

‎Detectors/Upgrades/ALICE3/ECal/base/src/Geometry.cxx‎

Lines changed: 30 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -73,11 +73,18 @@ double Geometry::getSamplingPhiMin()
7373
return (superModuleDeltaPhi - samplingDeltaPhi * mNSamplingModulesPhi) / 2.;
7474
}
7575

76+
double Geometry::getFrontFaceMaxEta(int i)
77+
{
78+
double theta = std::atan(mRMin / getFrontFaceZatMinR(i));
79+
return -std::log(std::tan(theta / 2.));
80+
}
81+
7682
//==============================================================================
7783
void Geometry::fillFrontFaceCenterCoordinates()
7884
{
79-
if (mFrontFaceCenterR.size() > 0)
85+
if (mFrontFaceCenterR.size() > 0) {
8086
return;
87+
}
8188
mFrontFaceCenterTheta.resize(mNCrystalModulesZ + mNSamplingModulesZ);
8289
mFrontFaceZatMinR.resize(mNCrystalModulesZ + mNSamplingModulesZ);
8390
mFrontFaceCenterR.resize(mNCrystalModulesZ + mNSamplingModulesZ);
@@ -153,7 +160,7 @@ int Geometry::getCellID(int moduleId, int sectorId, bool isCrystal)
153160
if (sectorId % 2 == 0) { // sampling at positive eta
154161
cellID = sectorId / 2 * mNModulesZ + moduleId + mNSamplingModulesZ + mNCrystalModulesZ * 2;
155162
} else { // sampling at negative eta
156-
cellID = sectorId / 2 * mNModulesZ - moduleId + mNSamplingModulesZ;
163+
cellID = sectorId / 2 * mNModulesZ - moduleId + mNSamplingModulesZ - 1;
157164
}
158165
}
159166
return cellID;
@@ -206,13 +213,15 @@ void Geometry::detIdToGlobalPosition(int detId, double& x, double& y, double& z)
206213
{
207214
int chamber, sector, iphi, iz;
208215
detIdToRelIndex(detId, chamber, sector, iphi, iz);
216+
double r = 0;
209217
if (iz < mNSamplingModulesZ + mNCrystalModulesZ) {
210218
z = -mFrontFaceCenterZ[mNSamplingModulesZ + mNCrystalModulesZ - iz - 1];
219+
r = mFrontFaceCenterR[mNSamplingModulesZ + mNCrystalModulesZ - iz - 1];
211220
} else {
212-
z = +mFrontFaceCenterZ[iz % (mNSamplingModulesZ + mNCrystalModulesZ)];
221+
z = mFrontFaceCenterZ[iz % (mNSamplingModulesZ + mNCrystalModulesZ)];
222+
r = mFrontFaceCenterR[iz % (mNSamplingModulesZ + mNCrystalModulesZ)];
213223
}
214224
double phi = chamber == 1 ? mFrontFaceCenterCrystalPhi[iphi] : mFrontFaceCenterSamplingPhi[iphi];
215-
double r = mFrontFaceCenterR[iz % (mNSamplingModulesZ + mNCrystalModulesZ)];
216225
x = r * std::cos(phi);
217226
y = r * std::sin(phi);
218227
}
@@ -224,40 +233,45 @@ int Geometry::areNeighboursVertex(int detId1, int detId2) const
224233
int ch2, sector2, iphi2, iz2;
225234
detIdToRelIndex(detId1, ch1, sector1, iphi1, iz1);
226235
detIdToRelIndex(detId2, ch2, sector2, iphi2, iz2);
227-
if (sector1 != sector2 || ch1 != ch2)
236+
if (sector1 != sector2 || ch1 != ch2) {
228237
return 0;
229-
if (std::abs(iphi1 - iphi2) <= 1 && std::abs(iz1 - iz2) <= 1)
238+
}
239+
if (std::abs(iphi1 - iphi2) <= 1 && std::abs(iz1 - iz2) <= 1) {
230240
return 1;
241+
}
231242
return 0;
232243
}
233244

234245
//==============================================================================
235246
bool Geometry::isAtTheEdge(int cellId)
236247
{
237248
auto [row, col] = globalRowColFromIndex(cellId);
238-
if (col == 0)
249+
if (col == 0) {
239250
return 1;
240-
if (col == mNSamplingModulesZ)
251+
} else if (col == mNSamplingModulesZ) {
241252
return 1;
242-
if (col == mNSamplingModulesZ - 1)
253+
} else if (col == mNSamplingModulesZ - 1) {
243254
return 1;
244-
if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ)
255+
} else if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ) {
245256
return 1;
246-
if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ - 1)
257+
} else if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ - 1) {
247258
return 1;
248-
if (col == mNModulesZ - 1)
259+
} else if (col == mNModulesZ - 1) {
249260
return 1;
261+
}
250262
for (int m = 0; m <= mNSuperModules; m++) {
251263
if (isCrystal(cellId)) {
252-
if (row == m * mNCrystalModulesPhi)
264+
if (row == m * mNCrystalModulesPhi) {
253265
return 1;
254-
if (row == m * mNCrystalModulesPhi - 1)
266+
} else if (row == m * mNCrystalModulesPhi - 1) {
255267
return 1;
268+
}
256269
} else {
257-
if (row == m * mNSamplingModulesPhi)
270+
if (row == m * mNSamplingModulesPhi) {
258271
return 1;
259-
if (row == m * mNSamplingModulesPhi - 1)
272+
} else if (row == m * mNSamplingModulesPhi - 1) {
260273
return 1;
274+
}
261275
}
262276
}
263277
return 0;

‎Detectors/Upgrades/ALICE3/ECal/reconstruction/include/ECalReconstruction/Clusterizer.h‎

Lines changed: 25 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -48,25 +48,33 @@ class Clusterizer
4848
void setClusteringThreshold(double threshold) { mClusteringThreshold = threshold; }
4949
void setCrystalDigitThreshold(double threshold) { mCrystalDigitThreshold = threshold; }
5050
void setSamplingDigitThreshold(double threshold) { mSamplingDigitThreshold = threshold; }
51+
void setCrystalEnergyCorrectionPars(std::vector<double> pars) { mCrystalEnergyCorrectionPars = pars; }
52+
void setSamplingEnergyCorrectionPars(std::vector<double> pars) { mSamplingEnergyCorrectionPars = pars; }
53+
void setCrystalZCorrectionPars(std::vector<double> pars) { mCrystalZCorrectionPars = pars; }
54+
void setSamplingZCorrectionPars(std::vector<double> pars) { mSamplingZCorrectionPars = pars; }
5155

5256
private:
53-
std::vector<std::vector<int>> mDigitIndices; // 2D map of digit indices used for recursive cluster finding
54-
bool mUnfoldClusters = true; // to perform cluster unfolding
55-
double mCrystalDigitThreshold = 0.040; // minimal energy of crystal digit
56-
double mSamplingDigitThreshold = 0.100; // minimal energy of sampling digit
57-
double mClusteringThreshold = 0.050; // minimal energy of digit to start clustering (GeV)
58-
double mClusteringTimeGate = 1e9; // maximal time difference between digits to be accepted to clusters (in ns)
59-
int mNLMMax = 30; // maximal number of local maxima in unfolding
60-
double mLogWeight = 4.; // cutoff used in log. weight calculation
61-
double mUnfogingEAccuracy = 1.e-4; // accuracy of energy calculation in unfoding prosedure (GeV)
62-
double mUnfogingXZAccuracy = 1.e-2; // accuracy of position calculation in unfolding procedure (cm)
63-
int mNMaxIterations = 100; // maximal number of iterations in unfolding procedure
64-
double mLocalMaximumCut = 0.015; // minimal height of local maximum over neighbours
65-
bool mApplyCorrectionZ = 1; // z-correction
66-
bool mApplyCorrectionE = 1; // energy-correction
67-
TF1* fCrystalShowerShape; //! Crystal shower shape
68-
TF1* fSamplingShowerShape; //! Sampling shower shape
69-
TF1* fCrystalRMS; //! Crystal RMS
57+
std::vector<std::vector<int>> mDigitIndices; // 2D map of digit indices used for recursive cluster finding
58+
bool mUnfoldClusters = true; // to perform cluster unfolding
59+
double mCrystalDigitThreshold = 0.040; // minimal energy of crystal digit
60+
double mSamplingDigitThreshold = 0.100; // minimal energy of sampling digit
61+
double mClusteringThreshold = 0.050; // minimal energy of digit to start clustering (GeV)
62+
double mClusteringTimeGate = 1e9; // maximal time difference between digits to be accepted to clusters (in ns)
63+
int mNLMMax = 30; // maximal number of local maxima in unfolding
64+
double mLogWeight = 4.; // cutoff used in log. weight calculation
65+
double mUnfogingEAccuracy = 1.e-4; // accuracy of energy calculation in unfoding prosedure (GeV)
66+
double mUnfogingXZAccuracy = 1.e-2; // accuracy of position calculation in unfolding procedure (cm)
67+
int mNMaxIterations = 100; // maximal number of iterations in unfolding procedure
68+
double mLocalMaximumCut = 0.015; // minimal height of local maximum over neighbours
69+
bool mApplyCorrectionZ = 1; // apply z-correction
70+
bool mApplyCorrectionE = 1; // apply energy-correction
71+
TF1* fCrystalShowerShape; //! Crystal shower shape
72+
TF1* fSamplingShowerShape; //! Sampling shower shape
73+
TF1* fCrystalRMS; //! Crystal RMS
74+
std::vector<double> mCrystalEnergyCorrectionPars; // crystal energy-correction parameters
75+
std::vector<double> mSamplingEnergyCorrectionPars; // sampling energy-correction parameters
76+
std::vector<double> mCrystalZCorrectionPars; // crystal z-correction parameters
77+
std::vector<double> mSamplingZCorrectionPars; // sampling z-correction parameters
7078
};
7179

7280
} // namespace ecal

‎Detectors/Upgrades/ALICE3/ECal/reconstruction/src/Clusterizer.cxx‎

Lines changed: 71 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,45 @@ Clusterizer::Clusterizer(bool applyCorrectionZ, bool applyCorrectionE)
3131
mDigitIndices.resize(geo.getNrows(), std::vector<int>(geo.getNcols(), -1));
3232
mApplyCorrectionZ = applyCorrectionZ;
3333
mApplyCorrectionE = applyCorrectionE;
34+
35+
mCrystalEnergyCorrectionPars.reserve(6);
36+
mCrystalEnergyCorrectionPars[0] = 0.00444;
37+
mCrystalEnergyCorrectionPars[1] = -1.322;
38+
mCrystalEnergyCorrectionPars[2] = 1.021;
39+
mCrystalEnergyCorrectionPars[3] = 0.0018;
40+
mCrystalEnergyCorrectionPars[4] = 0.;
41+
mCrystalEnergyCorrectionPars[5] = 0.;
42+
43+
mSamplingEnergyCorrectionPars.reserve(6);
44+
mSamplingEnergyCorrectionPars[0] = 0.0033;
45+
mSamplingEnergyCorrectionPars[1] = -2.09;
46+
mSamplingEnergyCorrectionPars[2] = 1.007;
47+
mSamplingEnergyCorrectionPars[3] = 0.0667;
48+
mSamplingEnergyCorrectionPars[4] = -0.108;
49+
mSamplingEnergyCorrectionPars[5] = 0.0566;
50+
51+
mCrystalZCorrectionPars.reserve(9);
52+
mCrystalZCorrectionPars[0] = -0.005187;
53+
mCrystalZCorrectionPars[1] = 0.7301;
54+
mCrystalZCorrectionPars[2] = -0.7382;
55+
mCrystalZCorrectionPars[3] = 0.;
56+
mCrystalZCorrectionPars[4] = 0.;
57+
mCrystalZCorrectionPars[5] = 0.;
58+
mCrystalZCorrectionPars[6] = 0.;
59+
mCrystalZCorrectionPars[7] = 0.;
60+
mCrystalZCorrectionPars[8] = 0.;
61+
62+
mSamplingZCorrectionPars.reserve(9);
63+
mSamplingZCorrectionPars[0] = -2.137;
64+
mSamplingZCorrectionPars[1] = 6.400;
65+
mSamplingZCorrectionPars[2] = -3.342;
66+
mSamplingZCorrectionPars[3] = -0.1364;
67+
mSamplingZCorrectionPars[4] = 0.4019;
68+
mSamplingZCorrectionPars[5] = -0.1969;
69+
mSamplingZCorrectionPars[6] = 0.008223;
70+
mSamplingZCorrectionPars[7] = -0.02425;
71+
mSamplingZCorrectionPars[8] = 0.01190;
72+
3473
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);
3574
double pc[12];
3675
pc[0] = 1. / 13.15;
@@ -94,22 +133,24 @@ void Clusterizer::findClusters(const gsl::span<const Digit>& digits, std::vector
94133
void Clusterizer::addDigitToCluster(Cluster& cluster, int row, int col, const gsl::span<const Digit>& digits)
95134
{
96135
auto& geo = Geometry::instance();
97-
if (row < 0 || row >= geo.getNrows() || col < 0 || col >= geo.getNcols())
136+
if (row < 0 || row >= geo.getNrows() || col < 0 || col >= geo.getNcols()) {
98137
return;
138+
}
99139
int digitIndex = mDigitIndices[row][col];
100140
LOGP(debug, " checking row={} and col={} digitIndex={}", row, col, digitIndex);
101-
if (digitIndex < 0)
141+
if (digitIndex < 0) {
102142
return;
103-
143+
}
104144
const Digit& digit = digits[digitIndex];
105145
if (cluster.getMultiplicity() > 0) {
106146
// check if new digit is in the same chamber and sector
107147
const Digit& digit2 = digits[cluster.getDigitIndex(0)];
108148
auto [sector1, ch1] = geo.getSectorChamber(digit.getTower());
109149
auto [sector2, ch2] = geo.getSectorChamber(digit2.getTower());
110150
LOGP(debug, " checking if sector and chamber are the same: ({},{}) ({},{})", sector1, ch1, sector2, ch2);
111-
if (sector1 != sector2 || ch1 != ch2)
151+
if (sector1 != sector2 || ch1 != ch2) {
112152
return;
153+
}
113154
}
114155

115156
mDigitIndices[row][col] = -1;
@@ -140,11 +181,13 @@ void Clusterizer::makeClusters(const gsl::span<const Digit>& digits, std::vector
140181
auto [row, col] = geo.globalRowColFromIndex(digit.getTower());
141182
bool isCrystal = geo.isCrystal(digit.getTower());
142183
if (isCrystal) {
143-
if (digit.getEnergy() < mCrystalDigitThreshold)
184+
if (digit.getEnergy() < mCrystalDigitThreshold) {
144185
continue;
186+
}
145187
} else {
146-
if (digit.getEnergy() < mSamplingDigitThreshold)
188+
if (digit.getEnergy() < mSamplingDigitThreshold) {
147189
continue;
190+
}
148191
}
149192
mDigitIndices[row][col] = i;
150193
}
@@ -153,10 +196,12 @@ void Clusterizer::makeClusters(const gsl::span<const Digit>& digits, std::vector
153196
for (int i = 0; i < nDigits; i++) {
154197
const Digit& digitSeed = digits[i];
155198
auto [row, col] = geo.globalRowColFromIndex(digitSeed.getTower());
156-
if (mDigitIndices[row][col] < 0)
199+
if (mDigitIndices[row][col] < 0) {
157200
continue; // digit was already added in one of the clusters
158-
if (digitSeed.getEnergy() < mClusteringThreshold)
201+
}
202+
if (digitSeed.getEnergy() < mClusteringThreshold) {
159203
continue;
204+
}
160205
LOGP(debug, " starting new cluster at row={} and col={}", row, col);
161206
auto& cluster = clusters.emplace_back();
162207
addDigitToCluster(cluster, row, col, digits);
@@ -343,8 +388,9 @@ void Clusterizer::evalClusters(std::vector<Cluster>& clusters)
343388
double xi, yi, zi;
344389
geo.detIdToGlobalPosition(towerId, xi, yi, zi);
345390
double r = std::sqrt((x - xi) * (x - xi) + (y - yi) * (y - yi) + (z - zi) * (z - zi));
346-
if (r > 2.2)
391+
if (r > 2.2) {
347392
continue;
393+
}
348394
double frac = fCrystalShowerShape->Eval(r);
349395
double rms = fCrystalRMS->Eval(r);
350396
chi2 += std::pow((energy / ee - frac) / rms, 2.);
@@ -354,38 +400,30 @@ void Clusterizer::evalClusters(std::vector<Cluster>& clusters)
354400

355401
// correct cluster energy and z position
356402
float eta = std::abs(cluster.getEta());
357-
float eCor = 1;
358-
float zCor = 0;
359403
bool isCrystal = geo.isCrystal(cluster.getDigitTowerId(0));
360-
if (isCrystal) {
361-
eCor = 0.00444 * std::pow(ee, -1.322) + (1.021 + 0.0018 * eta);
362-
if (mApplyCorrectionE)
363-
ee *= eCor;
364-
if (mApplyCorrectionZ)
365-
zCor = (-0.00518682 + 0.730052 * eta - 0.73817 * eta * eta);
366-
} else {
367-
eCor = 0.0033 * std::pow(ee, -2.09) + (1.007 + 0.0667 * eta - 0.108 * eta * eta + 0.0566 * eta * eta * eta);
368-
if (mApplyCorrectionE)
369-
ee *= eCor;
370-
if (mApplyCorrectionZ)
371-
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;
404+
if (mApplyCorrectionE) {
405+
std::vector<double>& pe = isCrystal ? mCrystalEnergyCorrectionPars : mSamplingEnergyCorrectionPars;
406+
ee *= pe[0] * std::pow(ee, pe[1]) + pe[2] + pe[3] * eta + pe[4] * eta * eta + pe[5] * eta * eta * eta;
407+
cluster.setE(ee);
408+
}
409+
if (mApplyCorrectionZ) {
410+
std::vector<double>& pz = isCrystal ? mCrystalZCorrectionPars : mSamplingZCorrectionPars;
411+
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;
412+
cluster.setZ(z > 0 ? z - zCor : z + zCor);
372413
}
373-
374-
cluster.setE(ee);
375-
cluster.setZ(cluster.getZ() - zCor);
376414

377415
// check if cluster is at the edge of detector module
378416
bool isEdge = 0;
379417
for (size_t i = 0; i < cluster.getMultiplicity(); i++) {
380418
int towerId = cluster.getDigitTowerId(i);
381-
if (!geo.isAtTheEdge(towerId))
382-
continue;
383-
isEdge = 1;
384-
break;
419+
if (geo.isAtTheEdge(towerId)) {
420+
isEdge = 1;
421+
break;
422+
}
385423
}
386424
cluster.setEdgeFlag(isEdge);
387425

388-
LOGF(debug, "Cluster coordinates: (%6.2f,%6.2f,%6.2f), eCor=%6.2f zCor=%6.2f", cluster.getX(), cluster.getY(), cluster.getZ(), eCor, zCor);
426+
LOGF(debug, "Cluster coordinates: (%6.2f,%6.2f,%6.2f)", cluster.getX(), cluster.getY(), cluster.getZ());
389427
}
390428
}
391429

@@ -403,8 +441,9 @@ int Clusterizer::getNumberOfLocalMax(Cluster& clu, int* maxAt, float* maxAtEnerg
403441
for (int i = 0; i < n; i++) {
404442
isLocalMax[i] = false;
405443
float en1 = clu.getDigitEnergy(i);
406-
if (en1 > mClusteringThreshold)
444+
if (en1 > mClusteringThreshold) {
407445
isLocalMax[i] = true;
446+
}
408447
}
409448

410449
for (int i = 0; i < n; i++) {

‎Detectors/Upgrades/ALICE3/ECal/simulation/src/Detector.cxx‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -358,9 +358,9 @@ bool Detector::ProcessHits(FairVolume* vol)
358358
return false;
359359
}
360360

361-
if (isCrystal)
361+
if (isCrystal) {
362362
LOGP(debug, "Processing crystal {}", volName.Data());
363-
else {
363+
} else {
364364
eloss *= mSamplingFactorTransportModel;
365365
LOGP(debug, "Processing scintillator {}", volName.Data());
366366
}

‎Detectors/Upgrades/ALICE3/ECal/simulation/src/Digitizer.cxx‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -57,16 +57,18 @@ void Digitizer::processHits(const std::vector<Hit>* hits, std::vector<Digit>& di
5757
bool isCrystal = geo.isCrystal(cellID);
5858
if (isCrystal) { // crystal
5959
double elossSmearedNpe = gRandom->Poisson(eloss * mCrystalPePerGeV) / mCrystalPePerGeV;
60-
if (mSmearCrystal)
60+
if (mSmearCrystal) {
6161
elossSmeared = elossSmearedNpe * gRandom->Gaus(1, 0.007); // light attenuation in crystals
62-
} else { // sampling
62+
}
63+
} else { // sampling
6364
elossSmeared *= mSamplingFraction;
6465
}
6566

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

7274
// Adding MC info

0 commit comments

Comments
 (0)