ClusterFactory.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
#include "EMCALBase/ClusterFactory.h"
#include "DataFormatsEMCAL/Cluster.h"
#include "DataFormatsEMCAL/Digit.h"
#include "DataFormatsEMCAL/Cell.h"
#include "DataFormatsEMCAL/AnalysisCluster.h"
#include "DataFormatsEMCAL/Constants.h"
#include "DataFormatsEMCAL/CellLabel.h"
#include "DataFormatsEMCAL/ClusterLabel.h"
#include "EMCALBase/Geometry.h"
// #include "MathUtils/Cartesian.h"
 
#include <Rtypes.h>
 
#include <array>
#include <span>
 
using namespace o2::emcal;
 
template <class InputType>
ClusterFactory<InputType>::ClusterFactory(std::span<const o2::emcal::Cluster> clustersContainer, std::span<const InputType> inputsContainer, std::span<const int> cellsIndices)
{
  setContainer(clustersContainer, inputsContainer, cellsIndices);
}
 
template <class InputType>
void ClusterFactory<InputType>::reset()
{
  mClustersContainer = std::span<const o2::emcal::Cluster>();
  mInputsContainer = std::span<const InputType>();
  mCellsIndices = std::span<int>();
  mLookUpInit = false;
  mCellLabelContainer = std::span<const o2::emcal::CellLabel>();
}
 
//____________________________________________________________________________
template <class InputType>
o2::emcal::AnalysisCluster ClusterFactory<InputType>::buildCluster(int clusterIndex, o2::emcal::ClusterLabel* clusterLabel) const
{
  if (clusterIndex >= mClustersContainer.size()) {
    throw ClusterRangeException(clusterIndex, mClustersContainer.size());
  }
  if (!mGeomPtr) {
    throw GeometryNotSetException();
  }
 
  o2::emcal::AnalysisCluster clusterAnalysis;
  clusterAnalysis.setID(clusterIndex);
 
  int firstCellIndex = mClustersContainer[clusterIndex].getCellIndexFirst();
  int nCells = mClustersContainer[clusterIndex].getNCells();
 
  std::span<const int> inputsIndices = std::span<const int>(&mCellsIndices[firstCellIndex], nCells);
 
  // First calculate the index of input with maximum amplitude and get
  // the supermodule number where it sits.
 
  auto [inputIndMax, inputEnergyMax, cellAmp, shared] = getMaximalEnergyIndex(inputsIndices);
 
  short towerId = mInputsContainer[inputIndMax].getTower();
 
  float exoticTime = mInputsContainer[inputIndMax].getTimeStamp();
 
  float fCross = 0.;
 
  try {
    clusterAnalysis.setIsExotic(isExoticCell(towerId, inputEnergyMax, exoticTime, fCross));
    clusterAnalysis.setFCross(fCross);
  } catch (UninitLookUpTableException& e) {
    LOG(error) << e.what();
  }
 
  clusterAnalysis.setIndMaxInput(inputIndMax);
 
  clusterAnalysis.setE(cellAmp);
 
  mSuperModuleNumber = mGeomPtr->GetSuperModuleNumber(towerId);
 
  clusterAnalysis.setNCells(inputsIndices.size());
 
  std::vector<unsigned short> cellsIdices;
 
  bool addClusterLabels = ((clusterLabel != nullptr) && (mCellLabelContainer.size() > 0));
  for (auto cellIndex : inputsIndices) {
    cellsIdices.push_back(cellIndex);
    if (addClusterLabels) {
      for (size_t iLabel = 0; iLabel < mCellLabelContainer[cellIndex].GetLabelSize(); iLabel++) {
        if (mCellLabelContainer[cellIndex].GetAmplitudeFraction(iLabel) <= 0.f) {
          continue; // skip 0 entries
        }
        clusterLabel->addValue(mCellLabelContainer[cellIndex].GetLabel(iLabel),
                               mCellLabelContainer[cellIndex].GetAmplitudeFraction(iLabel) * mInputsContainer[cellIndex].getEnergy());
      }
    }
  }
  if (addClusterLabels) {
    clusterLabel->orderLabels();
    clusterLabel->normalize(cellAmp);
  }
 
  clusterAnalysis.setCellsIndices(cellsIdices);
 
  // evaluate global and local position
  evalGlobalPosition(inputsIndices, clusterAnalysis);
  evalLocalPosition(inputsIndices, clusterAnalysis);
 
  // evaluate shower parameters
  evalElipsAxis(inputsIndices, clusterAnalysis);
  evalDispersion(inputsIndices, clusterAnalysis);
 
  // evaluate number of local maxima
  evalNExMax(inputsIndices, clusterAnalysis);
 
  evalCoreEnergy(inputsIndices, clusterAnalysis);
  evalTime(inputsIndices, clusterAnalysis);
 
  // TODO to be added at a later stage
  // evalPrimaries(inputsIndices, clusterAnalysis);
  // evalParents(inputsIndices, clusterAnalysis);
 
  // TODO to be added at a later stage
  // Called last because it sets the global position of the cluster?
  // Do not call it when recalculating clusters out of standard reconstruction
  // if (!mJustCluster)
  //  evalLocal2TrackingCSTransform();
 
  return clusterAnalysis;
}
 
