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     for (unsigned int i = 0; i < m_n_hits; ++i) {
0095       const Hit &h = hitv[i];
0096
0097       float phi = h.phi();
0098       float q = m_is_barrel ? h.z() : h.r();
0099
0100       m_binnor.register_entry_safe(phi, q);
0101
0102       if (Config::usePhiQArrays) {
0103         const float sqrt3 = std::sqrt(3);
0104         float half_length, qbar;
0105         if (m_is_barrel) {
0106           half_length = sqrt3 * std::sqrt(h.ezz());
0107           qbar = h.r();
0108         } else {
0109           half_length = sqrt3 * std::sqrt(h.exx() + h.eyy());
0110           qbar = h.z();
0111         }
0112         hinfos.emplace_back(HitInfo({phi, q, half_length, qbar}));
0113       }
0114     }
0115
0116     m_binnor.finalize_registration();
0117
0118     for (unsigned int i = 0; i < m_n_hits; ++i) {
0119       unsigned int j = m_binnor.m_ranks[i];
0120 #ifdef COPY_SORTED_HITS
0121       memcpy(&m_hits[i], &hitv[j], sizeof(Hit));
0122 #endif
0123       if (Config::usePhiQArrays) {
0124         m_hit_infos.emplace_back(hinfos[j]);
0125       }
0126     }
0127   }
0128
0129   //==============================================================================
0130
0131   void LayerOfHits::suckInDeads(const DeadVec &deadv) {
0132     m_dead_bins.assign(m_dead_bins.size(), false);
0133
0134     for (const auto &d : deadv) {
0135       bin_index_t q_bin_1 = qBinChecked(d.q1);
0136       bin_index_t q_bin_2 = qBinChecked(d.q2) + 1;
0137       bin_index_t phi_bin_1 = phiBin(d.phi1);
0138       bin_index_t phi_bin_2 = phiMaskApply(phiBin(d.phi2) + 1);
0139
0140       for (bin_index_t q_bin = q_bin_1; q_bin != q_bin_2; q_bin++) {
0141         const unsigned int qoff = q_bin * m_ax_phi.size_of_N();
0142         for (bin_index_t pb = phi_bin_1; pb != phi_bin_2; pb = phiMaskApply(pb + 1)) {
0143           m_dead_bins[qoff + pb] = true;
0144         }
0145       }
0146     }
0147   }
0148
0149   //==============================================================================
0150
0151   void LayerOfHits::beginRegistrationOfHits(const HitVec &hitv) {
0152     m_ext_hits = &hitv;
0153     m_n_hits = 0;
0154
0155     m_binnor.begin_registration(128);  // initial reserve for cons vectors
0156   }
0157
0158   void LayerOfHits::registerHit(unsigned int idx) {
0159     const Hit &h = (*m_ext_hits)[idx];
0160
0161     m_ext_idcs.push_back(idx);
0162     m_min_ext_idx = std::min(m_min_ext_idx, idx);
0163     m_max_ext_idx = std::max(m_max_ext_idx, idx);
0164
0165     float phi = h.phi();
0166     float q = m_is_barrel ? h.z() : h.r();
0167
0168     m_binnor.register_entry_safe(phi, q);
0169
0170     if (Config::usePhiQArrays) {
0171       const float sqrt3 = std::sqrt(3);
0172       float half_length, qbar;
0173       if (m_is_barrel) {
0174         half_length = sqrt3 * std::sqrt(h.ezz());
0175         qbar = h.r();
0176       } else {
0177         half_length = sqrt3 * std::sqrt(h.exx() + h.eyy());
0178         qbar = h.z();
0179       }
0180       m_hit_infos.emplace_back(HitInfo({phi, q, half_length, qbar}));
0181     }
0182   }
0183
0184   void LayerOfHits::endRegistrationOfHits(bool build_original_to_internal_map) {
0185     m_n_hits = m_ext_idcs.size();
0186     if (m_n_hits == 0)
0187       return;
0188
0189     m_binnor.finalize_registration();
0190
0191     // copy q/phi
0192
0193 #ifdef COPY_SORTED_HITS
0194     if (m_capacity < m_n_hits) {
0195       free_hits();
0196       alloc_hits(m_n_hits);
0197     }
0198 #endif
0199
0200     std::vector<HitInfo> hinfos;
0201     if (Config::usePhiQArrays) {
0202       hinfos.swap(m_hit_infos);
0203       m_hit_infos.reserve(m_n_hits);
0204     }
0205
0206     for (unsigned int i = 0; i < m_n_hits; ++i) {
0207       unsigned int j = m_binnor.m_ranks[i];  // index in intermediate
0208       unsigned int k = m_ext_idcs[j];        // index in external hit_vec
0209
0210 #ifdef COPY_SORTED_HITS
0211       memcpy(&m_hits[i], &hitv[k], sizeof(Hit));
0212 #endif
0213
0214       if (Config::usePhiQArrays) {
0215         m_hit_infos.emplace_back(hinfos[j]);
0216       }
0217
0218       // Redirect m_binnor.m_ranks[i] to point to external/original index.
0219       m_binnor.m_ranks[i] = k;
0220     }
0221
0222     if (build_original_to_internal_map) {
0223       if (m_max_ext_idx - m_min_ext_idx + 1 > 8 * m_n_hits) {
0224         // If this happens we might:
0225         // a) Use external indices for everything. -- *** We are now. ***
0226         // b) Build these maps for seeding layers only.
0227         // c) Have a flag in hit-on-track that tells us if the hit index has been remapped,
0228         //    essentially, if it is a seed hit. This might be smart anyway.
0229         //    One could use index < -256 or something similar.
0230
0231         printf(
0232             "LayerOfHits::endRegistrationOfHits() original_to_internal index map vector is largish: m_n_hits=%d, "
0233             "map_vector_size=%d\n",
0234             m_n_hits,
0235             m_max_ext_idx - m_min_ext_idx + 1);
0236       }
0237
0238       m_ext_idcs.resize(m_max_ext_idx - m_min_ext_idx + 1);
0239       for (unsigned int i = 0; i < m_n_hits; ++i) {
0240         m_ext_idcs[m_hit_ranks[i] - m_min_ext_idx] = i;
0241       }
0242     }
0243
0244     // We can release m_hit_infos and, if not used, also m_ext_idcs -- and realloc them
0245     // on next beginRegistration().
0246     // If binnor had keep_cons on we could use it for pre-selection in selectHitIndices()
0247     // instead of q and phi arrays -- assuming sufficient precision can be achieved..
0248   }
0249
0250   void LayerOfHits::printBins() {
0251     for (bin_index_t qb = 0; qb <= m_ax_eta.m_last_N_bin; ++qb) {
0252       printf("%c bin %d\n", is_barrel() ? 'Z' : 'R', qb);
0253       for (bin_index_t pb = 0; pb <= m_ax_phi.m_last_N_bin; ++pb) {
0254         if (pb % 8 == 0)
0255           printf(" Phi %4d: ", pb);
0256         auto content = m_binnor.get_content(pb, qb);
0257         printf("%5d,%4d   %s", content.first, content.count, ((pb + 1) % 8 == 0) ? "\n" : "");
0258       }
0259     }
0260   }
0261
0262   //==============================================================================
0263   // EventOfHits
0264   //==============================================================================
0265
0266   EventOfHits::EventOfHits(const TrackerInfo &trk_inf) : m_n_layers(trk_inf.n_layers()) {
0267     m_layers_of_hits.reserve(trk_inf.n_layers());
0268     for (int ii = 0; ii < trk_inf.n_layers(); ++ii) {
0269       m_layers_of_hits.emplace_back(LayerOfHits::Initializator(trk_inf.layer(ii)));
0270     }
0271   }
0272
0273 }  // end namespace mkfit