packages/trackbase/TrkrClusterv2.cc

0001 /**
0002  * @file trackbase/TrkrClusterv2.cc
0003  * @author J. Osborn
0004  * @date March 2021
0005  * @brief Implementation of TrkrClusterv2
0006  */
0007 #include "TrkrClusterv2.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 (TrkrClusterv2::*accessor)(unsigned int, unsigned int) const>
0034   float rotate(const TrkrClusterv2* 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 TrkrClusterv2::TrkrClusterv2()
0048   : m_cluskey(TrkrDefs::CLUSKEYMAX)
0049   , m_subsurfkey(TrkrDefs::SUBSURFKEYMAX)
0050   , m_isGlobal(true)
0051   , m_adc(0xFFFFFFFF)
0052 {
0053   for (float& m_po : m_pos)
0054   {
0055     m_po = NAN;
0056   }
0057
0058   for (int j = 0; j < 6; ++j)
0059   {
0060     m_err[j] = NAN;
0061     m_size[j] = NAN;
0062   }
0063   for (int i = 0; i < 2; i++)
0064   {
0065     m_local[i] = NAN;
0066     for (int j = 0; j < 2; j++)
0067     {
0068       m_actsLocalErr[i][j] = NAN;
0069     }
0070   }
0071 }
0072
0073 void TrkrClusterv2::identify(std::ostream& os) const
0074 {
0075   os << "---TrkrClusterv2--------------------" << std::endl;
0076
0077   os << " (x,y,z) =  (" << getPosition(0);
0078   os << ", " << getPosition(1) << ", ";
0079   os << getPosition(2) << ") cm";
0080   if (m_isGlobal)
0081   {
0082     os << " - global coordinates" << std::endl;
0083   }
0084   else
0085   {
0086     os << " - local coordinates" << std::endl;
0087   }
0088
0089   os << " adc = " << getAdc() << std::endl;
0090
0091   os << " size phi = " << getPhiSize();
0092   os << " cm, size z = " << getZSize() << " cm" << std::endl;
0093
0094   os << "         ( ";
0095   os << getSize(0, 0) << " , ";
0096   os << getSize(0, 1) << " , ";
0097   os << getSize(0, 2) << " )" << std::endl;
0098   os << "  size = ( ";
0099   os << getSize(1, 0) << " , ";
0100   os << getSize(1, 1) << " , ";
0101   os << getSize(1, 2) << " )" << std::endl;
0102   os << "         ( ";
0103   os << getSize(2, 0) << " , ";
0104   os << getSize(2, 1) << " , ";
0105   os << getSize(2, 2) << " )" << std::endl;
0106
0107   os << "         ( ";
0108   os << getError(0, 0) << " , ";
0109   os << getError(0, 1) << " , ";
0110   os << getError(0, 2) << " )" << std::endl;
0111   os << "  err  = ( ";
0112   os << getError(1, 0) << " , ";
0113   os << getError(1, 1) << " , ";
0114   os << getError(1, 2) << " )" << std::endl;
0115   os << "         ( ";
0116   os << getError(2, 0) << " , ";
0117   os << getError(2, 1) << " , ";
0118   os << getError(2, 2) << " )" << std::endl;
0119
0120   os << std::endl;
0121   os << "-----------------------------------------------" << std::endl;
0122
0123   return;
0124 }
0125
0126 int TrkrClusterv2::isValid() const
0127 {
0128   if (m_cluskey == TrkrDefs::CLUSKEYMAX)
0129   {
0130     return 0;
0131   }
0132   for (int i = 0; i < 3; ++i)
0133   {
0134     if (std::isnan(getPosition(i)))
0135     {
0136       return 0;
0137     }
0138   }
0139   if (m_adc == 0xFFFFFFFF)
0140   {
0141     return 0;
0142   }
0143   for (int j = 0; j < 3; ++j)
0144   {
0145     for (int i = j; i < 3; ++i)
0146     {
0147       if (std::isnan(getSize(i, j)))
0148       {
0149         return 0;
0150       }
0151       if (std::isnan(getError(i, j)))
0152       {
0153         return 0;
0154       }
0155     }
0156   }
0157
0158   return 1;
0159 }
0160
0161 void TrkrClusterv2::CopyFrom(const TrkrCluster& source)
0162 {
0163   // do nothing if copying onto oneself
0164   if (this == &source)
0165   {
0166     return;
0167   }
0168
0169   // parent class method
0170   TrkrCluster::CopyFrom(source);
0171
0172   setX(source.getX());
0173   setY(source.getY());
0174   setZ(source.getZ());
0175   m_isGlobal = source.isGlobal();
0176   setAdc(source.getAdc());
0177
0178   for (int j = 0; j < 3; ++j)
0179   {
0180     for (int i = 0; i < 3; ++i)
0181     {
0182       setSize(i, j, source.getSize(i, j));
0183       setError(i, j, source.getError(i, j));
0184     }
0185   }
0186
0187   setSubSurfKey(source.getSubSurfKey());
0188   setLocalX(source.getLocalX());
0189   setLocalY(source.getLocalY());
0190
0191   for (int j = 0; j < 2; ++j)
0192   {
0193     for (int i = 0; i < 2; ++i)
0194     {
0195       setActsLocalError(i, j, source.getActsLocalError(i, j));
0196     }
0197   }
0198 }
0199
0200 void TrkrClusterv2::setSize(unsigned int i, unsigned int j, float value)
0201 {
0202   m_size[covarIndex(i, j)] = value;
0203   return;
0204 }
0205
0206 float TrkrClusterv2::getSize(unsigned int i, unsigned int j) const
0207 {
0208   return m_size[covarIndex(i, j)];
0209 }
0210
0211 void TrkrClusterv2::setError(unsigned int i, unsigned int j, float value)
0212 {
0213   m_err[covarIndex(i, j)] = value;
0214   return;
0215 }
0216
0217 float TrkrClusterv2::getError(unsigned int i, unsigned int j) const
0218 {
0219   return m_err[covarIndex(i, j)];
0220 }
0221
0222 float TrkrClusterv2::getPhiSize() const
0223 {
0224   return 2 * std::sqrt(rotate<&TrkrClusterv2::getSize>(this));
0225 }
0226
0227 float TrkrClusterv2::getZSize() const
0228 {
0229   return 2. * std::sqrt(getSize(2, 2));
0230 }
0231
0232 float TrkrClusterv2::getPhiError() const
0233 {
0234   const float rad = std::sqrt(square(m_pos[0]) + square(m_pos[1]));
0235   if (rad > 0)
0236   {
0237     return getRPhiError() / rad;
0238   }
0239   return 0;
0240 }
0241
0242 float TrkrClusterv2::getRPhiError() const
0243 {
0244   return std::sqrt(rotate<&TrkrClusterv2::getError>(this));
0245 }
0246
0247 float TrkrClusterv2::getZError() const
0248 {
0249   return std::sqrt(getError(2, 2));
0250 }
0251
0252 void TrkrClusterv2::setActsLocalError(unsigned int i, unsigned int j,
0253                                       float value)
0254 {
0255   m_actsLocalErr[i][j] = value;
0256 }