packages/tpc2019/TpcPrototypeCluster.cc

0001 /**
0002  * @file trackbase/TpcPrototypeCluster.cc
0003  * @author D. McGlinchey
0004  * @date June 2018
0005  * @brief Implementation of TpcPrototypeCluster
0006  */
0007 #include "TpcPrototypeCluster.h"
0008
0009 #include <cmath>
0010 #include <utility>  // for swap
0011
0012 namespace
0013 {
0014 // square convenience function
0015 template <class T>
0016 inline constexpr T square(const T& x)
0017 {
0018   return x * x;
0019 }
0020
0021 // get unique index in cov. matrix array from i and j
0022 inline unsigned int covarIndex(unsigned int i, unsigned int j)
0023 {
0024   if (i > j) std::swap(i, j);
0025   return i + 1 + (j + 1) * (j) / 2 - 1;
0026 }
0027
0028 // rotate size or covariant matrix to polar coordinates and return the phi component
0029 template <float (TpcPrototypeCluster::*accessor)(unsigned int, unsigned int) const>
0030 float rotate(const TpcPrototypeCluster* cluster)
0031 {
0032   const auto phi = -std::atan2(cluster->getY(), cluster->getX());
0033   const auto cosphi = std::cos(phi);
0034   const auto sinphi = std::sin(phi);
0035
0036   return square(sinphi) * (cluster->*accessor)(0, 0) +
0037          square(cosphi) * (cluster->*accessor)(1, 1) +
0038          2. * cosphi * sinphi * (cluster->*accessor)(0, 1);
0039 }
0040
0041 }  // namespace
0042
0043 TpcPrototypeCluster::TpcPrototypeCluster()
0044   : m_cluskey(TrkrDefs::CLUSKEYMAX)
0045   , m_isGlobal(true)
0046   , m_adc(0xFFFFFFFF)
0047   , min_sample(-1)
0048   , max_sample(-1)
0049   , min_pad_azimuth(-1)
0050   , max_pad_azimuth(-1)
0051   , peak(NAN)
0052   , peak_sample(NAN)
0053   , pedstal(NAN)
0054   , avg_pad_radial(-1)
0055   , avg_pad_azimuth(NAN)
0056   , size_pad_radial(-1)
0057   , size_pad_azimuth(-1)
0058   , delta_azimuth_bin(NAN)
0059   , delta_z(NAN)
0060 {
0061   for (int i = 0; i < 3; ++i) m_pos[i] = NAN;
0062
0063   for (int j = 0; j < 3; ++j)
0064   {
0065     for (int i = j; i < 3; ++i)
0066     {
0067       TpcPrototypeCluster::setSize(i, j, NAN);
0068       TpcPrototypeCluster::setError(i, j, NAN);
0069     }
0070   }
0071 }
0072
0073 void TpcPrototypeCluster::identify(std::ostream& os) const
0074 {
0075   os << "---TpcPrototypeCluster--------------------" << std::endl;
0076   os << "clusid: " << getClusKey() << std::dec << std::endl;
0077
0078   os << " (x,y,z) =  (" << getPosition(0);
0079   os << ", " << getPosition(1) << ", ";
0080   os << getPosition(2) << ") cm";
0081   if (m_isGlobal)
0082     os << " - global coordinates" << std::endl;
0083   else
0084     os << " - local coordinates" << std::endl;
0085
0086   os << " adc = " << getAdc() << std::endl;
0087
0088   os << " size phi = " << getPhiSize();
0089   os << " cm, size z = " << getZSize() << " cm" << std::endl;
0090
0091   os << "         ( ";
0092   os << getSize(0, 0) << " , ";
0093   os << getSize(0, 1) << " , ";
0094   os << getSize(0, 2) << " )" << std::endl;
0095   os << "  size = ( ";
0096   os << getSize(1, 0) << " , ";
0097   os << getSize(1, 1) << " , ";
0098   os << getSize(1, 2) << " )" << std::endl;
0099   os << "         ( ";
0100   os << getSize(2, 0) << " , ";
0101   os << getSize(2, 1) << " , ";
0102   os << getSize(2, 2) << " )" << std::endl;
0103
0104   os << "         ( ";
0105   os << getError(0, 0) << " , ";
0106   os << getError(0, 1) << " , ";
0107   os << getError(0, 2) << " )" << std::endl;
0108   os << "  err  = ( ";
0109   os << getError(1, 0) << " , ";
0110   os << getError(1, 1) << " , ";
0111   os << getError(1, 2) << " )" << std::endl;
0112   os << "         ( ";
0113   os << getError(2, 0) << " , ";
0114   os << getError(2, 1) << " , ";
0115   os << getError(2, 2) << " )" << std::endl;
0116
0117   //  std::set<int> pad_radials;
0118   os << "  pad_radials  = ( ";
0119   for (auto& element : pad_radials) os << element << " ";
0120   os << ")" << std::endl;
0121
0122   //  std::set<int> pad_azimuths;
0123   os << "  pad_azimuths  = ( ";
0124   for (auto& element : pad_azimuths) os << element << " ";
0125   os << ")" << std::endl;
0126   //  std::set<int> samples;
0127   os << "  samples  = ( ";
0128   for (auto& element : samples) os << element << " ";
0129   os << ")" << std::endl;
0130   //
0131   //  std::map<int, std::vector<double>> pad_radial_samples;
0132   os << "  pad_radial_samples  = [ " << std::endl;
0133   for (auto& pair : pad_radial_samples)
0134   {
0135     os << pair.first << "->( ";
0136     for (auto& element : pair.second) os << element << " ";
0137     os << ")" << std::endl;
0138   }
0139   os << "]" << std::endl;
0140   //  std::map<int, std::vector<double>> pad_azimuth_samples;
0141   os << "  pad_azimuth_samples  = [ " << std::endl;
0142   for (auto& pair : pad_azimuth_samples)
0143   {
0144     os << pair.first << "->( ";
0145     for (auto& element : pair.second) os << element << " ";
0146     os << ")" << std::endl;
0147   }
0148   os << "]" << std::endl;
0149   //  std::vector<double> sum_samples;
0150   os << "  sum_samples  = ( ";
0151   for (auto& element : sum_samples) os << element << " ";
0152   os << ")" << std::endl;
0153   //
0154   //  int min_sample;
0155   os << "  min_sample  =  " << min_sample << std::endl;
0156   //  int max_sample;
0157   os << "  max_sample  =  " << max_sample << std::endl;
0158   //  int min_pad_azimuth;
0159   os << "  min_pad_azimuth  =  " << min_pad_azimuth << std::endl;
0160   //  int max_pad_azimuth;
0161   os << "  max_pad_azimuth  =  " << max_pad_azimuth << std::endl;
0162   //
0163   //  double peak;
0164   os << "  peak  =  " << peak << std::endl;
0165   //  double peak_sample;
0166   os << "  peak_sample  =  " << peak_sample << std::endl;
0167   //  double pedstal;
0168   os << "  pedstal  =  " << pedstal << std::endl;
0169   //
0170   //  //    std::map<int, double> pad_radial_peaks; // radial always have size = 1 in this analysis
0171   //  std::map<int, double> pad_azimuth_peaks;
0172   os << "  pad_azimuth_peaks  = ( ";
0173   for (auto& element : pad_azimuth_peaks) os << element.first << "->" << element.second << " ";
0174   os << ")" << std::endl;
0175   //
0176   //  //! pad coordinate
0177   //  int avg_pad_radial;
0178   os << "  avg_pad_radial  =  " << avg_pad_radial << std::endl;
0179   //  double avg_pad_azimuth;
0180   os << "  avg_pad_azimuth  =  " << avg_pad_azimuth << std::endl;
0181   //
0182   //  //! cluster size in units of pad bins
0183   //  int size_pad_radial;
0184   os << "  size_pad_radial  =  " << size_pad_radial << std::endl;
0185   //  int size_pad_azimuth;
0186   os << "  size_pad_azimuth  =  " << size_pad_azimuth << std::endl;
0187   //
0188   //  //! pad bin size
0189   //  //! phi size per pad in rad
0190   //  double delta_azimuth_bin;
0191   os << "  delta_azimuth_bin  =  " << delta_azimuth_bin << std::endl;
0192   //  //! z size per ADC sample bin
0193   //  double delta_z;
0194   os << "  delta_z  =  " << delta_z << std::endl;
0195
0196   os << std::endl;
0197   os << "-----------------------------------------------" << std::endl;
0198
0199   return;
0200 }
0201
0202 int TpcPrototypeCluster::isValid() const
0203 {
0204   if (m_cluskey == TrkrDefs::CLUSKEYMAX) return 0;
0205   for (int i = 0; i < 3; ++i)
0206   {
0207     if (isnan(getPosition(i))) return 0;
0208   }
0209   if (m_adc == 0xFFFFFFFF) return 0;
0210   for (int j = 0; j < 3; ++j)
0211   {
0212     for (int i = j; i < 3; ++i)
0213     {
0214       if (isnan(getSize(i, j))) return 0;
0215       if (isnan(getError(i, j))) return 0;
0216     }
0217   }
0218
0219   return 1;
0220 }
0221
0222 void TpcPrototypeCluster::setSize(unsigned int i, unsigned int j, float value)
0223 {
0224   m_size[covarIndex(i, j)] = value;
0225   return;
0226 }
0227
0228 float TpcPrototypeCluster::getSize(unsigned int i, unsigned int j) const
0229 {
0230   return m_size[covarIndex(i, j)];
0231 }
0232
0233 void TpcPrototypeCluster::setError(unsigned int i, unsigned int j, float value)
0234 {
0235   m_err[covarIndex(i, j)] = value;
0236   return;
0237 }
0238
0239 float TpcPrototypeCluster::getError(unsigned int i, unsigned int j) const
0240 {
0241   return m_err[covarIndex(i, j)];
0242 }
0243
0244 float TpcPrototypeCluster::getPhiSize() const
0245 {
0246   return 2 * std::sqrt(rotate<&TpcPrototypeCluster::getSize>(this));
0247 }
0248
0249 float TpcPrototypeCluster::getZSize() const
0250 {
0251   return 2. * sqrt(getSize(2, 2));
0252 }
0253
0254 float TpcPrototypeCluster::getPhiError() const
0255 {
0256   const float rad = std::sqrt(square(m_pos[0]) + square(m_pos[1]));
0257   if (rad > 0) return getRPhiError() / rad;
0258   return 0;
0259 }
0260
0261 float TpcPrototypeCluster::getRPhiError() const
0262 {
0263   return std::sqrt(rotate<&TpcPrototypeCluster::getError>(this));
0264 }
0265
0266 float TpcPrototypeCluster::getZError() const
0267 {
0268   return std::sqrt(getError(2, 2));
0269 }