//____________________________________________________________________________
template <class InputType>
void ClusterFactory<InputType>::evalDispersion(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
{
  double d = 0., wtot = 0.;
  int nstat = 0;
 
  // Calculates the dispersion in cell units
  double etaMean = 0.0, phiMean = 0.0;
 
  // Calculate mean values
  for (auto iInput : inputsIndices) {
 
    if (clusterAnalysis.E() > 0 && mInputsContainer[iInput].getEnergy() > 0) {
      auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
      auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta);
 
      // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
      // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
      if (mSharedCluster && nSupMod % 2) {
        ieta += EMCAL_COLS;
      }
 
      auto etai = static_cast<double>(ieta);
      auto phii = static_cast<double>(iphi);
      double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
 
      if (w > 0.0) {
        phiMean += phii * w;
        etaMean += etai * w;
        wtot += w;
      }
    }
  }
 
  if (wtot > 0) {
    phiMean /= wtot;
    etaMean /= wtot;
  } else {
    LOG(error) << Form("Wrong weight %f\n", wtot);
  }
 
  // Calculate dispersion
  for (auto iInput : inputsIndices) {
 
    if (clusterAnalysis.E() > 0 && mInputsContainer[iInput].getEnergy() > 0) {
      auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
      auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta);
 
      // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
      // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
      if (mSharedCluster && nSupMod % 2) {
        ieta += EMCAL_COLS;
      }
 
      auto etai = static_cast<double>(ieta);
      auto phii = static_cast<double>(iphi);
      double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
 
      if (w > 0.0) {
        nstat++;
        d += w * ((etai - etaMean) * (etai - etaMean) + (phii - phiMean) * (phii - phiMean));
      }
    }
  }
 
  if (wtot > 0 && nstat > 1) {
    d /= wtot;
  } else {
    d = 0.;
  }
 
  clusterAnalysis.setDispersion(TMath::Sqrt(d));
}
 
