MkFitCore/src/HitStructures.cc

0001 #include "RecoTracker/MkFitCore/interface/HitStructures.h"
0002
0003 #include "RecoTracker/MkFitCore/interface/IterationConfig.h"
0004 #include "Matriplex/Memory.h"
0005
0006 #include "Debug.h"
0007
0008 namespace mkfit {
0009
0010   void LayerOfHits::Initializator::setup(float qmin, float qmax, float dq) {
0011     assert(qmax > qmin);
0012     float extent = qmax - qmin;
0013     m_nq = std::ceil(extent / dq);
0014     float extra = 0.5f * (m_nq * dq - extent);
0015     m_qmin = qmin - extra;
0016     m_qmax = qmax + extra;
0017   }
0018
0019   LayerOfHits::Initializator::Initializator(const LayerInfo &li, float qmin, float qmax, unsigned int nq)
0020       : m_linfo(li), m_qmin(qmin), m_qmax(qmax), m_nq(nq) {}
0021
0022   LayerOfHits::Initializator::Initializator(const LayerInfo &li, float qmin, float qmax, float dq) : m_linfo(li) {
0023     setup(qmin, qmax, dq);
0024   }
0025
0026   LayerOfHits::Initializator::Initializator(const LayerInfo &li) : m_linfo(li) {
0027     if (li.is_barrel())
0028       setup(li.zmin(), li.zmax(), li.q_bin());
0029     else
0030       setup(li.rin(), li.rout(), li.q_bin());
0031   }
0032
0033   LayerOfHits::LayerOfHits(const LayerOfHits::Initializator &i)
0034       : m_ax_phi(-Const::PI, Const::PI),
0035         m_ax_eta(i.m_qmin, i.m_qmax, i.m_nq),
0036         m_binnor(m_ax_phi, m_ax_eta, true, false)  // yes-radix, no-keep-cons
0037   {
0038     m_layer_info = &i.m_linfo;
0039     m_is_barrel = m_layer_info->is_barrel();
0040
0041     m_dead_bins.resize(m_ax_eta.size_of_N() * m_ax_phi.size_of_N());
0042   }
0043
0044   LayerOfHits::~LayerOfHits() {
0045 #ifdef COPY_SORTED_HITS
0046     free_hits();
0047 #endif
0048   }
0049
0050 #ifdef COPY_SORTED_HITS
0051   void LayerOfHits::alloc_hits(int size) {
0052     m_hits = (Hit *)Matriplex::aligned_alloc64(sizeof(Hit) * size);
0053     m_capacity = size;
0054     for (int ihit = 0; ihit < m_capacity; ihit++) {
0055       m_hits[ihit] = Hit();
0056     }
0057   }
0058
0059   void LayerOfHits::free_hits() { std::free(m_hits); }
0060 #endif
0061
0062   //==============================================================================
0063
0064   void LayerOfHits::reset() {
0065     m_hit_infos.clear();
0066     m_ext_idcs.clear();
0067     m_min_ext_idx = std::numeric_limits<unsigned int>::max();
0068     m_max_ext_idx = std::numeric_limits<unsigned int>::min();
0069     m_n_hits = 0;
0070     m_binnor.reset_contents();
0071   }
0072
0073   //==============================================================================
0074
0075   void LayerOfHits::suckInHits(const HitVec &hitv) {
0076     m_ext_hits = &hitv;
0077     m_n_hits = hitv.size();
0078
0079     m_binnor.begin_registration(m_n_hits);
0080
0081 #ifdef COPY_SORTED_HITS
0082     if (m_capacity < m_n_hits) {
0083       free_hits();
0084       alloc_hits(m_n_hits);
0085     }
0086 #endif
0087
0088     std::vector<HitInfo> hinfos;
0089     if (Config::usePhiQArrays) {
0090       hinfos.reserve(m_n_hits);
0091       m_hit_infos.reserve(m_n_hits);
0092     }
0093
0094     // Factor to get from hit sigma to half-length in q direction.
0095     const float hl_fac = is_pixel() ? 3.0f : std::sqrt(3.0f);
0096
0097     for (unsigned int i = 0; i < m_n_hits; ++i) {
0098       const Hit &h = hitv[i];
0099
0100       float phi = h.phi();
0101       float q = m_is_barrel ? h.z() : h.r();
0102
0103       m_binnor.register_entry_safe(phi, q);
0104
0105       if (Config::usePhiQArrays) {
0106         float half_length, qbar;
0107         if (m_is_barrel) {
0108           half_length = hl_fac * std::sqrt(h.ezz());
0109           qbar = h.r();
0110         } else {
0111           half_length = hl_fac * std::sqrt(h.exx() + h.eyy());
0112           qbar = h.z();
0113         }
0114         hinfos.emplace_back(HitInfo({phi, q, half_length, qbar}));
0115       }
0116     }
0117
0118     m_binnor.finalize_registration();
0119
0120     for (unsigned int i = 0; i < m_n_hits; ++i) {
0121       unsigned int j = m_binnor.m_ranks[i];
0122 #ifdef COPY_SORTED_HITS
0123       memcpy(&m_hits[i], &hitv[j], sizeof(Hit));
0124 #endif
0125       if (Config::usePhiQArrays) {
0126         m_hit_infos.emplace_back(hinfos[j]);
0127       }
0128     }
0129   }
0130
0131   //==============================================================================
0132
0133   void LayerOfHits::suckInDeads(const DeadVec &deadv) {
0134     m_dead_bins.assign(m_dead_bins.size(), false);
0135
0136     for (const auto &d : deadv) {
0137       bin_index_t q_bin_1 = qBinChecked(d.q1);
0138       bin_index_t q_bin_2 = qBinChecked(d.q2) + 1;
0139       bin_index_t phi_bin_1 = phiBin(d.phi1);
0140       bin_index_t phi_bin_2 = phiMaskApply(phiBin(d.phi2) + 1);
0141
0142       for (bin_index_t q_bin = q_bin_1; q_bin != q_bin_2; q_bin++) {
0143         const unsigned int qoff = q_bin * m_ax_phi.size_of_N();
0144         for (bin_index_t pb = phi_bin_1; pb != phi_bin_2; pb = phiMaskApply(pb + 1)) {
0145           m_dead_bins[qoff + pb] = true;
0146         }
0147       }
0148     }
0149   }
0150
0151   //==============================================================================
0152
0153   void LayerOfHits::beginRegistrationOfHits(const HitVec &hitv) {
0154     m_ext_hits = &hitv;
0155     m_n_hits = 0;
0156
0157     m_binnor.begin_registration(128);  // initial reserve for cons vectors
0158   }
0159
0160   void LayerOfHits::registerHit(unsigned int idx) {
0161     const Hit &h = (*m_ext_hits)[idx];
0162
0163     m_ext_idcs.push_back(idx);
0164     m_min_ext_idx = std::min(m_min_ext_idx, idx);
0165     m_max_ext_idx = std::max(m_max_ext_idx, idx);
0166
0167     float phi = h.phi();
0168     float q = m_is_barrel ? h.z() : h.r();
0169
0170     m_binnor.register_entry_safe(phi, q);
0171
0172     if (Config::usePhiQArrays) {
0173       // Factor to get from hit sigma to half-length in q direction.
0174       const float hl_fac = is_pixel() ? 3.0f : std::sqrt(3.0f);
0175       float half_length, qbar;
0176       if (m_is_barrel) {
0177         half_length = hl_fac * std::sqrt(h.ezz());
0178         qbar = h.r();
0179       } else {
0180         half_length = hl_fac * std::sqrt(h.exx() + h.eyy());
0181         qbar = h.z();
0182       }
0183       m_hit_infos.emplace_back(HitInfo({phi, q, half_length, qbar}));
0184     }
0185   }
0186
0187   void LayerOfHits::endRegistrationOfHits(bool build_original_to_internal_map) {
0188     m_n_hits = m_ext_idcs.size();
0189     if (m_n_hits == 0)
0190       return;
0191
0192     m_binnor.finalize_registration();
0193
0194     // copy q/phi
0195
0196 #ifdef COPY_SORTED_HITS
0197     if (m_capacity < m_n_hits) {
0198       free_hits();
0199       alloc_hits(m_n_hits);
0200     }
0201 #endif
0202
0203     std::vector<HitInfo> hinfos;
0204     if (Config::usePhiQArrays) {
0205       hinfos.swap(m_hit_infos);
0206       m_hit_infos.reserve(m_n_hits);
0207     }
0208
0209     for (unsigned int i = 0; i < m_n_hits; ++i) {
0210       unsigned int j = m_binnor.m_ranks[i];  // index in intermediate
0211       unsigned int k = m_ext_idcs[j];        // index in external hit_vec
0212
0213 #ifdef COPY_SORTED_HITS
0214       memcpy(&m_hits[i], &hitv[k], sizeof(Hit));
0215 #endif
0216
0217       if (Config::usePhiQArrays) {
0218         m_hit_infos.emplace_back(hinfos[j]);
0219       }
0220
0221       // Redirect m_binnor.m_ranks[i] to point to external/original index.
0222       m_binnor.m_ranks[i] = k;
0223     }
0224
0225     if (build_original_to_internal_map) {
0226       if (m_max_ext_idx - m_min_ext_idx + 1 > 8 * m_n_hits) {
0227         // If this happens we might:
0228         // a) Use external indices for everything. -- *** We are now. ***
0229         // b) Build these maps for seeding layers only.
0230         // c) Have a flag in hit-on-track that tells us if the hit index has been remapped,
0231         //    essentially, if it is a seed hit. This might be smart anyway.
0232         //    One could use index < -256 or something similar.
0233
0234         printf(
0235             "LayerOfHits::endRegistrationOfHits() original_to_internal index map vector is largish: m_n_hits=%d, "
0236             "map_vector_size=%d\n",
0237             m_n_hits,
0238             m_max_ext_idx - m_min_ext_idx + 1);
0239       }
0240
0241       m_ext_idcs.resize(m_max_ext_idx - m_min_ext_idx + 1);
0242       for (unsigned int i = 0; i < m_n_hits; ++i) {
0243         m_ext_idcs[m_hit_ranks[i] - m_min_ext_idx] = i;
0244       }
0245     }
0246
0247     // We can release m_hit_infos and, if not used, also m_ext_idcs -- and realloc them
0248     // on next beginRegistration().
0249     // If binnor had keep_cons on we could use it for pre-selection in selectHitIndices()
0250     // instead of q and phi arrays -- assuming sufficient precision can be achieved..
0251   }
0252
0253   void LayerOfHits::printBins() {
0254     for (bin_index_t qb = 0; qb <= m_ax_eta.m_last_N_bin; ++qb) {
0255       printf("%c bin %d\n", is_barrel() ? 'Z' : 'R', qb);
0256       for (bin_index_t pb = 0; pb <= m_ax_phi.m_last_N_bin; ++pb) {
0257         if (pb % 8 == 0)
0258           printf(" Phi %4d: ", pb);
0259         auto content = m_binnor.get_content(pb, qb);
0260         printf("%5d,%4d   %s", content.first, content.count, ((pb + 1) % 8 == 0) ? "\n" : "");
0261       }
0262     }
0263   }
0264
0265   //==============================================================================
0266   // EventOfHits
0267   //==============================================================================
0268
0269   EventOfHits::EventOfHits(const TrackerInfo &trk_inf) : m_n_layers(trk_inf.n_layers()) {
0270     m_layers_of_hits.reserve(trk_inf.n_layers());
0271     for (int ii = 0; ii < trk_inf.n_layers(); ++ii) {
0272       m_layers_of_hits.emplace_back(LayerOfHits::Initializator(trk_inf.layer(ii)));
0273     }
0274   }
0275
0276 }  // end namespace mkfit