CaloAxis.h

Go to the documentation of this file.

/*
 * CaloAxis.h
 *
 *  Created on: 27 April 2020
 *      Author: Valerio Vagelli
 */

#ifndef HERD_CALOAXIS_H_
#define HERD_CALOAXIS_H_

// HerdSoftware headersCalo
#include "common/DefaultValues.h"
#include "common/Vec3D.h"
#include "dataobjects/CaloGeoParams.h"
#include "dataobjects/CaloHits.h"
#include "dataobjects/Hit.h"
#include "dataobjects/Point.h"

#ifdef HS_USE_ROOT
#include "Rtypes.h"
#endif

// C/C++ standard headers
#include <algorithm>
#include <array>
#include <map>
#include <numeric>
#include <vector>

namespace Herd {

enum class CaloAxisAlgoType {
  NONE, 
  PCA   
};

class CaloAxis {
public:
  CaloAxis()
      : DEFAULT_INIT(m_axisCOG), DEFAULT_INIT(m_axisEigenvalues), DEFAULT_INIT(m_axisEigenvectors),
        DEFAULT_INIT(m_axisDir), DEFAULT_INIT(m_axisSigma), DEFAULT_INIT(m_axisSkewness), DEFAULT_INIT(m_axisKurtosis),
        DEFAULT_INIT(m_axisEntryHit), DEFAULT_INIT(m_axisExitHit), DEFAULT_INIT(m_axisPathlengthHit){};

  CaloAxis(const CaloHits caloHits, float energyThreshold)
      : DEFAULT_INIT(m_axisCOG), DEFAULT_INIT(m_axisEigenvalues), DEFAULT_INIT(m_axisEigenvectors),
        DEFAULT_INIT(m_axisDir), DEFAULT_INIT(m_axisSigma), DEFAULT_INIT(m_axisSkewness), DEFAULT_INIT(m_axisKurtosis),
        DEFAULT_INIT(m_axisEntryHit), DEFAULT_INIT(m_axisExitHit),
        DEFAULT_INIT(m_axisPathlengthHit), m_caloHits{std::move(caloHits)}, m_energyThreshold{energyThreshold} {};

  virtual ~CaloAxis() = default;

  std::array<double, 3> AxisSigma(const CaloGeoParams &geoParams) {
    if (m_axisSigma == DefaultValue(m_axisSigma)) {
      ComputeMomenta(geoParams);
    }
    return m_axisSigma;
  };

  std::array<double, 3> AxisSkewness(const CaloGeoParams &geoParams) {
    if (m_axisSkewness == DefaultValue(m_axisSkewness)) {
      ComputeMomenta(geoParams);
    }
    return m_axisSkewness;
  };

  std::array<double, 3> AxisKurtosis(const CaloGeoParams &geoParams) {
    if (m_axisKurtosis == DefaultValue(m_axisKurtosis)) {
      ComputeMomenta(geoParams);
    }
    return m_axisKurtosis;
  };

  CaloAxisAlgoType AlgoType() const { return m_axisAlgo; };

  unsigned short AxisHits() const {
    return std::count_if(begin(m_caloHits), end(m_caloHits),
                         [this](const Hit &hit) { return hit.EDep() > m_energyThreshold; });
  };

  Point COG() const { return m_axisCOG; };

  std::array<double, 3> Eigenvalues() const { return m_axisEigenvalues; };

  std::array<Point, 3> Eigenvectors() const { return m_axisEigenvectors; };

  Line Axis() const { return Line(m_axisCOG, m_axisDir); };

  Point AxisEntryHit(const CaloGeoParams &geoParams) {
    if (m_axisEntryHit == DefaultValue(m_axisEntryHit)) {
      ComputePathlength(geoParams);
    }
    return m_axisEntryHit;
  }

  Point AxisExitHit(const CaloGeoParams &geoParams) {
    if (m_axisExitHit == DefaultValue(m_axisExitHit)) {
      ComputePathlength(geoParams);
    }
    return m_axisExitHit;
  }

  double AxisPathlengthHit(const CaloGeoParams &geoParams) {
    if (m_axisPathlengthHit == DefaultValue(m_axisPathlengthHit)) {
      ComputePathlength(geoParams);
    }
    return m_axisPathlengthHit;
  }

  double LateralEnergy(const CaloGeoParams &geoParams, double radmin = 0, double radmax = defaultValue<double>);

  double LongitudinalEntryEnergy(const CaloGeoParams &geoParams, double dist, double radmin = 0,
                                 double radmax = defaultValue<double>);

  double LongitudinalExitEnergy(const CaloGeoParams &geoParams, double dist, double radmin = 0,
                                double radmax = defaultValue<double>);

  double LocalEnergy(const CaloGeoParams &geoParams, Point center, double distneg, double distpos, double radmin = 0,
                     double radmax = defaultValue<double>);

  void SetAlgoType(CaloAxisAlgoType algoType) { m_axisAlgo = algoType; };

  void SetCOG(Point cog) { m_axisCOG = cog; };

  void SetDirection(Vec3D dir) { m_axisDir = dir; };

  void SetEigenvalues(std::array<double, 3> eig) { m_axisEigenvalues = eig; };

  void SetEigenvectors(std::array<Point, 3> eig) { m_axisEigenvectors = eig; };

  void Reset();

private:
  CaloAxisAlgoType m_axisAlgo;
  Point m_axisCOG;
  std::array<double, 3> m_axisEigenvalues;
  std::array<Point, 3> m_axisEigenvectors;
  Vec3D m_axisDir;
  std::array<double, 3> m_axisSigma;
  std::array<double, 3> m_axisSkewness;
  std::array<double, 3> m_axisKurtosis;
  Point m_axisEntryHit;
  Point m_axisExitHit;
  double m_axisPathlengthHit;

  CaloHits m_caloHits; // Internal reference to caloHits used for this Axis
  void ComputeMomenta(const CaloGeoParams &geoParams);

  void ComputePathlength(const CaloGeoParams &geoParams);

  float m_energyThreshold = 0;

#ifdef HS_USE_ROOT
  ClassDef(CaloAxis, 1)
#endif
};

} // namespace Herd

#endif /* HERD_CALOAXIS_H_ */