//____________________________________________________________________________
template <class InputType>
void ClusterFactory<InputType>::evalLocalPosition(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
{
 
  int nstat = 0;
 
  double dist = tMaxInCm(double(clusterAnalysis.E()));
 
  std::array<double, 3> clXYZ = {0., 0., 0.}, clRmsXYZ = {0., 0., 0.}, xyzi{};
  double wtot = 0., w = 0.;
 
  for (auto iInput : inputsIndices) {
 
    try {
      mGeomPtr->RelPosCellInSModule(mInputsContainer[iInput].getTower(), dist).GetCoordinates(xyzi[0], xyzi[1], xyzi[2]);
    } catch (InvalidCellIDException& e) {
      LOG(error) << e.what();
      continue;
    }
 
    // 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
    if (mSharedCluster && mSuperModuleNumber != mGeomPtr->GetSuperModuleNumber(mInputsContainer[iInput].getTower())) {
      xyzi[1] *= -1;
    }
 
    if (mLogWeight > 0.0) {
      w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
    } else {
      w = mInputsContainer[iInput].getEnergy(); // just energy
    }
 
    if (w > 0.0) {
      wtot += w;
      nstat++;
 
      for (int i = 0; i < 3; i++) {
        clXYZ[i] += (w * xyzi[i]);
        clRmsXYZ[i] += (w * xyzi[i] * xyzi[i]);
      }
    } // w > 0
  } // dig loop
 
  //  cout << " wtot " << wtot << endl;
 
  if (wtot > 0) {
    //    xRMS   = TMath::Sqrt(x2m - xMean*xMean);
    for (int i = 0; i < 3; i++) {
      clXYZ[i] /= wtot;
 
      if (nstat > 1) {
        clRmsXYZ[i] /= (wtot * wtot);
        clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i] * clXYZ[i];
 
        if (clRmsXYZ[i] > 0.0) {
          clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
        } else {
          clRmsXYZ[i] = 0;
        }
      } else {
        clRmsXYZ[i] = 0;
      }
    }
  } else {
    for (int i = 0; i < 3; i++) {
      clXYZ[i] = clRmsXYZ[i] = -1.;
    }
  }
 
  clusterAnalysis.setLocalPosition(math_utils::Point3D<float>(clXYZ[0], clXYZ[1], clXYZ[2]));
}
 
//____________________________________________________________________________
template <class InputType>
void ClusterFactory<InputType>::evalGlobalPosition(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
{
 
  int i = 0, nstat = 0;
 
  double dist = tMaxInCm(double(clusterAnalysis.E()));
 
  std::array<double, 3> clXYZ = {0., 0., 0.}, clRmsXYZ = {0., 0., 0.}, lxyzi{}, xyzi{};
  double wtot = 0., w = 0.;
 
  for (auto iInput : inputsIndices) {
 
    // get the local coordinates of the cell
    try {
      mGeomPtr->RelPosCellInSModule(mInputsContainer[iInput].getTower(), dist).GetCoordinates(lxyzi[0], lxyzi[1], lxyzi[2]);
    } catch (InvalidCellIDException& e) {
      LOG(error) << e.what();
      continue;
    }
 
    // Now get the global coordinate
    mGeomPtr->GetGlobal(lxyzi, xyzi, mGeomPtr->GetSuperModuleNumber(mInputsContainer[iInput].getTower()));
 
    if (mLogWeight > 0.0) {
      w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
    } else {
      w = mInputsContainer[iInput].getEnergy(); // just energy
    }
 
    if (w > 0.0) {
      wtot += w;
      nstat++;
 
      for (i = 0; i < 3; i++) {
        clXYZ[i] += (w * xyzi[i]);
        clRmsXYZ[i] += (w * xyzi[i] * xyzi[i]);
      }
    }
  }
 
  //  cout << " wtot " << wtot << endl;
 
  if (wtot > 0) {
    //    xRMS   = TMath::Sqrt(x2m - xMean*xMean);
    for (i = 0; i < 3; i++) {
      clXYZ[i] /= wtot;
 
      if (nstat > 1) {
        clRmsXYZ[i] /= (wtot * wtot);
        clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i] * clXYZ[i];
 
        if (clRmsXYZ[i] > 0.0) {
          clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
        } else {
          clRmsXYZ[i] = 0;
        }
      } else {
        clRmsXYZ[i] = 0;
      }
    }
  } else {
    for (i = 0; i < 3; i++) {
      clXYZ[i] = clRmsXYZ[i] = -1.;
    }
  }
 
  clusterAnalysis.setGlobalPosition(math_utils::Point3D<float>(clXYZ[0], clXYZ[1], clXYZ[2]));
}
 
