ClusterFactory.h

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.
#ifndef DETECTORS_EMCAL_BASE_INCLUDE_EMCALBASE_CLUSTERFACTORY_H_
#define DETECTORS_EMCAL_BASE_INCLUDE_EMCALBASE_CLUSTERFACTORY_H_
 
#include "EMCALBase/Geometry.h"
#include "DataFormatsEMCAL/AnalysisCluster.h"
#include "DataFormatsEMCAL/Cell.h" // IWYU pragma: keep
#include "DataFormatsEMCAL/CellLabel.h"
#include "DataFormatsEMCAL/Cluster.h"
#include "DataFormatsEMCAL/ClusterLabel.h"
#include "DataFormatsEMCAL/Digit.h" // IWYU pragma: keep
// #include "MathUtils/Cartesian.h"
 
#include <Rtypes.h>
 
#include <fmt/format.h>
 
#include <array>
#include <span>
#include <string>
#include <tuple>
#include <vector>
 
namespace o2::emcal
{
 
template <class InputType>
class ClusterFactory
{
 
 public:
  class ClusterRangeException final : public std::exception
  {
   public:
    ClusterRangeException(int clusterIndex, int maxClusters) : std::exception(),
                                                               mClusterID(clusterIndex),
                                                               mMaxClusters(maxClusters),
                                                               mErrorMessage(fmt::format("Cluster out of range: %d, max %d", mClusterID, mMaxClusters))
    {
    }
 
    ~ClusterRangeException() noexcept final = default;
 
    [[nodiscard]] const char* what() const noexcept final { return mErrorMessage.data(); }
 
    [[nodiscard]] int getClusterID() const { return mClusterID; }
 
    [[nodiscard]] int getMaxNumberOfClusters() const { return mMaxClusters; }
 
   private:
    int mClusterID = 0;        
    int mMaxClusters = 0;      
    std::string mErrorMessage; 
  };
 
  class CellIndexRangeException final : public std::exception
  {
   public:
    CellIndexRangeException(int cellIndex, int maxCellIndex) : std::exception(),
                                                               mCellIndex(cellIndex),
                                                               mMaxCellIndex(maxCellIndex),
                                                               mErrorMessage(fmt::format("Cell Index out of range: %d, max %d", mCellIndex, mMaxCellIndex))
    {
    }
 
    ~CellIndexRangeException() noexcept final = default;
 
    [[nodiscard]] const char* what() const noexcept final { return mErrorMessage.data(); }
 
    [[nodiscard]] int getCellIndex() const { return mCellIndex; }
 
    [[nodiscard]] int getMaxNumberOfCellIndexs() const { return mMaxCellIndex; }
 
   private:
    int mCellIndex = 0;        
    int mMaxCellIndex = 0;     
    std::string mErrorMessage; 
  };
 
  class GeometryNotSetException final : public std::exception
  {
   public:
    GeometryNotSetException() = default;
    ~GeometryNotSetException() noexcept final = default;
 
    [[nodiscard]] const char* what() const noexcept final { return "Geometry not set"; }
  };
 
  class ClusterIterator
  {
   public:
    ClusterIterator(const ClusterFactory& factory, int clusterIndex, bool forward);
 
    ~ClusterIterator() = default;
 
    bool operator==(const ClusterIterator& rhs) const;
 
    bool operator!=(const ClusterIterator& rhs) const { return !(*this == rhs); }
 
    ClusterIterator& operator++();
 
    ClusterIterator operator++(int);
 
    ClusterIterator& operator--();
 
    ClusterIterator operator--(int);
 
    AnalysisCluster* operator*() { return &mCurrentCluster; }
 
    [[nodiscard]] int current_index() const { return mClusterID; }
 
   private:
    const ClusterFactory* mClusterFactory; 
    AnalysisCluster mCurrentCluster;       
    int mClusterID = 0;                    
    bool mForward = true;                  
  };
 
  ClusterFactory() = default;
 
  ClusterFactory(std::span<const o2::emcal::Cluster> clustersContainer, std::span<const InputType> inputsContainer, std::span<const int> cellsIndices);
 
  ClusterFactory(const ClusterFactory& rp) = default;
 
  ClusterFactory& operator=(const ClusterFactory& cf) = default;
 
  ~ClusterFactory() = default;
 
  ClusterIterator begin() const { return ClusterIterator(*this, 0, true); }
 
  ClusterIterator end() const { return ClusterIterator(*this, getNumberOfClusters(), true); }
 
  ClusterIterator rbegin() const { return ClusterIterator(*this, getNumberOfClusters() - 1, false); }
 
  ClusterIterator rend() const { return ClusterIterator(*this, -1, false); }
 
  void reset();
 
  AnalysisCluster buildCluster(int index, o2::emcal::ClusterLabel* clusterLabel = nullptr) const;
 
  void SetECALogWeight(Float_t w) { mLogWeight = w; }
  float GetECALogWeight() const { return mLogWeight; }
 
  void doEvalLocal2tracking(bool justCluster)
  {
    mJustCluster = justCluster;
  }
 
  void evalLocalPosition(std::span<const int> inputsIndices, AnalysisCluster& cluster) const;
 
  void evalGlobalPosition(std::span<const int> inputsIndices, AnalysisCluster& cluster) const;
 
  void evalLocal2TrackingCSTransform() const;
 
