Project
Loading...
Searching...
No Matches
Clusterizer.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
16
17#include <TMath.h>
18#include <TF1.h>
19#include <fairlogger/Logger.h>
21#include <ECalBase/Geometry.h>
24
25using namespace o2::ecal;
26
27//==============================================================================
28Clusterizer::Clusterizer(bool applyCorrectionZ, bool applyCorrectionE)
29{
30 auto& geo = Geometry::instance();
31 mDigitIndices.resize(geo.getNrows(), std::vector<int>(geo.getNcols(), -1));
32 mApplyCorrectionZ = applyCorrectionZ;
33 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
73 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);
74 double pc[12];
75 pc[0] = 1. / 13.15;
76 pc[1] = 2.2;
77 pc[2] = 5;
78 pc[3] = 4.38969;
79 pc[4] = -5.15975;
80 pc[5] = 1.18978;
81 pc[6] = 1.48726;
82 pc[7] = -1.54621;
83 pc[8] = 0.0814617;
84 pc[9] = 0.0369055;
85 pc[10] = -0.174372;
86 pc[11] = -0.0455978;
87
88 fCrystalShowerShape->SetParameters(pc);
89
90 fSamplingShowerShape = new TF1("fSampling", "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);
91 double ps[12];
92 ps[0] = 1 / 35.6;
93 ps[1] = 3.2;
94 ps[2] = 6;
95 ps[3] = 3.06543;
96 ps[4] = -2.23235;
97 ps[5] = 0.325344;
98 ps[6] = 6.0733;
99 ps[7] = -1.62713;
100 ps[8] = 0.0965569;
101 ps[9] = 0.0765706;
102 ps[10] = -0.217398;
103 ps[11] = -0.0204646;
104 fSamplingShowerShape->SetParameters(ps);
105
106 fCrystalRMS = new TF1("fCrystalRMS", "[0]*x*exp([1]*x+[2]*x*x+[3]*x*x*x)", 0, 2.2);
107 double p[4];
108 p[0] = 1.39814;
109 p[1] = -6.05426;
110 p[2] = 6.26678;
111 p[3] = -1.97092;
112 fCrystalRMS->SetParameters(p);
113}
114
115//==============================================================================
116void Clusterizer::findClusters(const gsl::span<const Digit>& digits, std::vector<Cluster>& foundClusters, std::vector<Cluster>& unfoldedClusters)
117{
118 foundClusters.clear();
119 unfoldedClusters.clear();
120
121 // Collect list of clusters
122 makeClusters(digits, foundClusters);
123
124 // Split clusters with several local maxima if necessary
125 makeUnfoldings(foundClusters, unfoldedClusters);
126
127 // Evaluate cluster position, dispersion etc.
128 evalClusters(foundClusters);
129 evalClusters(unfoldedClusters);
130}
131
132//==============================================================================
133void Clusterizer::addDigitToCluster(Cluster& cluster, int row, int col, const gsl::span<const Digit>& digits)
134{
135 auto& geo = Geometry::instance();
136 if (row < 0 || row >= geo.getNrows() || col < 0 || col >= geo.getNcols()) {
137 return;
138 }
139 int digitIndex = mDigitIndices[row][col];
140 LOGP(debug, " checking row={} and col={} digitIndex={}", row, col, digitIndex);
141 if (digitIndex < 0) {
142 return;
143 }
144 const Digit& digit = digits[digitIndex];
145 if (cluster.getMultiplicity() > 0) {
146 // check if new digit is in the same chamber and sector
147 const Digit& digit2 = digits[cluster.getDigitIndex(0)];
148 auto [sector1, ch1] = geo.getSectorChamber(digit.getTower());
149 auto [sector2, ch2] = geo.getSectorChamber(digit2.getTower());
150 LOGP(debug, " checking if sector and chamber are the same: ({},{}) ({},{})", sector1, ch1, sector2, ch2);
151 if (sector1 != sector2 || ch1 != ch2) {
152 return;
153 }
154 }
155
156 mDigitIndices[row][col] = -1;
157 cluster.addDigit(digitIndex, digit.getTower(), digit.getEnergy());
158 LOGP(debug, " adding new digit at row={} and col={}", row, col);
159 addDigitToCluster(cluster, row - 1, col, digits);
160 addDigitToCluster(cluster, row + 1, col, digits);
161 addDigitToCluster(cluster, row, col - 1, digits);
162 addDigitToCluster(cluster, row, col + 1, digits);
163}
164
165//==============================================================================
166void Clusterizer::makeClusters(const gsl::span<const Digit>& digits, std::vector<Cluster>& clusters)
167{
168 // Combine digits into cluster
169
170 int nDigits = digits.size();
171
172 // reset mDigitIndices
173 for (auto& rows : mDigitIndices) {
174 rows.assign(rows.size(), -1);
175 }
176
177 // fill mDigitIndices
178 auto& geo = Geometry::instance();
179 for (int i = 0; i < nDigits; i++) {
180 const Digit& digit = digits[i];
181 auto [row, col] = geo.globalRowColFromIndex(digit.getTower());
182 bool isCrystal = geo.isCrystal(digit.getTower());
183 if (isCrystal) {
184 if (digit.getEnergy() < mCrystalDigitThreshold) {
185 continue;
186 }
187 } else {
188 if (digit.getEnergy() < mSamplingDigitThreshold) {
189 continue;
190 }
191 }
192 mDigitIndices[row][col] = i;
193 }
194
195 // add digit seeds to clusters and recursively add neighbours
196 for (int i = 0; i < nDigits; i++) {
197 const Digit& digitSeed = digits[i];
198 auto [row, col] = geo.globalRowColFromIndex(digitSeed.getTower());
199 if (mDigitIndices[row][col] < 0) {
200 continue; // digit was already added in one of the clusters
201 }
202 if (digitSeed.getEnergy() < mClusteringThreshold) {
203 continue;
204 }
205 LOGP(debug, " starting new cluster at row={} and col={}", row, col);
206 auto& cluster = clusters.emplace_back();
207 addDigitToCluster(cluster, row, col, digits);
208 }
209
210 LOGP(debug, "made {} clusters from {} digits", clusters.size(), nDigits);
211}
212
213//==============================================================================
214void Clusterizer::makeUnfoldings(std::vector<Cluster>& foundClusters, std::vector<Cluster>& unfoldedClusters)
215{
216 // Split cluster if several local maxima are found
217 if (!mUnfoldClusters) {
218 return;
219 }
220
221 int* maxAt = new int[mNLMMax];
222 float* maxAtEnergy = new float[mNLMMax];
223
224 for (auto& clu : foundClusters) {
225 int nMax = getNumberOfLocalMax(clu, maxAt, maxAtEnergy);
226 if (nMax > 1) {
227 unfoldOneCluster(&clu, nMax, maxAt, maxAtEnergy, unfoldedClusters);
228 } else {
229 clu.setNLM(1);
230 unfoldedClusters.emplace_back(clu);
231 }
232 }
233 delete[] maxAt;
234 delete[] maxAtEnergy;
235}
236
237//==============================================================================
238void Clusterizer::unfoldOneCluster(Cluster* iniClu, int nMax, int* digitId, float* maxAtEnergy, std::vector<Cluster>& unfoldedClusters)
239{
240 // Based on MpdEmcClusterizerKI::UnfoldOneCluster by D. Peresunko
241 // Performs the unfolding of a cluster with nMax overlapping showers
242 // Parameters: iniClu cluster to be unfolded
243 // nMax number of local maxima found (this is the number of new clusters)
244 // digitId: index of digits, corresponding to local maxima
245 // maxAtEnergy: energies of digits, corresponding to local maxima
246
247 // Take initial cluster and calculate local coordinates of digits
248 // To avoid multiple re-calculation of same parameters
249 int mult = iniClu->getMultiplicity();
250 std::vector<double> x(mult);
251 std::vector<double> y(mult);
252 std::vector<double> z(mult);
253 std::vector<double> e(mult);
254 std::vector<std::vector<double>> eInClusters(mult, std::vector<double>(nMax));
255
256 auto& geo = Geometry::instance();
257 bool isCrystal = geo.isCrystal(iniClu->getDigitTowerId(0));
258
259 for (int idig = 0; idig < mult; idig++) {
260 e[idig] = iniClu->getDigitEnergy(idig);
261 geo.detIdToGlobalPosition(iniClu->getDigitTowerId(idig), x[idig], y[idig], z[idig]);
262 }
263
264 // Coordinates of centers of clusters
265 std::vector<double> xMax(nMax);
266 std::vector<double> yMax(nMax);
267 std::vector<double> zMax(nMax);
268 std::vector<double> eMax(nMax);
269
270 for (int iclu = 0; iclu < nMax; iclu++) {
271 xMax[iclu] = x[digitId[iclu]];
272 yMax[iclu] = y[digitId[iclu]];
273 zMax[iclu] = z[digitId[iclu]];
274 eMax[iclu] = e[digitId[iclu]];
275 }
276
277 std::vector<double> prop(nMax); // proportion of clusters in the current digit
278
279 // Try to decompose cluster to contributions
280 int nIterations = 0;
281 bool insuficientAccuracy = true;
282
283 while (insuficientAccuracy && nIterations < mNMaxIterations) {
284 // Loop over all digits of parent cluster and split their energies between daughter clusters
285 // according to shower shape
286 for (int idig = 0; idig < mult; idig++) {
287 double eEstimated = 0;
288 for (int iclu = 0; iclu < nMax; iclu++) {
289 prop[iclu] = eMax[iclu] * showerShape(std::sqrt((x[idig] - xMax[iclu]) * (x[idig] - xMax[iclu]) +
290 (y[idig] - yMax[iclu]) * (y[idig] - yMax[iclu])),
291 z[idig] - zMax[iclu], isCrystal);
292 eEstimated += prop[iclu];
293 }
294 if (eEstimated == 0.) { // numerical accuracy
295 continue;
296 }
297 // Split energy of digit according to contributions
298 for (int iclu = 0; iclu < nMax; iclu++) {
299 eInClusters[idig][iclu] = e[idig] * prop[iclu] / eEstimated;
300 }
301 }
302 // Recalculate parameters of clusters and check relative variation of energy and absolute of position
303 insuficientAccuracy = false; // will be true if at least one parameter changed too much
304 for (int iclu = 0; iclu < nMax; iclu++) {
305 double oldX = xMax[iclu];
306 double oldY = yMax[iclu];
307 double oldZ = zMax[iclu];
308 double oldE = eMax[iclu];
309 // new energy, need for weight
310 eMax[iclu] = 0;
311 for (int idig = 0; idig < mult; idig++) {
312 eMax[iclu] += eInClusters[idig][iclu];
313 }
314 xMax[iclu] = 0;
315 yMax[iclu] = 0;
316 zMax[iclu] = 0;
317 double wtot = 0.;
318 for (int idig = 0; idig < mult; idig++) {
319 double w = std::max(std::log(eInClusters[idig][iclu] / eMax[iclu]) + mLogWeight, 0.);
320 xMax[iclu] += x[idig] * w;
321 yMax[iclu] += y[idig] * w;
322 zMax[iclu] += z[idig] * w;
323 wtot += w;
324 }
325 if (wtot > 0.) {
326 xMax[iclu] /= wtot;
327 yMax[iclu] /= wtot;
328 zMax[iclu] /= wtot;
329 }
330 // Compare variation of parameters
331 insuficientAccuracy += (std::abs(eMax[iclu] - oldE) > mUnfogingEAccuracy);
332 insuficientAccuracy += (std::abs(xMax[iclu] - oldX) > mUnfogingXZAccuracy);
333 insuficientAccuracy += (std::abs(yMax[iclu] - oldY) > mUnfogingXZAccuracy);
334 insuficientAccuracy += (std::abs(zMax[iclu] - oldZ) > mUnfogingXZAccuracy);
335 }
336 nIterations++;
337 }
338
339 // Iterations finished, add new clusters
340 for (int iclu = 0; iclu < nMax; iclu++) {
341 auto& clu = unfoldedClusters.emplace_back();
342 clu.setNLM(nMax);
343 for (int idig = 0; idig < mult; idig++) {
344 int jdigit = iniClu->getDigitIndex(idig);
345 int towerId = iniClu->getDigitTowerId(idig);
346 clu.addDigit(jdigit, towerId, eInClusters[idig][iclu]);
347 }
348 }
349}
350
351//==============================================================================
352void Clusterizer::evalClusters(std::vector<Cluster>& clusters)
353{
354 auto& geo = Geometry::instance();
355 for (auto& cluster : clusters) {
356 double x = 0;
357 double y = 0;
358 double z = 0;
359 double wtot = 0;
360 double etot = cluster.getE();
361 for (size_t i = 0; i < cluster.getMultiplicity(); i++) {
362 float energy = cluster.getDigitEnergy(i);
363 int towerId = cluster.getDigitTowerId(i);
364 double xi, yi, zi;
365 geo.detIdToGlobalPosition(towerId, xi, yi, zi);
366 double w = std::max(0., mLogWeight + std::log(energy / etot));
367 x += w * xi;
368 y += w * yi;
369 z += w * zi;
370 wtot += w;
371 }
372 if (wtot != 0) {
373 x /= wtot;
374 y /= wtot;
375 z /= wtot;
376 }
377 cluster.setX(x);
378 cluster.setY(y);
379 cluster.setZ(z);
380
381 // cluster shape
382 float chi2 = 0;
383 int ndf = 0;
384 float ee = cluster.getE();
385 for (size_t i = 0; i < cluster.getMultiplicity(); i++) {
386 float energy = cluster.getDigitEnergy(i);
387 int towerId = cluster.getDigitTowerId(i);
388 double xi, yi, zi;
389 geo.detIdToGlobalPosition(towerId, xi, yi, zi);
390 double r = std::sqrt((x - xi) * (x - xi) + (y - yi) * (y - yi) + (z - zi) * (z - zi));
391 if (r > 2.2) {
392 continue;
393 }
394 double frac = fCrystalShowerShape->Eval(r);
395 double rms = fCrystalRMS->Eval(r);
396 chi2 += std::pow((energy / ee - frac) / rms, 2.);
397 ndf++;
398 }
399 cluster.setChi2(chi2 / ndf);
400
401 // correct cluster energy and z position
402 float eta = std::abs(cluster.getEta());
403 bool isCrystal = geo.isCrystal(cluster.getDigitTowerId(0));
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);
413 }
414
415 // check if cluster is at the edge of detector module
416 bool isEdge = 0;
417 for (size_t i = 0; i < cluster.getMultiplicity(); i++) {
418 int towerId = cluster.getDigitTowerId(i);
419 if (geo.isAtTheEdge(towerId)) {
420 isEdge = 1;
421 break;
422 }
423 }
424 cluster.setEdgeFlag(isEdge);
425
426 LOGF(debug, "Cluster coordinates: (%6.2f,%6.2f,%6.2f)", cluster.getX(), cluster.getY(), cluster.getZ());
427 }
428}
429
430//==============================================================================
431int Clusterizer::getNumberOfLocalMax(Cluster& clu, int* maxAt, float* maxAtEnergy)
432{
433 // Based on MpdEmcClusterizerKI::GetNumberOfLocalMax by D. Peresunko
434 // Calculates the number of local maxima in the cluster using LocalMaxCut as the minimum
435 // energy difference between maximum and surrounding digits
436 auto& geo = Geometry::instance();
437 int n = clu.getMultiplicity();
438 bool isCrystal = geo.isCrystal(clu.getDigitTowerId(0));
439 bool* isLocalMax = new bool[n];
440
441 for (int i = 0; i < n; i++) {
442 isLocalMax[i] = false;
443 float en1 = clu.getDigitEnergy(i);
444 if (en1 > mClusteringThreshold) {
445 isLocalMax[i] = true;
446 }
447 }
448
449 for (int i = 0; i < n; i++) {
450 int detId1 = clu.getDigitTowerId(i);
451 float en1 = clu.getDigitEnergy(i);
452 for (int j = i + 1; j < n; j++) {
453 int detId2 = clu.getDigitTowerId(j);
454 float en2 = clu.getDigitEnergy(j);
455 if (geo.areNeighboursVertex(detId1, detId2) == 1) {
456 if (en1 > en2) {
457 isLocalMax[j] = false;
458 // but may be digit too is not local max ?
459 if (en2 > en1 - mLocalMaximumCut) {
460 isLocalMax[i] = false;
461 }
462 } else {
463 isLocalMax[i] = false;
464 // but may be digitN is not local max too?
465 if (en1 > en2 - mLocalMaximumCut) {
466 isLocalMax[j] = false;
467 }
468 }
469 } // if neighbours
470 } // digit j
471 } // digit i
472
473 int iDigitN = 0;
474 for (int i = 0; i < n; i++) {
475 if (isLocalMax[i]) {
476 maxAt[iDigitN] = i;
477 maxAtEnergy[iDigitN] = clu.getDigitEnergy(i);
478 iDigitN++;
479 if (iDigitN >= mNLMMax) { // Note that size of output arrays is limited:
480 LOGP(error, "Too many local maxima, cluster multiplicity {} region={}", n, isCrystal ? "crystal" : "sampling");
481 return 0;
482 }
483 }
484 }
485 delete[] isLocalMax;
486 return iDigitN;
487}
488
489//==============================================================================
490double Clusterizer::showerShape(double dx, double dz, bool isCrystal)
491{
492 double x = std::sqrt(dx * dx + dz * dz);
493 return isCrystal ? fCrystalShowerShape->Eval(x) : fSamplingShowerShape->Eval(x);
494}
Definition of ECal cluster class.
Definition of ECal digit class.
int32_t i
float float float & zMax
float & yMax
uint32_t j
Definition RawData.h:0
uint32_t col
Definition RawData.h:4
Geometry helper class.
Class for cluster finding and unfolding.
std::ostringstream debug
float getDigitEnergy(int i) const
Definition Cluster.h:51
void addDigit(int digitIndex, int towerId, double energy)
Definition Cluster.cxx:27
int getMultiplicity() const
Definition Cluster.h:48
int getDigitTowerId(int i) const
Definition Cluster.h:50
int getDigitIndex(int i) const
Definition Cluster.h:49
void setNLM(int nMax)
Definition Cluster.h:37
Clusterizer(bool applyCorrectionZ=1, bool applyCorrectionE=1)
void makeUnfoldings(std::vector< Cluster > &foundClusters, std::vector< Cluster > &unfoldedClusters)
void unfoldOneCluster(Cluster *iniClu, int nMax, int *digitId, float *maxAtEnergy, std::vector< Cluster > &unfoldedClusters)
void findClusters(const gsl::span< const Digit > &digits, std::vector< Cluster > &foundClusters, std::vector< Cluster > &unfoldedClusters)
void makeClusters(const gsl::span< const Digit > &digits, std::vector< Cluster > &clusters)
double showerShape(double dx, double dz, bool isCrystal)
void addDigitToCluster(Cluster &cluster, int row, int col, const gsl::span< const Digit > &digits)
void evalClusters(std::vector< Cluster > &clusters)
int getNumberOfLocalMax(Cluster &clu, int *maxAt, float *maxAtEnergy)
int getTower() const
Definition Digit.h:41
double getEnergy() const
Definition Digit.h:43
static Geometry & instance()
Definition Geometry.h:31
GLdouble n
Definition glcorearb.h:1982
GLint GLenum GLint x
Definition glcorearb.h:403
GLint y
Definition glcorearb.h:270
GLboolean r
Definition glcorearb.h:1233
GLubyte GLubyte GLubyte GLubyte w
Definition glcorearb.h:852
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
Cluster clu
std::vector< Cluster > clusters
std::vector< Digit > digits
std::vector< int > row
std::vector< ReadoutWindowData > rows