//____________________________________________________________________________
template <class InputType>
void ClusterFactory<InputType>::evalLocalPositionFit(double deff, double mLogWeight,
                                                     double phiSlope, std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
{
  int i = 0, nstat = 0;
  std::array<double, 3> clXYZ = {0., 0., 0.}, clRmsXYZ = {0., 0., 0.}, xyzi{};
  double wtot = 0., w = 0.;
 
  for (auto iInput : inputsIndices) {
 
    try {
      mGeomPtr->RelPosCellInSModule(mInputsContainer[iInput].getTower(), deff).GetCoordinates(xyzi[0], xyzi[1], xyzi[2]);
    } catch (InvalidCellIDException& e) {
      LOG(error) << e.what();
      continue;
    }
 
    if (mLogWeight > 0.0) {
      w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
    } else {
      w = mInputsContainer[iInput].getEnergy(); // just energy
    }
 
    if (w > 0.0) {
      wtot += w;
      nstat++;
 
      for (i = 0; i < 3; i++) {
        clXYZ[i] += (w * xyzi[i]);
        clRmsXYZ[i] += (w * xyzi[i] * xyzi[i]);
      }
    }
  } // loop
 
  //  cout << " wtot " << wtot << endl;
 
  if (wtot > 0) {
    //    xRMS   = TMath::Sqrt(x2m - xMean*xMean);
    for (i = 0; i < 3; i++) {
      clXYZ[i] /= wtot;
 
      if (nstat > 1) {
        clRmsXYZ[i] /= (wtot * wtot);
        clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i] * clXYZ[i];
 
        if (clRmsXYZ[i] > 0.0) {
          clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
        } else {
          clRmsXYZ[i] = 0;
        }
      } else {
        clRmsXYZ[i] = 0;
      }
    }
  } else {
    for (i = 0; i < 3; i++) {
      clXYZ[i] = clRmsXYZ[i] = -1.;
    }
  }
 
  // clRmsXYZ[i] ??
 
  if (phiSlope != 0.0 && mLogWeight > 0.0 && wtot != 0.0) {
    // Correction in phi direction (y - coords here); Aug 16;
    // May be put to global level or seperate method
    double ycorr = clXYZ[1] * (1. + phiSlope);
 
    // printf(" y %f : ycorr %f : slope %f \n", clXYZ[1], ycorr, phiSlope);
    clXYZ[1] = ycorr;
  }
 
  clusterAnalysis.setLocalPosition(math_utils::Point3D<float>(clXYZ[0], clXYZ[1], clXYZ[2]));
}
 
//_____________________________________________________________________________
template <class InputType>
void ClusterFactory<InputType>::getDeffW0(const double esum, double& deff, double& w0)
{
  double e = 0.0;
  const double kdp0 = 9.25147, kdp1 = 1.16700; // Hard coded now
  const double kwp0 = 4.83713, kwp1 = -2.77970e-01, kwp2 = 4.41116;
 
  // No extrapolation here
  e = esum < 0.5 ? 0.5 : esum;
  e = e > 100. ? 100. : e;
 
  deff = kdp0 + kdp1 * TMath::Log(e);
  w0 = kwp0 / (1. + TMath::Exp(kwp1 * (e + kwp2)));
}
 
