Project
Loading...
Searching...
No Matches
ClusterFactory.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
13#include <array>
14#include <gsl/span>
15#include "Rtypes.h"
23#include "EMCALBase/Geometry.h"
24#include "MathUtils/Cartesian.h"
25
27
28using namespace o2::emcal;
29
30template <class InputType>
31ClusterFactory<InputType>::ClusterFactory(gsl::span<const o2::emcal::Cluster> clustersContainer, gsl::span<const InputType> inputsContainer, gsl::span<const int> cellsIndices)
32{
33 setContainer(clustersContainer, inputsContainer, cellsIndices);
34}
35
36template <class InputType>
38{
39 mClustersContainer = gsl::span<const o2::emcal::Cluster>();
40 mInputsContainer = gsl::span<const InputType>();
41 mCellsIndices = gsl::span<int>();
42 mLookUpInit = false;
43 mCellLabelContainer = gsl::span<const o2::emcal::CellLabel>();
44}
45
48//____________________________________________________________________________
49template <class InputType>
51{
52 if (clusterIndex >= mClustersContainer.size()) {
53 throw ClusterRangeException(clusterIndex, mClustersContainer.size());
54 }
55 if (!mGeomPtr) {
57 }
58
59 o2::emcal::AnalysisCluster clusterAnalysis;
60 clusterAnalysis.setID(clusterIndex);
61
62 int firstCellIndex = mClustersContainer[clusterIndex].getCellIndexFirst();
63 int nCells = mClustersContainer[clusterIndex].getNCells();
64
65 gsl::span<const int> inputsIndices = gsl::span<const int>(&mCellsIndices[firstCellIndex], nCells);
66
67 // First calculate the index of input with maximum amplitude and get
68 // the supermodule number where it sits.
69
70 auto [inputIndMax, inputEnergyMax, cellAmp, shared] = getMaximalEnergyIndex(inputsIndices);
71
72 short towerId = mInputsContainer[inputIndMax].getTower();
73
74 float exoticTime = mInputsContainer[inputIndMax].getTimeStamp();
75
76 float fCross = 0.;
77
78 try {
79 clusterAnalysis.setIsExotic(isExoticCell(towerId, inputEnergyMax, exoticTime, fCross));
80 clusterAnalysis.setFCross(fCross);
81 } catch (UninitLookUpTableException& e) {
82 LOG(error) << e.what();
83 }
84
85 clusterAnalysis.setIndMaxInput(inputIndMax);
86
87 clusterAnalysis.setE(cellAmp);
88
89 mSuperModuleNumber = mGeomPtr->GetSuperModuleNumber(towerId);
90
91 clusterAnalysis.setNCells(inputsIndices.size());
92
93 std::vector<unsigned short> cellsIdices;
94
95 bool addClusterLabels = ((clusterLabel != nullptr) && (mCellLabelContainer.size() > 0));
96 for (auto cellIndex : inputsIndices) {
97 cellsIdices.push_back(cellIndex);
98 if (addClusterLabels) {
99 for (size_t iLabel = 0; iLabel < mCellLabelContainer[cellIndex].GetLabelSize(); iLabel++) {
100 if (mCellLabelContainer[cellIndex].GetAmplitudeFraction(iLabel) <= 0.f) {
101 continue; // skip 0 entries
102 }
103 clusterLabel->addValue(mCellLabelContainer[cellIndex].GetLabel(iLabel),
104 mCellLabelContainer[cellIndex].GetAmplitudeFraction(iLabel) * mInputsContainer[cellIndex].getEnergy());
105 }
106 }
107 }
108 if (addClusterLabels) {
109 clusterLabel->orderLabels();
110 clusterLabel->normalize(cellAmp);
111 }
112
113 clusterAnalysis.setCellsIndices(cellsIdices);
114
115 // evaluate global and local position
116 evalGlobalPosition(inputsIndices, clusterAnalysis);
117 evalLocalPosition(inputsIndices, clusterAnalysis);
118
119 // evaluate shower parameters
120 evalElipsAxis(inputsIndices, clusterAnalysis);
121 evalDispersion(inputsIndices, clusterAnalysis);
122
123 evalCoreEnergy(inputsIndices, clusterAnalysis);
124 evalTime(inputsIndices, clusterAnalysis);
125
126 // TODO to be added at a later stage
127 // evalPrimaries(inputsIndices, clusterAnalysis);
128 // evalParents(inputsIndices, clusterAnalysis);
129
130 // TODO to be added at a later stage
131 // Called last because it sets the global position of the cluster?
132 // Do not call it when recalculating clusters out of standard reconstruction
133 // if (!mJustCluster)
134 // evalLocal2TrackingCSTransform();
135
136 return clusterAnalysis;
137}
138
142//____________________________________________________________________________
143template <class InputType>
144void ClusterFactory<InputType>::evalDispersion(gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
145{
146 double d = 0., wtot = 0.;
147 int nstat = 0;
148
149 // Calculates the dispersion in cell units
150 double etaMean = 0.0, phiMean = 0.0;
151
152 // Calculate mean values
153 for (auto iInput : inputsIndices) {
154
155 if (clusterAnalysis.E() > 0 && mInputsContainer[iInput].getEnergy() > 0) {
156 auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
157 auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta);
158
159 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
160 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
161 if (mSharedCluster && nSupMod % 2) {
162 ieta += EMCAL_COLS;
163 }
164
165 double etai = (double)ieta;
166 double phii = (double)iphi;
167 double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
168
169 if (w > 0.0) {
170 phiMean += phii * w;
171 etaMean += etai * w;
172 wtot += w;
173 }
174 }
175 }
176
177 if (wtot > 0) {
178 phiMean /= wtot;
179 etaMean /= wtot;
180 } else {
181 LOG(error) << Form("Wrong weight %f\n", wtot);
182 }
183
184 // Calculate dispersion
185 for (auto iInput : inputsIndices) {
186
187 if (clusterAnalysis.E() > 0 && mInputsContainer[iInput].getEnergy() > 0) {
188 auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
189 auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta);
190
191 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
192 // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
193 if (mSharedCluster && nSupMod % 2) {
194 ieta += EMCAL_COLS;
195 }
196
197 double etai = (double)ieta;
198 double phii = (double)iphi;
199 double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
200
201 if (w > 0.0) {
202 nstat++;
203 d += w * ((etai - etaMean) * (etai - etaMean) + (phii - phiMean) * (phii - phiMean));
204 }
205 }
206 }
207
208 if (wtot > 0 && nstat > 1) {
209 d /= wtot;
210 } else {
211 d = 0.;
212 }
213
214 clusterAnalysis.setDispersion(TMath::Sqrt(d));
215}
216
219//____________________________________________________________________________
220template <class InputType>
221void ClusterFactory<InputType>::evalLocalPosition(gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
222{
223
224 int nstat = 0;
225
226 double dist = tMaxInCm(double(clusterAnalysis.E()));
227
228 double clXYZ[3] = {0., 0., 0.}, clRmsXYZ[3] = {0., 0., 0.}, xyzi[3], wtot = 0., w = 0.;
229
230 for (auto iInput : inputsIndices) {
231
232 try {
233 mGeomPtr->RelPosCellInSModule(mInputsContainer[iInput].getTower(), dist).GetCoordinates(xyzi[0], xyzi[1], xyzi[2]);
234 } catch (InvalidCellIDException& e) {
235 LOG(error) << e.what();
236 continue;
237 }
238
239 // Temporal patch, due to mapping problem, need to swap "y" in one of the 2 SM, although no effect in position calculation. GCB 05/2010
240 if (mSharedCluster && mSuperModuleNumber != mGeomPtr->GetSuperModuleNumber(mInputsContainer[iInput].getTower())) {
241 xyzi[1] *= -1;
242 }
243
244 if (mLogWeight > 0.0) {
245 w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
246 } else {
247 w = mInputsContainer[iInput].getEnergy(); // just energy
248 }
249
250 if (w > 0.0) {
251 wtot += w;
252 nstat++;
253
254 for (int i = 0; i < 3; i++) {
255 clXYZ[i] += (w * xyzi[i]);
256 clRmsXYZ[i] += (w * xyzi[i] * xyzi[i]);
257 }
258 } // w > 0
259 } // dig loop
260
261 // cout << " wtot " << wtot << endl;
262
263 if (wtot > 0) {
264 // xRMS = TMath::Sqrt(x2m - xMean*xMean);
265 for (int i = 0; i < 3; i++) {
266 clXYZ[i] /= wtot;
267
268 if (nstat > 1) {
269 clRmsXYZ[i] /= (wtot * wtot);
270 clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i] * clXYZ[i];
271
272 if (clRmsXYZ[i] > 0.0) {
273 clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
274 } else {
275 clRmsXYZ[i] = 0;
276 }
277 } else {
278 clRmsXYZ[i] = 0;
279 }
280 }
281 } else {
282 for (int i = 0; i < 3; i++) {
283 clXYZ[i] = clRmsXYZ[i] = -1.;
284 }
285 }
286
287 clusterAnalysis.setLocalPosition(math_utils::Point3D<float>(clXYZ[0], clXYZ[1], clXYZ[2]));
288}
289
292//____________________________________________________________________________
293template <class InputType>
294void ClusterFactory<InputType>::evalGlobalPosition(gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
295{
296
297 int i = 0, nstat = 0;
298
299 double dist = tMaxInCm(double(clusterAnalysis.E()));
300
301 double clXYZ[3] = {0., 0., 0.}, clRmsXYZ[3] = {0., 0., 0.}, lxyzi[3], xyzi[3], wtot = 0., w = 0.;
302
303 for (auto iInput : inputsIndices) {
304
305 // get the local coordinates of the cell
306 try {
307 mGeomPtr->RelPosCellInSModule(mInputsContainer[iInput].getTower(), dist).GetCoordinates(lxyzi[0], lxyzi[1], lxyzi[2]);
308 } catch (InvalidCellIDException& e) {
309 LOG(error) << e.what();
310 continue;
311 }
312
313 // Now get the global coordinate
314 mGeomPtr->GetGlobal(lxyzi, xyzi, mGeomPtr->GetSuperModuleNumber(mInputsContainer[iInput].getTower()));
315
316 if (mLogWeight > 0.0) {
317 w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
318 } else {
319 w = mInputsContainer[iInput].getEnergy(); // just energy
320 }
321
322 if (w > 0.0) {
323 wtot += w;
324 nstat++;
325
326 for (i = 0; i < 3; i++) {
327 clXYZ[i] += (w * xyzi[i]);
328 clRmsXYZ[i] += (w * xyzi[i] * xyzi[i]);
329 }
330 }
331 }
332
333 // cout << " wtot " << wtot << endl;
334
335 if (wtot > 0) {
336 // xRMS = TMath::Sqrt(x2m - xMean*xMean);
337 for (i = 0; i < 3; i++) {
338 clXYZ[i] /= wtot;
339
340 if (nstat > 1) {
341 clRmsXYZ[i] /= (wtot * wtot);
342 clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i] * clXYZ[i];
343
344 if (clRmsXYZ[i] > 0.0) {
345 clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
346 } else {
347 clRmsXYZ[i] = 0;
348 }
349 } else {
350 clRmsXYZ[i] = 0;
351 }
352 }
353 } else {
354 for (i = 0; i < 3; i++) {
355 clXYZ[i] = clRmsXYZ[i] = -1.;
356 }
357 }
358
359 clusterAnalysis.setGlobalPosition(math_utils::Point3D<float>(clXYZ[0], clXYZ[1], clXYZ[2]));
360}
361
364//____________________________________________________________________________
365template <class InputType>
366void ClusterFactory<InputType>::evalLocalPositionFit(double deff, double mLogWeight,
367 double phiSlope, gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
368{
369 int i = 0, nstat = 0;
370 double clXYZ[3] = {0., 0., 0.}, clRmsXYZ[3] = {0., 0., 0.}, xyzi[3], wtot = 0., w = 0.;
371
372 for (auto iInput : inputsIndices) {
373
374 try {
375 mGeomPtr->RelPosCellInSModule(mInputsContainer[iInput].getTower(), deff).GetCoordinates(xyzi[0], xyzi[1], xyzi[2]);
376 } catch (InvalidCellIDException& e) {
377 LOG(error) << e.what();
378 continue;
379 }
380
381 if (mLogWeight > 0.0) {
382 w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
383 } else {
384 w = mInputsContainer[iInput].getEnergy(); // just energy
385 }
386
387 if (w > 0.0) {
388 wtot += w;
389 nstat++;
390
391 for (i = 0; i < 3; i++) {
392 clXYZ[i] += (w * xyzi[i]);
393 clRmsXYZ[i] += (w * xyzi[i] * xyzi[i]);
394 }
395 }
396 } // loop
397
398 // cout << " wtot " << wtot << endl;
399
400 if (wtot > 0) {
401 // xRMS = TMath::Sqrt(x2m - xMean*xMean);
402 for (i = 0; i < 3; i++) {
403 clXYZ[i] /= wtot;
404
405 if (nstat > 1) {
406 clRmsXYZ[i] /= (wtot * wtot);
407 clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i] * clXYZ[i];
408
409 if (clRmsXYZ[i] > 0.0) {
410 clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
411 } else {
412 clRmsXYZ[i] = 0;
413 }
414 } else {
415 clRmsXYZ[i] = 0;
416 }
417 }
418 } else {
419 for (i = 0; i < 3; i++) {
420 clXYZ[i] = clRmsXYZ[i] = -1.;
421 }
422 }
423
424 // clRmsXYZ[i] ??
425
426 if (phiSlope != 0.0 && mLogWeight > 0.0 && wtot) {
427 // Correction in phi direction (y - coords here); Aug 16;
428 // May be put to global level or seperate method
429 double ycorr = clXYZ[1] * (1. + phiSlope);
430
431 // printf(" y %f : ycorr %f : slope %f \n", clXYZ[1], ycorr, phiSlope);
432 clXYZ[1] = ycorr;
433 }
434
435 clusterAnalysis.setLocalPosition(math_utils::Point3D<float>(clXYZ[0], clXYZ[1], clXYZ[2]));
436}
437
442//_____________________________________________________________________________
443template <class InputType>
444void ClusterFactory<InputType>::getDeffW0(const double esum, double& deff, double& w0)
445{
446 double e = 0.0;
447 const double kdp0 = 9.25147, kdp1 = 1.16700; // Hard coded now
448 const double kwp0 = 4.83713, kwp1 = -2.77970e-01, kwp2 = 4.41116;
449
450 // No extrapolation here
451 e = esum < 0.5 ? 0.5 : esum;
452 e = e > 100. ? 100. : e;
453
454 deff = kdp0 + kdp1 * TMath::Log(e);
455 w0 = kwp0 / (1. + TMath::Exp(kwp1 * (e + kwp2)));
456}
457
465//______________________________________________________________________________
466template <class InputType>
467void ClusterFactory<InputType>::evalCoreEnergy(gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
468{
469
470 float coreEnergy = 0.;
471
472 if (!clusterAnalysis.getLocalPosition().Mag2()) {
473 evalLocalPosition(inputsIndices, clusterAnalysis);
474 }
475
476 double phiPoint = clusterAnalysis.getLocalPosition().Phi();
477 double etaPoint = clusterAnalysis.getLocalPosition().Eta();
478 for (auto iInput : inputsIndices) {
479
480 auto [eta, phi] = mGeomPtr->EtaPhiFromIndex(mInputsContainer[iInput].getTower());
481 phi = phi * TMath::DegToRad();
482
483 double distance = TMath::Sqrt((eta - etaPoint) * (eta - etaPoint) + (phi - phiPoint) * (phi - phiPoint));
484
485 if (distance < mCoreRadius) {
486 coreEnergy += mInputsContainer[iInput].getEnergy();
487 }
488 }
489 clusterAnalysis.setCoreEnergy(coreEnergy);
490}
491
495//____________________________________________________________________________
496template <class InputType>
497void ClusterFactory<InputType>::evalElipsAxis(gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
498{
499 double wtot = 0.;
500 double x = 0.;
501 double z = 0.;
502 double dxx = 0.;
503 double dzz = 0.;
504 double dxz = 0.;
505
506 std::array<float, 2> lambda;
507
508 for (auto iInput : inputsIndices) {
509
510 auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
511 auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta);
512
513 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
514 // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
515 if (mSharedCluster && nSupMod % 2) {
516 ieta += EMCAL_COLS;
517 }
518
519 double etai = (double)ieta;
520 double phii = (double)iphi;
521
522 double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
523 // clusterAnalysis.E() summed amplitude of inputs, i.e. energy of cluster
524 // Gives smaller value of lambda than log weight
525 // w = mEnergyList[iInput] / clusterAnalysis.E(); // Nov 16, 2006 - try just energy
526
527 dxx += w * etai * etai;
528 x += w * etai;
529 dzz += w * phii * phii;
530 z += w * phii;
531
532 dxz += w * etai * phii;
533
534 wtot += w;
535 }
536
537 if (wtot > 0) {
538 dxx /= wtot;
539 x /= wtot;
540 dxx -= x * x;
541 dzz /= wtot;
542 z /= wtot;
543 dzz -= z * z;
544 dxz /= wtot;
545 dxz -= x * z;
546
547 lambda[0] = 0.5 * (dxx + dzz) + TMath::Sqrt(0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz);
548
549 if (lambda[0] > 0) {
550 lambda[0] = TMath::Sqrt(lambda[0]);
551 } else {
552 lambda[0] = 0;
553 }
554
555 lambda[1] = 0.5 * (dxx + dzz) - TMath::Sqrt(0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz);
556
557 if (lambda[1] > 0) { // To avoid exception if numerical errors lead to negative lambda.
558 lambda[1] = TMath::Sqrt(lambda[1]);
559 } else {
560 lambda[1] = 0.;
561 }
562 } else {
563 lambda[0] = 0.;
564 lambda[1] = 0.;
565 }
566
567 clusterAnalysis.setM02(lambda[0] * lambda[0]);
568 clusterAnalysis.setM20(lambda[1] * lambda[1]);
569}
570
573//____________________________________________________________________________
574template <class InputType>
575std::tuple<int, float, float, bool> ClusterFactory<InputType>::getMaximalEnergyIndex(gsl::span<const int> inputsIndices) const
576{
577
578 float energy = 0.;
579 int mid = 0;
580 float cellAmp = 0;
581 int iSupMod0 = -1;
582 bool shared = false;
583 for (auto iInput : inputsIndices) {
584 if (iInput >= mInputsContainer.size()) {
585 throw CellIndexRangeException(iInput, mInputsContainer.size());
586 }
587 cellAmp += mInputsContainer[iInput].getEnergy();
588 if (iSupMod0 == -1) {
589 iSupMod0 = mGeomPtr->GetSuperModuleNumber(mInputsContainer[iInput].getTower());
590 } else if (iSupMod0 != mGeomPtr->GetSuperModuleNumber(mInputsContainer[iInput].getTower())) {
591 shared = true;
592 }
593 if (mInputsContainer[iInput].getEnergy() > energy) {
594 energy = mInputsContainer[iInput].getEnergy();
595 mid = iInput;
596 }
597 } // loop on cluster inputs
598
599 return std::make_tuple(mid, energy, cellAmp, shared);
600}
601
604//____________________________________________________________________________
605template <class InputType>
606bool ClusterFactory<InputType>::isExoticCell(short towerId, float ecell, float const exoticTime, float& fCross) const
607{
608 if (ecell < mExoticCellMinAmplitude) {
609 return false; // do not reject low energy cells
610 }
611
612 // if the look up table is not set yet (mostly due to a reset call) then set it up now.
613 if (!getLookUpInit()) {
615 }
616
617 float eCross = getECross(towerId, ecell, exoticTime);
618 fCross = 1.f - eCross / ecell;
619
620 if (fCross > mExoticCellFraction) {
621 LOG(debug) << "EXOTIC CELL id " << towerId << ", eCell " << ecell << ", eCross " << eCross << ", 1-eCross/eCell " << 1 - eCross / ecell;
622 return true;
623 }
624
625 return false;
626}
627
630//____________________________________________________________________________
631template <class InputType>
632float ClusterFactory<InputType>::getECross(short towerId, float energy, float const exoticTime) const
633{
634 auto [iSM, iMod, iIphi, iIeta] = mGeomPtr->GetCellIndex(towerId);
635 auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(iSM, iMod, iIphi, iIeta);
636
637 // Get close cells index, energy and time, not in corners
638
639 short towerId1 = -1;
640 short towerId2 = -1;
641
642 if (iphi < o2::emcal::EMCAL_ROWS - 1) {
643 try {
644 towerId1 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi + 1, ieta);
645 } catch (InvalidCellIDException& e) {
646 towerId1 = -1 * e.getCellID();
647 }
648 }
649 if (iphi > 0) {
650 try {
651 towerId2 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi - 1, ieta);
652 } catch (InvalidCellIDException& e) {
653 towerId2 = -1 * e.getCellID();
654 }
655 }
656
657 // In case of cell in eta = 0 border, depending on SM shift the cross cell index
658
659 short towerId3 = -1;
660 short towerId4 = -1;
661
662 if (ieta == o2::emcal::EMCAL_COLS - 1 && !(iSM % 2)) {
663 try {
664 towerId3 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM + 1, iphi, 0);
665 } catch (InvalidCellIDException& e) {
666 towerId3 = -1 * e.getCellID();
667 }
668 try {
669 towerId4 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta - 1);
670 } catch (InvalidCellIDException& e) {
671 towerId4 = -1 * e.getCellID();
672 }
673 } else if (ieta == 0 && iSM % 2) {
674 try {
675 towerId3 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta + 1);
676 } catch (InvalidCellIDException& e) {
677 towerId3 = -1 * e.getCellID();
678 }
679 try {
680 towerId4 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM - 1, iphi, o2::emcal::EMCAL_COLS - 1);
681 } catch (InvalidCellIDException& e) {
682 towerId4 = -1 * e.getCellID();
683 }
684 } else {
685 if (ieta < o2::emcal::EMCAL_COLS - 1) {
686 try {
687 towerId3 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta + 1);
688 } catch (InvalidCellIDException& e) {
689 towerId3 = -1 * e.getCellID();
690 }
691 }
692 if (ieta > 0) {
693 try {
694 towerId4 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta - 1);
695 } catch (InvalidCellIDException& e) {
696 towerId4 = -1 * e.getCellID();
697 }
698 }
699 }
700
701 LOG(debug) << "iSM " << iSM << ", towerId " << towerId << ", a " << towerId1 << ", b " << towerId2 << ", c " << towerId3 << ", e " << towerId3;
702
703 short index1 = (towerId1 > -1) ? mLoolUpTowerToIndex.at(towerId1) : -1;
704 short index2 = (towerId2 > -1) ? mLoolUpTowerToIndex.at(towerId2) : -1;
705 short index3 = (towerId3 > -1) ? mLoolUpTowerToIndex.at(towerId3) : -1;
706 short index4 = (towerId4 > -1) ? mLoolUpTowerToIndex.at(towerId4) : -1;
707
708 std::array<std::pair<float, float>, 4> cellData = {
709 {{(index1 > -1) ? mInputsContainer[index1].getEnergy() : 0., (index1 > -1) ? mInputsContainer[index1].getTimeStamp() : 0.},
710 {(index2 > -1) ? mInputsContainer[index2].getEnergy() : 0., (index2 > -1) ? mInputsContainer[index2].getTimeStamp() : 0.},
711 {(index3 > -1) ? mInputsContainer[index3].getEnergy() : 0., (index3 > -1) ? mInputsContainer[index3].getTimeStamp() : 0.},
712 {(index4 > -1) ? mInputsContainer[index4].getEnergy() : 0., (index4 > -1) ? mInputsContainer[index4].getTimeStamp() : 0.}}};
713
714 for (auto& cell : cellData) {
715 if (std::abs(exoticTime - cell.second) > mExoticCellDiffTime) {
716 cell.first = 0;
717 }
718 }
719
720 float w1 = 1, w2 = 1, w3 = 1, w4 = 1;
721 if (mUseWeightExotic) {
722 w1 = GetCellWeight(cellData[0].first, energy);
723 w2 = GetCellWeight(cellData[1].first, energy);
724 w3 = GetCellWeight(cellData[2].first, energy);
725 w4 = GetCellWeight(cellData[3].first, energy);
726 }
727
728 if (cellData[0].first < mExoticCellInCrossMinAmplitude || w1 <= 0) {
729 cellData[0].first = 0;
730 }
731 if (cellData[1].first < mExoticCellInCrossMinAmplitude || w2 <= 0) {
732 cellData[1].first = 0;
733 }
734 if (cellData[2].first < mExoticCellInCrossMinAmplitude || w3 <= 0) {
735 cellData[2].first = 0;
736 }
737 if (cellData[3].first < mExoticCellInCrossMinAmplitude || w4 <= 0) {
738 cellData[3].first = 0;
739 }
740
741 return cellData[0].first + cellData[1].first + cellData[2].first + cellData[3].first;
742}
743
746//____________________________________________________________________________
747template <class InputType>
748float ClusterFactory<InputType>::GetCellWeight(float eCell, float eCluster) const
749{
750 if (eCell > 0 && eCluster > 0) {
751 if (mLogWeight > 0) {
752 return std::max(0.f, mLogWeight + std::log(eCell / eCluster));
753 } else {
754 return std::log(eCluster / eCell);
755 }
756 } else {
757 return 0.;
758 }
759}
760
763//____________________________________________________________________________
764template <class InputType>
765int ClusterFactory<InputType>::getMultiplicityAtLevel(float H, gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
766{
767 int multipl = 0;
768 for (auto iInput : inputsIndices) {
769 if (mInputsContainer[iInput].getEnergy() > H * clusterAnalysis.E()) {
770 multipl++;
771 }
772 }
773
774 return multipl;
775}
776
779//____________________________________________________________________________
780template <class InputType>
781void ClusterFactory<InputType>::evalTime(gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
782{
783 float maxE = 0;
784 unsigned short maxAt = 0;
785 for (auto iInput : inputsIndices) {
786 if (mInputsContainer[iInput].getEnergy() > maxE) {
787 maxE = mInputsContainer[iInput].getEnergy();
788 maxAt = iInput;
789 }
790 }
791
792 clusterAnalysis.setClusterTime(mInputsContainer[maxAt].getTimeStamp());
793}
794
798//_____________________________________________________________________
799template <class InputType>
800double ClusterFactory<InputType>::tMaxInCm(const double e, const int key) const
801{
802 const double ca = 4.82; // shower max parameter - first guess; ca=TMath::Log(1000./8.07)
803 double tmax = 0.; // position of electromagnetic shower max in cm
804
805 const double x0 = 1.31; // radiation lenght (cm)
806
807 if (e > 0.1) {
808 tmax = TMath::Log(e) + ca;
809 if (key == 0) {
810 tmax += 0.5;
811 } else {
812 tmax -= 0.5;
813 }
814 tmax *= x0; // convert to cm
815 }
816
817 return tmax;
818}
819
822//______________________________________________________________________________
823template <class InputType>
825{
826 return (2. * TMath::ATan(TMath::Exp(-arg)));
827}
828
831//______________________________________________________________________________
832template <class InputType>
834{
835 return (-1 * TMath::Log(TMath::Tan(0.5 * arg)));
836}
837
838template <class InputType>
839ClusterFactory<InputType>::ClusterIterator::ClusterIterator(const ClusterFactory& factory, int clusterIndex, bool forward) : mClusterFactory(factory),
840 mCurrentCluster(),
841 mClusterID(clusterIndex),
842 mForward(forward)
843{
844 mCurrentCluster = mClusterFactory.buildCluster(mClusterID);
845}
846
847template <class InputType>
849{
850 return &mClusterFactory == &rhs.mClusterFactory && mClusterID == rhs.mClusterID && mForward == rhs.mForward;
851}
852
853template <class InputType>
855{
856 if (mForward) {
857 mClusterID++;
858 } else {
859 mClusterID--;
860 }
861 mCurrentCluster = mClusterFactory.buildCluster(mClusterID);
862 return *this;
863}
864
865template <class InputType>
867{
868 auto tmp = *this;
869 ++(*this);
870 return tmp;
871}
872
873template <class InputType>
875{
876 if (mForward) {
877 mClusterID--;
878 } else {
879 mClusterID++;
880 }
881 mCurrentCluster = mClusterFactory.buildCluster(mClusterID);
882 return *this;
883}
884
885template <class InputType>
887{
888 auto tmp = *this;
889 --(*this);
890 return tmp;
891}
892
int32_t i
std::ostringstream debug
StringRef key
Cluster class for kinematic cluster parametersported from AliVCluster in AliRoot.
void setGlobalPosition(math_utils::Point3D< float > x)
Set the cluster global position.
void setFCross(float fCross)
void setCoreEnergy(float energy)
math_utils::Point3D< float > getLocalPosition() const
void setCellsIndices(const std::vector< unsigned short > &array)
Set the array of cell indices.
void setIndMaxInput(const int ind)
void setClusterTime(float time)
void setLocalPosition(math_utils::Point3D< float > x)
ClusterIterator & operator--()
Prefix decrementation operator.
ClusterIterator(const ClusterFactory &factory, int clusterIndex, bool forward)
Constructor, initializing the iterator.
bool operator==(const ClusterIterator &rhs) const
Check for equalness.
ClusterIterator & operator++()
Prefix incrementation operator.
Exception handling uninitialized look up table.
const char * what() const noexcept final
Access to error message of the exception.
EMCal clusters factory Ported from class AliEMCALcluster.
void reset()
Reset containers.
std::tuple< int, float, float, bool > getMaximalEnergyIndex(gsl::span< const int > inputsIndices) const
Finds the maximum energy in the cluster and computes the Summed amplitude of digits/cells.
float GetCellWeight(float eCell, float eCluster) const
return weight of cell for shower shape calculation
float thetaToEta(float arg) const
Converts Theta (Radians) to Eta (Radians)
float getECross(short absID, float energy, float const exoticTime) const
Calculate the energy in the cross around the energy of a given cell.
void evalGlobalPosition(gsl::span< const int > inputsIndices, AnalysisCluster &cluster) const
Calculates the center of gravity in the global ALICE coordinates.
void evalLocalPosition(gsl::span< const int > inputsIndices, AnalysisCluster &cluster) const
Calculates the center of gravity in the local EMCAL-module coordinates.
Double_t tMaxInCm(const Double_t e=0.0, const int key=0) const
static void getDeffW0(const Double_t esum, Double_t &deff, Double_t &w0)
void evalTime(gsl::span< const int > inputsIndices, AnalysisCluster &clusterAnalysis) const
Time is set to the time of the digit with the maximum energy.
int getMultiplicityAtLevel(float level, gsl::span< const int > inputsIndices, AnalysisCluster &clusterAnalysis) const
Calculates the multiplicity of digits/cells with energy larger than level*energy.
ClusterFactory()=default
Dummy constructor.
void evalLocalPositionFit(Double_t deff, Double_t w0, Double_t phiSlope, gsl::span< const int > inputsIndices, AnalysisCluster &cluster) const
evaluates local position of clusters in SM
void evalElipsAxis(gsl::span< const int > inputsIndices, AnalysisCluster &clusterAnalysis) const
void evalCoreEnergy(gsl::span< const int > inputsIndices, AnalysisCluster &clusterAnalysis) const
void evalDispersion(gsl::span< const int > inputsIndices, AnalysisCluster &clusterAnalysis) const
float etaToTheta(float arg) const
Converts Eta (Radians) to Theta (Radians)
bool isExoticCell(short towerId, float ecell, float const exoticTime, float &fCross) const
Look to cell neighbourhood and reject if it seems exotic.
AnalysisCluster buildCluster(int index, o2::emcal::ClusterLabel *clusterLabel=nullptr) const
evaluates cluster parameters: position, shower shape, primaries ...
cluster class for MC particle IDs and their respective energy fraction
void orderLabels()
Sort the labels and energy fraction in descending order (largest energy fraction to smallest)
void normalize(float factor)
Normalize the energy fraction.
void addValue(int label, float energyFraction)
Add label and energy fraction to the.
Exception handling non-existing cell IDs.
int getCellID() const noexcept
Access to cell ID raising the exception.
const char * what() const noexcept final
Access to error message of the exception.
GLint GLenum GLint x
Definition glcorearb.h:403
GLsizei GLsizei GLfloat distance
Definition glcorearb.h:5506
GLuint GLfloat x0
Definition glcorearb.h:5034
GLubyte GLubyte GLubyte GLubyte w
Definition glcorearb.h:852
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
@ EMCAL_ROWS
Number of rows per module for EMCAL.
Definition Constants.h:25
@ EMCAL_COLS
Number of columns per module for EMCAL.
Definition Constants.h:26
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"