Rich.cpp

Go to the documentation of this file.

#include "Containers/Rich.h"
 
namespace NAIA {
 
const std::string RichBase::BranchName = "RichBaseData";
const std::string RichPlus::BranchName = "RichPlusData";
 
void RichBaseData::Clear() { m_beta.clear(); }
void RichBaseData::Dump() const {}
 
bool RichBaseData::HasBetaLIP() const {
  if (!m_beta.empty()) {
    return m_beta[0].hasBetaLIP;
  }
 
  return false;
}
 
bool RichBaseData::IsNaF() const {
  if (!m_beta.empty()) {
    return m_beta[0].isNaF;
  }
 
  return false;
}
 
RichBetaVariable<float> RichBaseData::GetBeta() const {
  if (!m_beta.empty()) {
    return m_beta[0].beta;
  }
 
  return {};
}
 
RichBetaVariable<float> RichBaseData::GetBetaError() const {
  if (!m_beta.empty()) {
    return m_beta[0].ebeta;
  }
 
  return {};
}
 
std::array<float, 2> RichBaseData::GetTrackEmissionPoint() const {
  if (!m_beta.empty()) {
    return m_beta[0].Pos;
  }
 
  return {};
}
 
double RichBaseData::BetaConsistency() const {
  if (!m_beta.empty() && m_beta[0].beta.size() > 1) {
    return std::fabs(m_beta[0].beta.at(Rich::BetaType::CIEMAT) - m_beta[0].beta.at(Rich::BetaType::LIP));
  }
 
  return std::numeric_limits<double>::max();
}
 
double RichBaseData::BetaConsistencyTof(const TofBaseData &tofData) const {
  if (!m_beta.empty() && ContainsKeys(m_beta[0].beta, Rich::BetaType::CIEMAT) &&
      ContainsKeys(tofData.Beta, Tof::BetaType::BetaH)) {
    return std::fabs(m_beta[0].beta.at(Rich::BetaType::CIEMAT) - tofData.Beta.at(Tof::BetaType::BetaH));
  }
 
  return std::numeric_limits<double>::max();
}
 
int RichBaseData::GetTileIndex() const {
  if (m_beta.empty())
    return -1;
  float x = m_beta[0].Pos[0];
  float y = m_beta[0].Pos[1];
  constexpr int grid_side_length = 11;
  constexpr int n_tiles = grid_side_length * grid_side_length;
  constexpr double tile_width = 0.1 + 11.5;
 
  if (std::fabs(x) < tile_width * 1.5 && std::fabs(y) < tile_width * 1.5)
    return 121;
 
  int nx = int(x / tile_width + 5.5);
  int ny = int(y / tile_width + 5.5);
  int t = ny * grid_side_length + nx;
  return t > n_tiles ? -1 : t;
}
 
bool RichBaseData::RingGeomTest(RingGeom geom) const {
 
  if (m_beta.empty())
    return false;
 
  float x = m_beta[0].Pos[0];
  float y = m_beta[0].Pos[1];
 
  static constexpr int nbadtiles = 7; // tiles with bad beta reconstruction
  static constexpr std::array<int, nbadtiles> kbadtile = {3, 7, 12, 20, 87, 100, 108};
 
  // mdb: Cut particles that are on the outer edge of the aerogel, esp if they fall on small tiles (?)
  //  aerogel external border (r**2)
  static constexpr float cut_aerogelexternalborder = 3422.25;
  static constexpr std::array<float, 2> cut_aerogelborder = {28.5, 29.5};
 
  // mdb: Also cut particles that fall close to the NaF/aerogel border
  //  aerogel/NaF border ([0]=NaF, [1]=aerogel)
  static constexpr std::array<float, 2> cut_aerogel_nafborder = {17., 19.};
 
  // clang-format off
  // Bad aerogel regions (from TWiki)
  static constexpr int npos = 9;
  static constexpr float xpos[npos] = { -2.,  20., -29., 30, 16, 30., -30., -10., 12.};
  static constexpr float xw[npos]   = {  2.,   2.,   1., 2., 2.,  2.,   2.,   2.,  2.};
  static constexpr float ypos[npos] = {-60., -42., -18., 16, 22, 38.,  40.,  56., 56.};
  static constexpr float yw[npos]   = {  2.,   2.,   2., 2., 2., 1.5,   2.,  1.5,  2.};
  // clang-format on
 
  if (geom == RingGeom::Borders) {
    if (IsNaF()) {
      if (std::max(std::fabs(x), std::fabs(y)) > cut_aerogel_nafborder[0])
        return false;
    } else {
      if (std::max(std::fabs(x), std::fabs(y)) < cut_aerogel_nafborder[1])
        return false;
      if (x * x + y * y > cut_aerogelexternalborder)
        return false;
    }
  } else if (geom == RingGeom::BadAglTiles) {
    for (int kbad = 0; kbad < nbadtiles; kbad++) {
      if (GetTileIndex() == kbadtile[kbad])
        return false;
    }
  } else if (geom == RingGeom::AglTilesEdge) {
    if (std::max(std::fabs(x), std::fabs(y)) > cut_aerogelborder[0] &&
        std::max(std::fabs(x), std::fabs(y)) < cut_aerogelborder[1])
      return false;
  } else if (geom == RingGeom::BadAglRegions) {
    for (int i = 0; i < npos; i++) {
      bool isx = (x > (xpos[i]) && x < (xpos[i] + xw[i]));
      bool isy = (y > (ypos[i]) && y < (ypos[i] + yw[i]));
      if (isx && isy)
        return false;
    }
  }
  return true;
}
 
//=======================================
void RichPlusData::Clear() { m_data.clear(); }
void RichPlusData::Dump() const {}
 
} // namespace NAIA