packages/trackbase/TrkrClusterv1.cc

0001 /**
0002  * @file trackbase/TrkrClusterv1.cc
0003  * @author D. McGlinchey
0004  * @date June 2018
0005  * @brief Implementation of TrkrClusterv1
0006  */
0007 #include "TrkrClusterv1.h"
0008
0009 #include <cmath>
0010 #include <utility>  // for swap
0011
0012 namespace
0013 {
0014
0015   // square convenience function
0016   template <class T>
0017   inline constexpr T square(const T& x)
0018   {
0019     return x * x;
0020   }
0021
0022   // get unique index in cov. matrix array from i and j
0023   inline unsigned int covarIndex(unsigned int i, unsigned int j)
0024   {
0025     if (i > j)
0026     {
0027       std::swap(i, j);
0028     }
0029     return i + 1 + (j + 1) * (j) / 2 - 1;
0030   }
0031
0032   // rotate size or covariant matrix to polar coordinates and return the phi component
0033   template <float (TrkrClusterv1::*accessor)(unsigned int, unsigned int) const>
0034   float rotate(const TrkrClusterv1* cluster)
0035   {
0036     const auto phi = -std::atan2(cluster->getY(), cluster->getX());
0037     const auto cosphi = std::cos(phi);
0038     const auto sinphi = std::sin(phi);
0039
0040     return square(sinphi) * (cluster->*accessor)(0, 0) +
0041            square(cosphi) * (cluster->*accessor)(1, 1) +
0042            2. * cosphi * sinphi * (cluster->*accessor)(0, 1);
0043   }
0044
0045 }  // namespace
0046
0047 TrkrClusterv1::TrkrClusterv1()
0048   : m_cluskey(TrkrDefs::CLUSKEYMAX)
0049   , m_isGlobal(true)
0050   , m_adc(0xFFFFFFFF)
0051 {
0052   for (float& m_po : m_pos)
0053   {
0054     m_po = NAN;
0055   }
0056
0057   for (int j = 0; j < 3; ++j)
0058   {
0059     for (int i = 0; i < 3; ++i)
0060     {
0061       setSize(i, j, NAN);
0062       setError(i, j, NAN);
0063     }
0064   }
0065 }
0066
0067 void TrkrClusterv1::identify(std::ostream& os) const
0068 {
0069   os << "---TrkrClusterv1--------------------" << std::endl;
0070
0071   os << " (x,y,z) =  (" << getPosition(0);
0072   os << ", " << getPosition(1) << ", ";
0073   os << getPosition(2) << ") cm";
0074   if (m_isGlobal)
0075   {
0076     os << " - global coordinates" << std::endl;
0077   }
0078   else
0079   {
0080     os << " - local coordinates" << std::endl;
0081   }
0082
0083   os << " adc = " << getAdc() << std::endl;
0084
0085   os << " size phi = " << getPhiSize();
0086   os << " cm, size z = " << getZSize() << " cm" << std::endl;
0087
0088   os << "         ( ";
0089   os << getSize(0, 0) << " , ";
0090   os << getSize(0, 1) << " , ";
0091   os << getSize(0, 2) << " )" << std::endl;
0092   os << "  size = ( ";
0093   os << getSize(1, 0) << " , ";
0094   os << getSize(1, 1) << " , ";
0095   os << getSize(1, 2) << " )" << std::endl;
0096   os << "         ( ";
0097   os << getSize(2, 0) << " , ";
0098   os << getSize(2, 1) << " , ";
0099   os << getSize(2, 2) << " )" << std::endl;
0100
0101   os << "         ( ";
0102   os << getError(0, 0) << " , ";
0103   os << getError(0, 1) << " , ";
0104   os << getError(0, 2) << " )" << std::endl;
0105   os << "  err  = ( ";
0106   os << getError(1, 0) << " , ";
0107   os << getError(1, 1) << " , ";
0108   os << getError(1, 2) << " )" << std::endl;
0109   os << "         ( ";
0110   os << getError(2, 0) << " , ";
0111   os << getError(2, 1) << " , ";
0112   os << getError(2, 2) << " )" << std::endl;
0113
0114   os << std::endl;
0115   os << "-----------------------------------------------" << std::endl;
0116
0117   return;
0118 }
0119
0120 int TrkrClusterv1::isValid() const
0121 {
0122   if (m_cluskey == TrkrDefs::CLUSKEYMAX)
0123   {
0124     return 0;
0125   }
0126   for (int i = 0; i < 3; ++i)
0127   {
0128     if (std::isnan(getPosition(i)))
0129     {
0130       return 0;
0131     }
0132   }
0133   if (m_adc == 0xFFFFFFFF)
0134   {
0135     return 0;
0136   }
0137   for (int j = 0; j < 3; ++j)
0138   {
0139     for (int i = j; i < 3; ++i)
0140     {
0141       if (std::isnan(getSize(i, j)))
0142       {
0143         return 0;
0144       }
0145       if (std::isnan(getError(i, j)))
0146       {
0147         return 0;
0148       }
0149     }
0150   }
0151
0152   return 1;
0153 }
0154
0155 void TrkrClusterv1::CopyFrom(const TrkrCluster& source)
0156 {
0157   // do nothing if copying onto oneself
0158   if (this == &source)
0159   {
0160     return;
0161   }
0162
0163   // parent class method
0164   TrkrCluster::CopyFrom(source);
0165
0166   setX(source.getX());
0167   setY(source.getY());
0168   setZ(source.getZ());
0169   m_isGlobal = source.isGlobal();
0170   setAdc(source.getAdc());
0171
0172   for (int j = 0; j < 3; ++j)
0173   {
0174     for (int i = 0; i < 3; ++i)
0175     {
0176       setSize(i, j, source.getSize(i, j));
0177       setError(i, j, source.getError(i, j));
0178     }
0179   }
0180 }
0181
0182 void TrkrClusterv1::setSize(unsigned int i, unsigned int j, float value)
0183 {
0184   m_size[covarIndex(i, j)] = value;
0185   return;
0186 }
0187
0188 float TrkrClusterv1::getSize(unsigned int i, unsigned int j) const
0189 {
0190   return m_size[covarIndex(i, j)];
0191 }
0192
0193 void TrkrClusterv1::setError(unsigned int i, unsigned int j, float value)
0194 {
0195   m_err[covarIndex(i, j)] = value;
0196   return;
0197 }
0198
0199 float TrkrClusterv1::getError(unsigned int i, unsigned int j) const
0200 {
0201   return m_err[covarIndex(i, j)];
0202 }
0203
0204 float TrkrClusterv1::getPhiSize() const
0205 {
0206   return 2 * std::sqrt(rotate<&TrkrClusterv1::getSize>(this));
0207 }
0208
0209 float TrkrClusterv1::getZSize() const
0210 {
0211   return 2. * std::sqrt(getSize(2, 2));
0212 }
0213
0214 float TrkrClusterv1::getPhiError() const
0215 {
0216   const float rad = std::sqrt(square(m_pos[0]) + square(m_pos[1]));
0217   if (rad > 0)
0218   {
0219     return getRPhiError() / rad;
0220   }
0221   return 0;
0222 }
0223
0224 float TrkrClusterv1::getRPhiError() const
0225 {
0226   return std::sqrt(rotate<&TrkrClusterv1::getError>(this));
0227 }
0228
0229 float TrkrClusterv1::getZError() const
0230 {
0231   return std::sqrt(getError(2, 2));
0232 }