packages/tpctrackreco/IdealPadMap.cc

0001 #include "IdealPadMap.h"
0002
0003 #include <cdbobjects/CDBTTree.h>
0004 #include <ffamodules/CDBInterface.h>
0005
0006 #include <algorithm>
0007 #include <cmath>
0008 #include <iostream>
0009 #include <limits>
0010
0011 IdealPadMap::IdealPadMap()
0012 {
0013   m_radius_cm_by_layer.fill(std::numeric_limits<double>::quiet_NaN());
0014   m_layer_by_key.fill(-1);
0015 }
0016
0017 IdealPadMap::~IdealPadMap() = default;
0018
0019 int IdealPadMap::layer_index(const unsigned int layer) const
0020 {
0021   if (layer < FIRST_LAYER || layer > LAST_LAYER)
0022   {
0023     return -1;
0024   }
0025   return static_cast<int>(layer - FIRST_LAYER);
0026 }
0027
0028 int IdealPadMap::get_region(const unsigned int layer) const
0029 {
0030   const int ilayer = layer_index(layer);
0031   if (ilayer < 0)
0032   {
0033     return -1;
0034   }
0035   return ilayer / 16;
0036 }
0037
0038 double IdealPadMap::wrap_phi(const double phi) const
0039 {
0040   double out = phi;
0041   while (out <= -M_PI)
0042   {
0043     out += 2.0 * M_PI;
0044   }
0045   while (out > M_PI)
0046   {
0047     out -= 2.0 * M_PI;
0048   }
0049   return out;
0050 }
0051
0052 int IdealPadMap::load_from_cdb(const int /*verbosity*/)
0053 {
0054   m_is_loaded = false;
0055
0056   for (auto& v : m_cdb_phi_by_layer)
0057   {
0058     v.clear();
0059   }
0060   m_radius_cm_by_layer.fill(std::numeric_limits<double>::quiet_NaN());
0061   m_layer_by_key.fill(-1);
0062
0063   CDBInterface* cdb = CDBInterface::instance();
0064   const std::string calibdir = cdb->getUrl("TPC_FEE_CHANNEL_MAP");
0065
0066   if (calibdir.empty())
0067   {
0068     std::cout << "IdealPadMap::load_from_cdb - no TPC_FEE_CHANNEL_MAP found" << std::endl;
0069     return -1;
0070   }
0071
0072   CDBTTree* cdbttree = new CDBTTree(calibdir);
0073   cdbttree->LoadCalibrations();
0074
0075   std::array<std::vector<double>, N_LAYERS> radius_mm_by_layer;
0076
0077   for (unsigned int fee = 0; fee < N_FEE; ++fee)
0078   {
0079     for (unsigned int ch = 0; ch < N_CH; ++ch)
0080     {
0081       const unsigned int key = 256U * fee + ch;
0082       const int layer = cdbttree->GetIntValue(key, "layer");
0083       m_layer_by_key[key] = layer;
0084
0085       const int ilayer = layer_index(static_cast<unsigned int>(layer));
0086       if (ilayer < 0)
0087       {
0088         continue;
0089       }
0090
0091       m_cdb_phi_by_layer[ilayer].push_back(cdbttree->GetDoubleValue(key, "phi"));
0092       radius_mm_by_layer[ilayer].push_back(cdbttree->GetDoubleValue(key, "R"));
0093     }
0094   }
0095
0096   delete cdbttree;
0097   cdbttree = nullptr;
0098
0099   for (unsigned int ilayer = 0; ilayer < N_LAYERS; ++ilayer)
0100   {
0101     if (m_cdb_phi_by_layer[ilayer].empty() || radius_mm_by_layer[ilayer].empty())
0102     {
0103       std::cout << "IdealPadMap::load_from_cdb - missing CDB entries for layer "
0104                 << ilayer + FIRST_LAYER << std::endl;
0105       return -1;
0106     }
0107
0108     std::sort(m_cdb_phi_by_layer[ilayer].begin(), m_cdb_phi_by_layer[ilayer].end());
0109
0110     double radius_mm = 0.0;
0111     for (const double r : radius_mm_by_layer[ilayer])
0112     {
0113       radius_mm += r;
0114     }
0115     radius_mm /= static_cast<double>(radius_mm_by_layer[ilayer].size());
0116
0117     // CDB R is in mm.  Most TPC tracking/display code uses cm.
0118     m_radius_cm_by_layer[ilayer] = radius_mm / 10.0;
0119   }
0120
0121   m_is_loaded = true;
0122
0123   // if (verbosity > 0)
0124   {
0125     std::cout << "IdealPadMap::load_from_cdb - loaded " << calibdir << std::endl;
0126     for (unsigned int region = 0; region < N_REGIONS; ++region)
0127     {
0128       const unsigned int first_layer = FIRST_LAYER + 16U * region;
0129       const unsigned int pads_per_sector = get_pads_per_sector(region);
0130       std::cout << "  region " << region
0131                 << " first_layer " << first_layer
0132                 << " pads_per_sector " << pads_per_sector
0133                 << " total_phibins " << get_total_phibins(first_layer)
0134                 << std::endl;
0135       if (pads_per_sector == 0U)
0136       {
0137         continue;
0138       }
0139
0140       const unsigned int first_pad = 0U;
0141       const unsigned int last_pad = pads_per_sector - 1U;
0142       for (unsigned int side = 0; side < N_SIDES; ++side)
0143       {
0144         for (unsigned int sector = 0; sector < N_SECTORS; ++sector)
0145         {
0146           std::cout << "    side " << side
0147                     << " sector " << sector
0148                     << " first_pad " << first_pad
0149                     << " first_phi " << get_phi(side, sector, first_layer, first_pad)
0150                     << " last_pad " << last_pad
0151                     << " last_phi " << get_phi(side, sector, first_layer, last_pad)
0152                     << std::endl;
0153         }
0154       }
0155     }
0156   }
0157
0158   return 0;
0159 }
0160
0161 unsigned int IdealPadMap::get_pads_per_sector(const unsigned int region) const
0162 {
0163   if (region >= N_REGIONS)
0164   {
0165     return 0U;
0166   }
0167
0168   const unsigned int layer = FIRST_LAYER + 16U * region;
0169   return get_pads_per_sector_for_layer(layer);
0170 }
0171
0172 unsigned int IdealPadMap::get_pads_per_sector_for_layer(const unsigned int layer) const
0173 {
0174   const int ilayer = layer_index(layer);
0175   if (ilayer < 0)
0176   {
0177     return 0U;
0178   }
0179   return static_cast<unsigned int>(m_cdb_phi_by_layer[ilayer].size());
0180 }
0181
0182 unsigned int IdealPadMap::get_total_phibins(const unsigned int layer) const
0183 {
0184   return N_SECTORS * get_pads_per_sector_for_layer(layer);
0185 }
0186
0187 double IdealPadMap::get_radius(const unsigned int layer) const
0188 {
0189   const int ilayer = layer_index(layer);
0190   if (ilayer < 0)
0191   {
0192     return std::numeric_limits<double>::quiet_NaN();
0193   }
0194   return m_radius_cm_by_layer[ilayer];
0195 }
0196
0197 double IdealPadMap::get_layer_thickness(const unsigned int layer) const
0198 {
0199   const int ilayer = layer_index(layer);
0200   if (ilayer < 0)
0201   {
0202     return std::numeric_limits<double>::quiet_NaN();
0203   }
0204
0205   if (ilayer == 0)
0206   {
0207     return m_radius_cm_by_layer[1] - m_radius_cm_by_layer[0];
0208   }
0209   if (ilayer == static_cast<int>(N_LAYERS - 1))
0210   {
0211     return m_radius_cm_by_layer[N_LAYERS - 1] - m_radius_cm_by_layer[N_LAYERS - 2];
0212   }
0213
0214   return 0.5 * (m_radius_cm_by_layer[ilayer + 1] - m_radius_cm_by_layer[ilayer - 1]);
0215 }
0216
0217 double IdealPadMap::get_cdb_local_phi(const unsigned int layer,
0218                                       const unsigned int local_phibin) const
0219 {
0220   const int ilayer = layer_index(layer);
0221   if (ilayer < 0)
0222   {
0223     return std::numeric_limits<double>::quiet_NaN();
0224   }
0225
0226   const auto& phi_vec = m_cdb_phi_by_layer[ilayer];
0227   if (local_phibin >= phi_vec.size())
0228   {
0229     return std::numeric_limits<double>::quiet_NaN();
0230   }
0231
0232   return phi_vec[local_phibin];
0233 }
0234
0235 double IdealPadMap::get_phi(const unsigned int side,
0236                             const unsigned int layer,
0237                             const unsigned int phibin) const
0238 {
0239   const unsigned int pads_per_sector = get_pads_per_sector_for_layer(layer);
0240   if (pads_per_sector == 0U)
0241   {
0242     return std::numeric_limits<double>::quiet_NaN();
0243   }
0244
0245   const unsigned int sector = phibin / pads_per_sector;
0246   const unsigned int local_phibin = phibin % pads_per_sector;
0247
0248   return get_phi(side, sector, layer, local_phibin);
0249 }
0250
0251 double IdealPadMap::get_phi(const unsigned int side,
0252                             const unsigned int sector,
0253                             const unsigned int layer,
0254                             const unsigned int local_phibin) const
0255 {
0256   if (side >= N_SIDES)
0257   {
0258     return std::numeric_limits<double>::quiet_NaN();
0259   }
0260   if (sector >= N_SECTORS)
0261   {
0262     return std::numeric_limits<double>::quiet_NaN();
0263   }
0264
0265   unsigned int lookup_phibin = local_phibin;
0266   if (side == 1U)
0267   {
0268     const unsigned int pads_per_sector = get_pads_per_sector_for_layer(layer);
0269     if (pads_per_sector == 0U)
0270     {
0271       return std::numeric_limits<double>::quiet_NaN();
0272     }
0273     if (local_phibin >= pads_per_sector)
0274     {
0275       return std::numeric_limits<double>::quiet_NaN();
0276     }
0277     lookup_phibin = pads_per_sector - 1U - local_phibin;
0278   }
0279
0280   const double cdb_phi = get_cdb_local_phi(layer, lookup_phibin);
0281   if (!std::isfinite(cdb_phi))
0282   {
0283     return std::numeric_limits<double>::quiet_NaN();
0284   }
0285
0286   const unsigned int mapped_sector = (5U + N_SECTORS - (sector % N_SECTORS)) % N_SECTORS;
0287
0288   const double phi = ((side == 1U ? 1.0 : -1.0) * (cdb_phi - M_PI / 2.0)) + (static_cast<double>(mapped_sector) * M_PI / 6.0);
0289
0290   return wrap_phi(phi);
0291 }
0292 int IdealPadMap::get_layer_from_fee_channel(const unsigned int fee,
0293                                             const unsigned int channel) const
0294 {
0295   if (fee >= N_FEE || channel >= N_CH)
0296   {
0297     return -1;
0298   }
0299   const unsigned int key = 256U * fee + channel;
0300   return m_layer_by_key[key];
0301 }