  void evalLocalPositionFit(double deff, double mLogWeight, double phiSlope, std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const;
 
  static void getDeffW0(const Double_t esum, Double_t& deff, Double_t& w0);
 
  std::tuple<int, float, float, bool> getMaximalEnergyIndex(std::span<const int> inputsIndices) const;
 
  bool isExoticCell(int16_t towerId, float ecell, float const exoticTime, float& fCross) const;
 
  float getECross(int16_t absID, float energy, float const exoticTime) const;
 
  float GetCellWeight(float eCell, float eCluster) const;
 
  int getMultiplicityAtLevel(float level, std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const;
 
  int getSuperModuleNumber() const { return mSuperModuleNumber; }
 
  // searches for the local maxima
  // energy above relative level
  // int getNumberOfLocalMax(int nInputMult,
  //                        float locMaxCut, std::span<InputType> inputs) const;
 
  // int getNumberOfLocalMax(std::vector<InputType>& maxAt, std::vector<float>& maxAtEnergy,
  //                         float locMaxCut, std::span<InputType> inputs) const;
 
  bool sharedCluster() const { return mSharedCluster; }
  void setSharedCluster(bool s) { mSharedCluster = s; }
 
  Double_t tMaxInCm(const Double_t e = 0.0, const int key = 0) const;
 
  bool getLookUpInit() const { return mLookUpInit; }
 
  float getCoreRadius() const { return mCoreRadius; }
  void setCoreRadius(float radius) { mCoreRadius = radius; }
 
  float getExoticCellFraction() const { return mExoticCellFraction; }
  void setExoticCellFraction(float exoticCellFraction) { mExoticCellFraction = exoticCellFraction; }
 
  float getExoticCellDiffTime() const { return mExoticCellDiffTime; }
  void setExoticCellDiffTime(float exoticCellDiffTime) { mExoticCellDiffTime = exoticCellDiffTime; }
 
  float getExoticCellMinAmplitude() const { return mExoticCellMinAmplitude; }
  void setExoticCellMinAmplitude(float exoticCellMinAmplitude) { mExoticCellMinAmplitude = exoticCellMinAmplitude; }
 
  float getExoticCellInCrossMinAmplitude() const { return mExoticCellInCrossMinAmplitude; }
  void setExoticCellInCrossMinAmplitude(float exoticCellInCrossMinAmplitude) { mExoticCellInCrossMinAmplitude = exoticCellInCrossMinAmplitude; }
 
  bool getUseWeightExotic() const { return mUseWeightExotic; }
  void setUseWeightExotic(bool useWeightExotic) { mUseWeightExotic = useWeightExotic; }
 
  void setContainer(std::span<const o2::emcal::Cluster> clusterContainer, std::span<const InputType> cellContainer, std::span<const int> indicesContainer, std::span<const o2::emcal::CellLabel> cellLabelContainer = {})
  {
    mClustersContainer = clusterContainer;
    mInputsContainer = cellContainer;
    mCellsIndices = indicesContainer;
    if (!getLookUpInit()) {
      setLookUpTable();
    }
    if (!cellLabelContainer.empty()) {
      mCellLabelContainer = cellLabelContainer;
    }
  }
 
  void setLookUpTable()
  {
    mLoolUpTowerToIndex.fill(-1);
    for (auto iCellIndex : mCellsIndices) {
      mLoolUpTowerToIndex[mInputsContainer[iCellIndex].getTower()] = iCellIndex;
    }
    mLookUpInit = true;
  }
 
  int getNumberOfClusters() const
  {
    return mClustersContainer.size();
  }
 
  void setGeometry(o2::emcal::Geometry* geometry) { mGeomPtr = geometry; }
 
  class UninitLookUpTableException final : public std::exception
  {
   public:
    UninitLookUpTableException() = default;
 
    ~UninitLookUpTableException() noexcept final = default;
 
    [[nodiscard]] const char* what() const noexcept final { return "Lookup table not initialized, exotics evaluation not possible!"; }
  };
 
 protected:
  void evalCoreEnergy(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const;
 
  void evalDispersion(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const;
 
  void evalElipsAxis(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const;
 
  void evalNExMax(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const;
 
  void evalTime(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const;
 
  float thetaToEta(float theta) const;
 
  float etaToTheta(float eta) const;
 
 private:
  struct CellGeomInfo {
    int8_t nSupMod;
    int8_t iphi;
    int8_t ieta;
    int8_t ietaShared;
    int16_t row;
    int16_t col;
  };
 
  o2::emcal::Geometry* mGeomPtr = nullptr;
 
  float mCoreRadius = 10; 
 
  float mLogWeight = 4.5; 
 
  bool mJustCluster = kFALSE; 
  bool mLookUpInit = false;   
 
  mutable int mSuperModuleNumber = 0;         
  float mDistToBadTower = -1;                 
  mutable bool mSharedCluster = false;        
  float mExoticCellFraction = 0.97;           
  float mExoticCellDiffTime = 1e6;            
  float mExoticCellMinAmplitude = 4.;         
  float mExoticCellInCrossMinAmplitude = 0.1; 
  bool mUseWeightExotic = false;              
 
  mutable std::vector<CellGeomInfo> mCellGeomBuffer;         
  std::span<const o2::emcal::Cluster> mClustersContainer;    
  std::span<const InputType> mInputsContainer;               
  std::span<const int> mCellsIndices;                        
  std::array<int16_t, 17664> mLoolUpTowerToIndex{};          
  std::span<const o2::emcal::CellLabel> mCellLabelContainer; 
 
  ClassDefNV(ClusterFactory, 3);
};
 
} // namespace o2::emcal
#endif // DETECTORS_EMCAL_BASE_INCLUDE_EMCALBASE_CLUSTERFACTORY_H_