//______________________________________________________________________________
template <class InputType>
void ClusterFactory<InputType>::evalCoreEnergy(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
{
 
  float coreEnergy = 0.;
 
  if (clusterAnalysis.getLocalPosition().Mag2() > 0.) {
    evalLocalPosition(inputsIndices, clusterAnalysis);
  }
 
  double phiPoint = clusterAnalysis.getLocalPosition().Phi();
  double etaPoint = clusterAnalysis.getLocalPosition().Eta();
  for (auto iInput : inputsIndices) {
 
    auto [eta, phi] = mGeomPtr->EtaPhiFromIndex(mInputsContainer[iInput].getTower());
    phi = phi * TMath::DegToRad();
 
    double distance = TMath::Sqrt((eta - etaPoint) * (eta - etaPoint) + (phi - phiPoint) * (phi - phiPoint));
 
    if (distance < mCoreRadius) {
      coreEnergy += mInputsContainer[iInput].getEnergy();
    }
  }
  clusterAnalysis.setCoreEnergy(coreEnergy);
}
 
//____________________________________________________________________________
template <class InputType>
void ClusterFactory<InputType>::evalNExMax(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
{
  // Pre-compute cell indices and energies for all cells in cluster to avoid multiple expensive geometry lookups
  const size_t n = inputsIndices.size();
  std::vector<short> rows;
  std::vector<short> columns;
  std::vector<double> energies;
 
  rows.reserve(n);
  columns.reserve(n);
  energies.reserve(n);
 
  for (auto iInput : inputsIndices) {
    auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
 
    // get a nice topological indexing that is done in exactly the same way as used by the clusterizer
    // this way we can handle the shared cluster cases correctly
    const auto [row, column] = mGeomPtr->GetTopologicalRowColumn(nSupMod, nModule, nIphi, nIeta);
 
    rows.push_back(row);
    columns.push_back(column);
    energies.push_back(mInputsContainer[iInput].getEnergy());
  }
 
  // Now find local maxima using pre-computed data
  int nExMax = 0;
  for (size_t i = 0; i < n; i++) {
    // this cell is assumed to be local maximum unless we find a higher energy cell in the neighborhood
    bool isExMax = true;
 
    // loop over all other cells in cluster
    for (size_t j = 0; j < n; j++) {
      if (i == j) {
        continue;
      }
 
      // adjacent cell is any cell with adjacent phi or eta index
      if (std::abs(rows[i] - rows[j]) <= 1 &&
          std::abs(columns[i] - columns[j]) <= 1) {
 
        // if there is a cell with higher energy than the current cell, it is not a local maximum
        if (energies[j] > energies[i]) {
          isExMax = false;
          break;
        }
      }
    }
    if (isExMax) {
      nExMax++;
    }
  }
  clusterAnalysis.setNExMax(nExMax);
}
 
//____________________________________________________________________________
template <class InputType>
void ClusterFactory<InputType>::evalElipsAxis(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
{
  double wtot = 0.;
  double x = 0.;
  double z = 0.;
  double dxx = 0.;
  double dzz = 0.;
  double dxz = 0.;
 
  std::array<float, 2> lambda{};
 
  for (auto iInput : inputsIndices) {
 
    auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
    auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta);
 
    // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
    // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
    if (mSharedCluster && nSupMod % 2) {
      ieta += EMCAL_COLS;
    }
 
    auto etai = static_cast<double>(ieta);
    auto phii = static_cast<double>(iphi);
 
    double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
    // clusterAnalysis.E() summed amplitude of inputs, i.e. energy of cluster
    // Gives smaller value of lambda than log weight
    // w = mEnergyList[iInput] / clusterAnalysis.E(); // Nov 16, 2006 - try just energy
 
    dxx += w * etai * etai;
    x += w * etai;
    dzz += w * phii * phii;
    z += w * phii;
 
    dxz += w * etai * phii;
 
    wtot += w;
  }
 
  if (wtot > 0) {
    dxx /= wtot;
    x /= wtot;
    dxx -= x * x;
    dzz /= wtot;
    z /= wtot;
    dzz -= z * z;
    dxz /= wtot;
    dxz -= x * z;
 
    lambda[0] = 0.5 * (dxx + dzz) + TMath::Sqrt(0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz);
 
    if (lambda[0] > 0) {
      lambda[0] = TMath::Sqrt(lambda[0]);
    } else {
      lambda[0] = 0;
    }
 
    lambda[1] = 0.5 * (dxx + dzz) - TMath::Sqrt(0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz);
 
    if (lambda[1] > 0) { // To avoid exception if numerical errors lead to negative lambda.
      lambda[1] = TMath::Sqrt(lambda[1]);
    } else {
      lambda[1] = 0.;
    }
  } else {
    lambda[0] = 0.;
    lambda[1] = 0.;
  }
 
  clusterAnalysis.setM02(lambda[0] * lambda[0]);
  clusterAnalysis.setM20(lambda[1] * lambda[1]);
}
 
//____________________________________________________________________________
template <class InputType>
std::tuple<int, float, float, bool> ClusterFactory<InputType>::getMaximalEnergyIndex(std::span<const int> inputsIndices) const
{
 
  float energy = 0.;
  int mid = 0;
  float cellAmp = 0;
  int iSupMod0 = -1;
  bool shared = false;
  for (auto iInput : inputsIndices) {
    if (iInput >= mInputsContainer.size()) {
      throw CellIndexRangeException(iInput, mInputsContainer.size());
    }
    cellAmp += mInputsContainer[iInput].getEnergy();
    if (iSupMod0 == -1) {
      iSupMod0 = mGeomPtr->GetSuperModuleNumber(mInputsContainer[iInput].getTower());
    } else if (iSupMod0 != mGeomPtr->GetSuperModuleNumber(mInputsContainer[iInput].getTower())) {
      shared = true;
    }
    if (mInputsContainer[iInput].getEnergy() > energy) {
      energy = mInputsContainer[iInput].getEnergy();
      mid = iInput;
    }
  } // loop on cluster inputs
 
  return std::make_tuple(mid, energy, cellAmp, shared);
}
 
//____________________________________________________________________________
template <class InputType>
bool ClusterFactory<InputType>::isExoticCell(short towerId, float ecell, float const exoticTime, float& fCross) const
{
  if (ecell < mExoticCellMinAmplitude) {
    return false; // do not reject low energy cells
  }
 
  // if the look up table is not set yet (mostly due to a reset call) then set it up now.
  if (!getLookUpInit()) {
    throw UninitLookUpTableException();
  }
 
  float eCross = getECross(towerId, ecell, exoticTime);
  fCross = 1.f - eCross / ecell;
 
  if (fCross > mExoticCellFraction) {
    LOG(debug) << "EXOTIC CELL id " << towerId << ", eCell " << ecell << ", eCross " << eCross << ", 1-eCross/eCell " << 1 - eCross / ecell;
    return true;
  }
 
  return false;
}
 
//____________________________________________________________________________
template <class InputType>
float ClusterFactory<InputType>::getECross(short absID, float energy, float const exoticTime) const
{
  auto [iSM, iMod, iIphi, iIeta] = mGeomPtr->GetCellIndex(absID);
  auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(iSM, iMod, iIphi, iIeta);
 
  // Get close cells index, energy and time, not in corners
 
  short towerId1 = -1;
  short towerId2 = -1;
 
  if (iphi < o2::emcal::EMCAL_ROWS - 1) {
    try {
      towerId1 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi + 1, ieta);
    } catch (InvalidCellIDException& e) {
      towerId1 = -1 * e.getCellID();
    }
  }
  if (iphi > 0) {
    try {
      towerId2 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi - 1, ieta);
    } catch (InvalidCellIDException& e) {
      towerId2 = -1 * e.getCellID();
    }
  }
 
  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
 
  short towerId3 = -1;
  short towerId4 = -1;
 
  if (ieta == o2::emcal::EMCAL_COLS - 1 && (iSM % 2) == 0) {
    try {
      towerId3 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM + 1, iphi, 0);
    } catch (InvalidCellIDException& e) {
      towerId3 = -1 * e.getCellID();
    }
    try {
      towerId4 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta - 1);
    } catch (InvalidCellIDException& e) {
      towerId4 = -1 * e.getCellID();
    }
  } else if (ieta == 0 && (iSM % 2) != 0) {
    try {
      towerId3 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta + 1);
    } catch (InvalidCellIDException& e) {
      towerId3 = -1 * e.getCellID();
    }
    try {
      towerId4 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM - 1, iphi, o2::emcal::EMCAL_COLS - 1);
    } catch (InvalidCellIDException& e) {
      towerId4 = -1 * e.getCellID();
    }
  } else {
    if (ieta < o2::emcal::EMCAL_COLS - 1) {
      try {
        towerId3 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta + 1);
      } catch (InvalidCellIDException& e) {
        towerId3 = -1 * e.getCellID();
      }
    }
    if (ieta > 0) {
      try {
        towerId4 = mGeomPtr->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta - 1);
      } catch (InvalidCellIDException& e) {
        towerId4 = -1 * e.getCellID();
      }
    }
  }
 
  LOG(debug) << "iSM " << iSM << ", absID " << absID << ", a " << towerId1 << ", b " << towerId2 << ", c " << towerId3 << ", e " << towerId3;
 
  short index1 = (towerId1 > -1) ? mLoolUpTowerToIndex.at(towerId1) : -1;
  short index2 = (towerId2 > -1) ? mLoolUpTowerToIndex.at(towerId2) : -1;
  short index3 = (towerId3 > -1) ? mLoolUpTowerToIndex.at(towerId3) : -1;
  short index4 = (towerId4 > -1) ? mLoolUpTowerToIndex.at(towerId4) : -1;
 
  std::array<std::pair<float, float>, 4> cellData = {
    {{(index1 > -1) ? mInputsContainer[index1].getEnergy() : 0., (index1 > -1) ? mInputsContainer[index1].getTimeStamp() : 0.},
     {(index2 > -1) ? mInputsContainer[index2].getEnergy() : 0., (index2 > -1) ? mInputsContainer[index2].getTimeStamp() : 0.},
     {(index3 > -1) ? mInputsContainer[index3].getEnergy() : 0., (index3 > -1) ? mInputsContainer[index3].getTimeStamp() : 0.},
     {(index4 > -1) ? mInputsContainer[index4].getEnergy() : 0., (index4 > -1) ? mInputsContainer[index4].getTimeStamp() : 0.}}};
 
  for (auto& cell : cellData) {
    if (std::abs(exoticTime - cell.second) > mExoticCellDiffTime) {
      cell.first = 0;
    }
  }
 
  float w1 = 1, w2 = 1, w3 = 1, w4 = 1;
  if (mUseWeightExotic) {
    w1 = GetCellWeight(cellData[0].first, energy);
    w2 = GetCellWeight(cellData[1].first, energy);
    w3 = GetCellWeight(cellData[2].first, energy);
    w4 = GetCellWeight(cellData[3].first, energy);
  }
 
  if (cellData[0].first < mExoticCellInCrossMinAmplitude || w1 <= 0) {
    cellData[0].first = 0;
  }
  if (cellData[1].first < mExoticCellInCrossMinAmplitude || w2 <= 0) {
    cellData[1].first = 0;
  }
  if (cellData[2].first < mExoticCellInCrossMinAmplitude || w3 <= 0) {
    cellData[2].first = 0;
  }
  if (cellData[3].first < mExoticCellInCrossMinAmplitude || w4 <= 0) {
    cellData[3].first = 0;
  }
 
  return cellData[0].first + cellData[1].first + cellData[2].first + cellData[3].first;
}
 
