packages/tpccalib/TpcSpaceChargeMatrixContainerv1.cc

0001
0002 /**
0003  * @file tpccalib/TpcSpaceChargeMatrixContainerv1.cc
0004  * @author Hugo Pereira Da Costa
0005  * @date June 2018
0006  * @brief Contains matrices needed for space charge trackbase reconstruction
0007  */
0008
0009 #include "TpcSpaceChargeMatrixContainerv1.h"
0010
0011 //___________________________________________________________
0012 TpcSpaceChargeMatrixContainerv1::TpcSpaceChargeMatrixContainerv1()
0013 {
0014   // reset all matrix arrays
0015   TpcSpaceChargeMatrixContainerv1::Reset();
0016 }
0017
0018 //___________________________________________________________
0019 void TpcSpaceChargeMatrixContainerv1::identify(std::ostream& out) const
0020 {
0021   out << "TpcSpaceChargeMatrixContainerv1" << std::endl;
0022   out << "  phibins: " << m_phibins << std::endl;
0023   out << "  rbins: " << m_rbins << std::endl;
0024   out << "  zbins: " << m_zbins << std::endl;
0025 }
0026
0027 //___________________________________________________________
0028 void TpcSpaceChargeMatrixContainerv1::get_grid_dimensions(int& phibins, int& rbins, int& zbins) const
0029 {
0030   phibins = m_phibins;
0031   rbins = m_rbins;
0032   zbins = m_zbins;
0033 }
0034
0035 //___________________________________________________________
0036 int TpcSpaceChargeMatrixContainerv1::get_grid_size() const
0037 {
0038   return m_phibins * m_rbins * m_zbins;
0039 }
0040
0041 //___________________________________________________________
0042 int TpcSpaceChargeMatrixContainerv1::get_cell_index(int iphi, int ir, int iz) const
0043 {
0044   if (iphi < 0 || iphi >= m_phibins)
0045   {
0046     return -1;
0047   }
0048   if (ir < 0 || ir >= m_rbins)
0049   {
0050     return -1;
0051   }
0052   if (iz < 0 || iz >= m_zbins)
0053   {
0054     return -1;
0055   }
0056   return iz + m_zbins * (ir + m_rbins * iphi);
0057 }
0058
0059 //___________________________________________________________
0060 int TpcSpaceChargeMatrixContainerv1::get_entries(int cell_index) const
0061 {
0062   // bound check
0063   if (!bound_check(cell_index))
0064   {
0065     return 0;
0066   }
0067   return m_entries[cell_index];
0068 }
0069
0070 //___________________________________________________________
0071 float TpcSpaceChargeMatrixContainerv1::get_lhs(int cell_index, int i, int j) const
0072 {
0073   // bound check
0074   if (!bound_check(cell_index, i, j))
0075   {
0076     return 0;
0077   }
0078   return m_lhs[cell_index][get_flat_index(i, j)];
0079 }
0080
0081 //___________________________________________________________
0082 float TpcSpaceChargeMatrixContainerv1::get_rhs(int cell_index, int i) const
0083 {
0084   // bound check
0085   if (!bound_check(cell_index, i))
0086   {
0087     return 0;
0088   }
0089   return m_rhs[cell_index][i];
0090 }
0091
0092 //___________________________________________________________
0093 void TpcSpaceChargeMatrixContainerv1::Reset()
0094 {
0095   // reset total number of bins
0096   const int totalbins = m_phibins * m_rbins * m_zbins;
0097
0098   // reset arrays
0099   m_entries = std::vector<int>(totalbins, 0);
0100   m_lhs = std::vector<matrix_t>(totalbins, {{}});
0101   m_rhs = std::vector<column_t>(totalbins, {{}});
0102 }
0103
0104 //___________________________________________________________
0105 void TpcSpaceChargeMatrixContainerv1::set_grid_dimensions(int phibins, int rbins, int zbins)
0106 {
0107   m_phibins = phibins;
0108   m_rbins = rbins;
0109   m_zbins = zbins;
0110   Reset();
0111 }
0112
0113 //___________________________________________________________
0114 void TpcSpaceChargeMatrixContainerv1::add_to_entries(int cell_index, int value)
0115 {
0116   if (bound_check(cell_index))
0117   {
0118     m_entries[cell_index] += value;
0119   }
0120 }
0121
0122 //___________________________________________________________
0123 void TpcSpaceChargeMatrixContainerv1::add_to_lhs(int cell_index, int i, int j, float value)
0124 {
0125   if (bound_check(cell_index, i, j))
0126   {
0127     m_lhs[cell_index][get_flat_index(i, j)] += value;
0128   }
0129 }
0130
0131 //___________________________________________________________
0132 void TpcSpaceChargeMatrixContainerv1::add_to_rhs(int cell_index, int i, float value)
0133 {
0134   if (bound_check(cell_index, i))
0135   {
0136     m_rhs[cell_index][i] += value;
0137   }
0138 }
0139
0140 //___________________________________________________________
0141 bool TpcSpaceChargeMatrixContainerv1::add(const TpcSpaceChargeMatrixContainer& other)
0142 {
0143   // check dimensions
0144   int phibins = 0;
0145   int rbins = 0;
0146   int zbins = 0;
0147   other.get_grid_dimensions(phibins, rbins, zbins);
0148   if ((m_phibins != phibins) || (m_rbins != rbins) || (m_zbins != zbins))
0149   {
0150     std::cout << "TpcSpaceChargeMatrixContainerv1::add - inconsistent grid sizes" << std::endl;
0151     return false;
0152   }
0153
0154   // increment cell entries
0155   for (size_t cell_index = 0; cell_index < m_lhs.size(); ++cell_index)
0156   {
0157     add_to_entries(cell_index, other.get_entries(cell_index));
0158   }
0159
0160   // increment left hand side matrices
0161   for (size_t cell_index = 0; cell_index < m_lhs.size(); ++cell_index)
0162   {
0163     for (int i = 0; i < m_ncoord; ++i)
0164     {
0165       for (int j = 0; j < m_ncoord; ++j)
0166       {
0167         add_to_lhs(cell_index, i, j, other.get_lhs(cell_index, i, j));
0168       }
0169     }
0170   }
0171
0172   // increment right hand side matrices
0173   for (size_t cell_index = 0; cell_index < m_lhs.size(); ++cell_index)
0174   {
0175     for (int i = 0; i < m_ncoord; ++i)
0176     {
0177       add_to_rhs(cell_index, i, other.get_rhs(cell_index, i));
0178     }
0179   }
0180
0181   return true;
0182 }
0183
0184 //___________________________________________________________
0185 bool TpcSpaceChargeMatrixContainerv1::bound_check(int cell_index) const
0186 {
0187   if (cell_index < 0 || cell_index >= (int) m_rhs.size())
0188   {
0189     return false;
0190   }
0191   return true;
0192 }
0193
0194 //___________________________________________________________
0195 bool TpcSpaceChargeMatrixContainerv1::bound_check(int cell_index, int i) const
0196 {
0197   if (cell_index < 0 || cell_index >= (int) m_rhs.size())
0198   {
0199     return false;
0200   }
0201   if (i < 0 || i >= m_ncoord)
0202   {
0203     return false;
0204   }
0205   return true;
0206 }
0207
0208 //___________________________________________________________
0209 bool TpcSpaceChargeMatrixContainerv1::bound_check(int cell_index, int i, int j) const
0210 {
0211   if (cell_index < 0 || cell_index >= (int) m_lhs.size())
0212   {
0213     return false;
0214   }
0215   if (i < 0 || i >= m_ncoord)
0216   {
0217     return false;
0218   }
0219   if (j < 0 || j >= m_ncoord)
0220   {
0221     return false;
0222   }
0223   return true;
0224 }
0225
0226 //___________________________________________________________
0227 int TpcSpaceChargeMatrixContainerv1::get_flat_index(int i, int j) const
0228 {
0229   return j + i * m_ncoord;
0230 }