30template <
class InputType>
33 setContainer(clustersContainer, inputsContainer, cellsIndices);
36template <
class InputType>
39 mClustersContainer = gsl::span<const o2::emcal::Cluster>();
40 mInputsContainer = gsl::span<const InputType>();
41 mCellsIndices = gsl::span<int>();
43 mCellLabelContainer = gsl::span<const o2::emcal::CellLabel>();
49template <
class InputType>
52 if (clusterIndex >= mClustersContainer.size()) {
60 clusterAnalysis.
setID(clusterIndex);
62 int firstCellIndex = mClustersContainer[clusterIndex].getCellIndexFirst();
63 int nCells = mClustersContainer[clusterIndex].getNCells();
65 gsl::span<const int> inputsIndices = gsl::span<const int>(&mCellsIndices[firstCellIndex], nCells);
70 auto [inputIndMax, inputEnergyMax, cellAmp, shared] = getMaximalEnergyIndex(inputsIndices);
72 short towerId = mInputsContainer[inputIndMax].getTower();
74 float exoticTime = mInputsContainer[inputIndMax].getTimeStamp();
79 clusterAnalysis.
setIsExotic(isExoticCell(towerId, inputEnergyMax, exoticTime, fCross));
87 clusterAnalysis.
setE(cellAmp);
89 mSuperModuleNumber = mGeomPtr->GetSuperModuleNumber(towerId);
91 clusterAnalysis.
setNCells(inputsIndices.size());
93 std::vector<unsigned short> cellsIdices;
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) {
103 clusterLabel->
addValue(mCellLabelContainer[cellIndex].GetLabel(iLabel),
104 mCellLabelContainer[cellIndex].GetAmplitudeFraction(iLabel) * mInputsContainer[cellIndex].getEnergy());
108 if (addClusterLabels) {
116 evalGlobalPosition(inputsIndices, clusterAnalysis);
117 evalLocalPosition(inputsIndices, clusterAnalysis);
120 evalElipsAxis(inputsIndices, clusterAnalysis);
121 evalDispersion(inputsIndices, clusterAnalysis);
123 evalCoreEnergy(inputsIndices, clusterAnalysis);
124 evalTime(inputsIndices, clusterAnalysis);
136 return clusterAnalysis;
143template <
class InputType>
146 double d = 0., wtot = 0.;
150 double etaMean = 0.0, phiMean = 0.0;
153 for (
auto iInput : inputsIndices) {
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);
161 if (mSharedCluster && nSupMod % 2) {
165 double etai = (double)ieta;
166 double phii = (double)iphi;
167 double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.
E()));
181 LOG(error) << Form(
"Wrong weight %f\n", wtot);
185 for (
auto iInput : inputsIndices) {
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);
193 if (mSharedCluster && nSupMod % 2) {
197 double etai = (double)ieta;
198 double phii = (double)iphi;
199 double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.
E()));
203 d +=
w * ((etai - etaMean) * (etai - etaMean) + (phii - phiMean) * (phii - phiMean));
208 if (wtot > 0 && nstat > 1) {
220template <
class InputType>
226 double dist = tMaxInCm(
double(clusterAnalysis.
E()));
228 double clXYZ[3] = {0., 0., 0.}, clRmsXYZ[3] = {0., 0., 0.}, xyzi[3], wtot = 0.,
w = 0.;
230 for (
auto iInput : inputsIndices) {
233 mGeomPtr->RelPosCellInSModule(mInputsContainer[iInput].getTower(), dist).GetCoordinates(xyzi[0], xyzi[1], xyzi[2]);
240 if (mSharedCluster && mSuperModuleNumber != mGeomPtr->GetSuperModuleNumber(mInputsContainer[iInput].getTower())) {
244 if (mLogWeight > 0.0) {
245 w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.
E()));
247 w = mInputsContainer[iInput].getEnergy();
254 for (
int i = 0;
i < 3;
i++) {
255 clXYZ[
i] += (
w * xyzi[
i]);
256 clRmsXYZ[
i] += (
w * xyzi[
i] * xyzi[
i]);
265 for (
int i = 0;
i < 3;
i++) {
269 clRmsXYZ[
i] /= (wtot * wtot);
270 clRmsXYZ[
i] = clRmsXYZ[
i] - clXYZ[
i] * clXYZ[
i];
272 if (clRmsXYZ[
i] > 0.0) {
273 clRmsXYZ[
i] = TMath::Sqrt(clRmsXYZ[
i]);
282 for (
int i = 0;
i < 3;
i++) {
283 clXYZ[
i] = clRmsXYZ[
i] = -1.;
293template <
class InputType>
297 int i = 0, nstat = 0;
299 double dist = tMaxInCm(
double(clusterAnalysis.
E()));
301 double clXYZ[3] = {0., 0., 0.}, clRmsXYZ[3] = {0., 0., 0.}, lxyzi[3], xyzi[3], wtot = 0.,
w = 0.;
303 for (
auto iInput : inputsIndices) {
307 mGeomPtr->RelPosCellInSModule(mInputsContainer[iInput].getTower(), dist).GetCoordinates(lxyzi[0], lxyzi[1], lxyzi[2]);
314 mGeomPtr->GetGlobal(lxyzi, xyzi, mGeomPtr->GetSuperModuleNumber(mInputsContainer[iInput].getTower()));
316 if (mLogWeight > 0.0) {
317 w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.
E()));
319 w = mInputsContainer[iInput].getEnergy();
326 for (
i = 0;
i < 3;
i++) {
327 clXYZ[
i] += (
w * xyzi[
i]);
328 clRmsXYZ[
i] += (
w * xyzi[
i] * xyzi[
i]);
337 for (
i = 0;
i < 3;
i++) {
341 clRmsXYZ[
i] /= (wtot * wtot);
342 clRmsXYZ[
i] = clRmsXYZ[
i] - clXYZ[
i] * clXYZ[
i];
344 if (clRmsXYZ[
i] > 0.0) {
345 clRmsXYZ[
i] = TMath::Sqrt(clRmsXYZ[
i]);
354 for (
i = 0;
i < 3;
i++) {
355 clXYZ[
i] = clRmsXYZ[
i] = -1.;
365template <
class InputType>
367 double phiSlope, gsl::span<const int> inputsIndices,
AnalysisCluster& clusterAnalysis)
const
369 int i = 0, nstat = 0;
370 double clXYZ[3] = {0., 0., 0.}, clRmsXYZ[3] = {0., 0., 0.}, xyzi[3], wtot = 0.,
w = 0.;
372 for (
auto iInput : inputsIndices) {
375 mGeomPtr->RelPosCellInSModule(mInputsContainer[iInput].getTower(), deff).GetCoordinates(xyzi[0], xyzi[1], xyzi[2]);
381 if (mLogWeight > 0.0) {
382 w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.
E()));
384 w = mInputsContainer[iInput].getEnergy();
391 for (
i = 0;
i < 3;
i++) {
392 clXYZ[
i] += (
w * xyzi[
i]);
393 clRmsXYZ[
i] += (
w * xyzi[
i] * xyzi[
i]);
402 for (
i = 0;
i < 3;
i++) {
406 clRmsXYZ[
i] /= (wtot * wtot);
407 clRmsXYZ[
i] = clRmsXYZ[
i] - clXYZ[
i] * clXYZ[
i];
409 if (clRmsXYZ[
i] > 0.0) {
410 clRmsXYZ[
i] = TMath::Sqrt(clRmsXYZ[
i]);
419 for (
i = 0;
i < 3;
i++) {
420 clXYZ[
i] = clRmsXYZ[
i] = -1.;
426 if (phiSlope != 0.0 && mLogWeight > 0.0 && wtot) {
429 double ycorr = clXYZ[1] * (1. + phiSlope);
443template <
class InputType>
447 const double kdp0 = 9.25147, kdp1 = 1.16700;
448 const double kwp0 = 4.83713, kwp1 = -2.77970e-01, kwp2 = 4.41116;
451 e = esum < 0.5 ? 0.5 : esum;
452 e = e > 100. ? 100. : e;
454 deff = kdp0 + kdp1 * TMath::Log(e);
455 w0 = kwp0 / (1. + TMath::Exp(kwp1 * (e + kwp2)));
466template <
class InputType>
470 float coreEnergy = 0.;
473 evalLocalPosition(inputsIndices, clusterAnalysis);
478 for (
auto iInput : inputsIndices) {
480 auto [eta, phi] = mGeomPtr->EtaPhiFromIndex(mInputsContainer[iInput].getTower());
481 phi = phi * TMath::DegToRad();
483 double distance = TMath::Sqrt((eta - etaPoint) * (eta - etaPoint) + (phi - phiPoint) * (phi - phiPoint));
486 coreEnergy += mInputsContainer[iInput].getEnergy();
496template <
class InputType>
506 std::array<float, 2> lambda;
508 for (
auto iInput : inputsIndices) {
510 auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
511 auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta);
515 if (mSharedCluster && nSupMod % 2) {
519 double etai = (double)ieta;
520 double phii = (double)iphi;
522 double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.
E()));
527 dxx +=
w * etai * etai;
529 dzz +=
w * phii * phii;
532 dxz +=
w * etai * phii;
547 lambda[0] = 0.5 * (dxx + dzz) + TMath::Sqrt(0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz);
550 lambda[0] = TMath::Sqrt(lambda[0]);
555 lambda[1] = 0.5 * (dxx + dzz) - TMath::Sqrt(0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz);
558 lambda[1] = TMath::Sqrt(lambda[1]);
567 clusterAnalysis.
setM02(lambda[0] * lambda[0]);
568 clusterAnalysis.
setM20(lambda[1] * lambda[1]);
574template <
class InputType>
583 for (
auto iInput : inputsIndices) {
584 if (iInput >= mInputsContainer.size()) {
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())) {
593 if (mInputsContainer[iInput].getEnergy() > energy) {
594 energy = mInputsContainer[iInput].getEnergy();
599 return std::make_tuple(mid, energy, cellAmp, shared);
605template <
class InputType>
608 if (ecell < mExoticCellMinAmplitude) {
613 if (!getLookUpInit()) {
617 float eCross = getECross(towerId, ecell, exoticTime);
618 fCross = 1.f - eCross / ecell;
620 if (fCross > mExoticCellFraction) {
621 LOG(
debug) <<
"EXOTIC CELL id " << towerId <<
", eCell " << ecell <<
", eCross " << eCross <<
", 1-eCross/eCell " << 1 - eCross / ecell;
631template <
class InputType>
634 auto [iSM, iMod, iIphi, iIeta] = mGeomPtr->GetCellIndex(towerId);
635 auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(iSM, iMod, iIphi, iIeta);
644 towerId1 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi + 1, ieta);
651 towerId2 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi - 1, ieta);
664 towerId3 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM + 1, iphi, 0);
669 towerId4 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta - 1);
673 }
else if (ieta == 0 && iSM % 2) {
675 towerId3 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta + 1);
687 towerId3 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta + 1);
694 towerId4 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta - 1);
701 LOG(
debug) <<
"iSM " << iSM <<
", towerId " << towerId <<
", a " << towerId1 <<
", b " << towerId2 <<
", c " << towerId3 <<
", e " << towerId3;
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;
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.}}};
714 for (
auto& cell : cellData) {
715 if (std::abs(exoticTime - cell.second) > mExoticCellDiffTime) {
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);
728 if (cellData[0].
first < mExoticCellInCrossMinAmplitude || w1 <= 0) {
729 cellData[0].first = 0;
731 if (cellData[1].
first < mExoticCellInCrossMinAmplitude || w2 <= 0) {
732 cellData[1].first = 0;
734 if (cellData[2].
first < mExoticCellInCrossMinAmplitude || w3 <= 0) {
735 cellData[2].first = 0;
737 if (cellData[3].
first < mExoticCellInCrossMinAmplitude || w4 <= 0) {
738 cellData[3].first = 0;
741 return cellData[0].first + cellData[1].first + cellData[2].first + cellData[3].first;
747template <
class InputType>
750 if (eCell > 0 && eCluster > 0) {
751 if (mLogWeight > 0) {
752 return std::max(0.f, mLogWeight + std::log(eCell / eCluster));
754 return std::log(eCluster / eCell);
764template <
class InputType>
768 for (
auto iInput : inputsIndices) {
769 if (mInputsContainer[iInput].getEnergy() > H * clusterAnalysis.
E()) {
780template <
class InputType>
784 unsigned short maxAt = 0;
785 for (
auto iInput : inputsIndices) {
786 if (mInputsContainer[iInput].getEnergy() > maxE) {
787 maxE = mInputsContainer[iInput].getEnergy();
792 clusterAnalysis.
setClusterTime(mInputsContainer[maxAt].getTimeStamp());
799template <
class InputType>
802 const double ca = 4.82;
805 const double x0 = 1.31;
808 tmax = TMath::Log(e) + ca;
823template <
class InputType>
826 return (2. * TMath::ATan(TMath::Exp(-arg)));
832template <
class InputType>
835 return (-1 * TMath::Log(TMath::Tan(0.5 * arg)));
838template <
class InputType>
841 mClusterID(clusterIndex),
844 mCurrentCluster = mClusterFactory.
buildCluster(mClusterID);
847template <
class InputType>
850 return &mClusterFactory == &rhs.mClusterFactory && mClusterID == rhs.mClusterID && mForward == rhs.mForward;
853template <
class InputType>
861 mCurrentCluster = mClusterFactory.buildCluster(mClusterID);
865template <
class InputType>
873template <
class InputType>
881 mCurrentCluster = mClusterFactory.buildCluster(mClusterID);
885template <
class InputType>
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 setDispersion(float disp)
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.
GLsizei GLsizei GLfloat distance
GLubyte GLubyte GLubyte GLubyte w
GLdouble GLdouble GLdouble z
@ EMCAL_ROWS
Number of rows per module for EMCAL.
@ EMCAL_COLS
Number of columns per module for EMCAL.
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"