//____________________________________________________________________________
template <class InputType>
float ClusterFactory<InputType>::GetCellWeight(float eCell, float eCluster) const
{
  if (eCell > 0 && eCluster > 0) {
    if (mLogWeight > 0) {
      return std::max(0.f, mLogWeight + std::log(eCell / eCluster));
    }
    return std::log(eCluster / eCell);
  }
  return 0.;
}
 
//____________________________________________________________________________
template <class InputType>
int ClusterFactory<InputType>::getMultiplicityAtLevel(float level, std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
{
  int multipl = 0;
  for (auto iInput : inputsIndices) {
    if (mInputsContainer[iInput].getEnergy() > level * clusterAnalysis.E()) {
      multipl++;
    }
  }
 
  return multipl;
}
 
//____________________________________________________________________________
template <class InputType>
void ClusterFactory<InputType>::evalTime(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
{
  float maxE = 0;
  unsigned short maxAt = 0;
  for (auto iInput : inputsIndices) {
    if (mInputsContainer[iInput].getEnergy() > maxE) {
      maxE = mInputsContainer[iInput].getEnergy();
      maxAt = iInput;
    }
  }
 
  clusterAnalysis.setClusterTime(mInputsContainer[maxAt].getTimeStamp());
}
 
//_____________________________________________________________________
template <class InputType>
double ClusterFactory<InputType>::tMaxInCm(const double e, const int key) const
{
  const double ca = 4.82; // shower max parameter - first guess; ca=TMath::Log(1000./8.07)
  double tmax = 0.;       // position of electromagnetic shower max in cm
 
  const double x0 = 1.31; // radiation lenght (cm)
 
  if (e > 0.1) {
    tmax = TMath::Log(e) + ca;
    if (key == 0) {
      tmax += 0.5;
    } else {
      tmax -= 0.5;
    }
    tmax *= x0; // convert to cm
  }
 
  return tmax;
}
 
//______________________________________________________________________________
template <class InputType>
float ClusterFactory<InputType>::etaToTheta(float arg) const
{
  return (2. * TMath::ATan(TMath::Exp(-arg)));
}
 
//______________________________________________________________________________
template <class InputType>
float ClusterFactory<InputType>::thetaToEta(float arg) const
{
  return (-1 * TMath::Log(TMath::Tan(0.5 * arg)));
}
 
template <class InputType>
ClusterFactory<InputType>::ClusterIterator::ClusterIterator(const ClusterFactory& factory, int clusterIndex, bool forward) : mClusterFactory(&factory),
                                                                                                                             mCurrentCluster(mClusterFactory->buildCluster(clusterIndex)),
                                                                                                                             mClusterID(clusterIndex),
                                                                                                                             mForward(forward)
{
}
 
template <class InputType>
bool ClusterFactory<InputType>::ClusterIterator::operator==(const ClusterFactory<InputType>::ClusterIterator& rhs) const
{
  return mClusterFactory == rhs.mClusterFactory && mClusterID == rhs.mClusterID && mForward == rhs.mForward;
}
 
template <class InputType>
typename ClusterFactory<InputType>::ClusterIterator& ClusterFactory<InputType>::ClusterIterator::operator++()
{
  if (mForward) {
    mClusterID++;
  } else {
    mClusterID--;
  }
  mCurrentCluster = mClusterFactory->buildCluster(mClusterID);
  return *this;
}
 
template <class InputType>
typename ClusterFactory<InputType>::ClusterIterator ClusterFactory<InputType>::ClusterIterator::operator++(int)
{
  auto tmp = *this;
  ++(*this);
  return tmp;
}
 
template <class InputType>
typename ClusterFactory<InputType>::ClusterIterator& ClusterFactory<InputType>::ClusterIterator::operator--()
{
  if (mForward) {
    mClusterID--;
  } else {
    mClusterID++;
  }
  mCurrentCluster = mClusterFactory->buildCluster(mClusterID);
  return *this;
}
 
template <class InputType>
typename ClusterFactory<InputType>::ClusterIterator ClusterFactory<InputType>::ClusterIterator::operator--(int)
{
  auto tmp = *this;
  --(*this);
  return tmp;
}
 
template class o2::emcal::ClusterFactory<o2::emcal::Cell>;
template class o2::emcal::ClusterFactory<o2::emcal::Digit>;