diff --git a/mkFit/HitStructures.cc b/mkFit/HitStructures.cc index ef362347..0299d83c 100644 --- a/mkFit/HitStructures.cc +++ b/mkFit/HitStructures.cc @@ -39,39 +39,45 @@ void LayerOfHits::SetupLayer(const LayerInfo &li) m_layer_info = &li; - if (is_barrel()) setup_bins(li.m_zmin, li.m_zmax, li.m_q_bin); + m_is_barrel = m_layer_info->is_barrel(); + + if (m_is_barrel) setup_bins(li.m_zmin, li.m_zmax, li.m_q_bin); else setup_bins(li.m_rin, li.m_rout, li.m_q_bin); } //============================================================================== -void LayerOfHits::SuckInHits(const HitVec &hitv) +/* +void detect_q_min_max(const HitVec &hitv) { - // This is now in SetupLayer() - // // should be layer dependant - // float dq = 0.5; - // // should know this from geom. - // //m_qmin = 1000; - // //m_qmax = -1000; - // for (auto const &h : hitv) - // { - // if (h.q() < m_qmin) m_qmin = h.q(); - // if (h.q() > m_qmax) m_qmax = h.q(); - // } - // printf("LoH::SuckInHits qmin=%f, qmax=%f", m_qmin, m_qmax); - // float nmin = std::floor(m_qmin / dq); - // float nmax = std::ceil (m_qmax / dq); - // m_qmin = dq * nmin; - // m_qmax = dq * nmax; - // int nq = nmax - nmin; - // int nqh = nq / 2; - // m_fq = 1.0f / dq; // qbin = (qhit - m_qmin) * m_fq; - // printf(" -> qmin=%f, qmax=%f, nq=%d, fq=%f\n", m_qmin, m_qmax, nq, m_fq); + float dq = 0.5; + // should know this from geom. + //m_qmin = 1000; + //m_qmax = -1000; + for (auto const &h : hitv) + { + if (h.q() < m_qmin) m_qmin = h.q(); + if (h.q() > m_qmax) m_qmax = h.q(); + } + printf("LoH::SuckInHits qmin=%f, qmax=%f", m_qmin, m_qmax); + float nmin = std::floor(m_qmin / dq); + float nmax = std::ceil (m_qmax / dq); + m_qmin = dq * nmin; + m_qmax = dq * nmax; + int nq = nmax - nmin; + int nqh = nq / 2; + m_fq = 1.0f / dq; // qbin = (qhit - m_qmin) * m_fq; + printf(" -> qmin=%f, qmax=%f, nq=%d, fq=%f\n", m_qmin, m_qmax, nq, m_fq); +} +*/ +void LayerOfHits::SuckInHits(const HitVec &hitv) +{ assert (m_nq > 0 && "SetupLayer() was not called."); - const int size = hitv.size(); - const bool is_brl = is_barrel(); + const int size = hitv.size(); + + m_ext_hits = & hitv; #ifdef COPY_SORTED_HITS if (m_capacity < size) @@ -83,48 +89,39 @@ void LayerOfHits::SuckInHits(const HitVec &hitv) if (Config::usePhiQArrays) { + m_hit_phis.resize(size); m_hit_qs.resize(size); + m_hit_infos.resize(size); } - m_hit_phis.resize(size); + m_qphifines.resize(size); - struct HitInfo + for (int i = 0; i < size; ++i) { - float phi; - float q; - uint16_t qbin; - uint16_t phibin; - }; - - std::vector ha(size); - std::vector hit_qphiFines(size); + const Hit &h = hitv[i]; - { - for (int i = 0; i < size; ++i) + HitInfo hi = { h.phi(), m_is_barrel ? h.z() : h.r() }; + + m_qphifines[i] = GetPhiBinFine(hi.phi) + (GetQBinChecked(hi.q) << 16); + + if (Config::usePhiQArrays) { - auto const& h = hitv[i]; - - HitInfo &hi = ha[i]; - // N.1.a For phi in [-pi, pi): squashPhiMinimal(h.phi()); Apparently atan2 can round the wrong way. - hi.phi = h.phi(); - hi.q = is_brl ? h.z() : h.r(); - hi.qbin = std::max(std::min(static_cast((hi.q - m_qmin) * m_fq), m_nq - 1), 0); - // N.1.b ha[j].phi is returned by atan2 and can be rounded the wrong way, resulting in bin -1 or Config::m_nphi - hi.phibin = GetPhiBin(hi.phi) & m_phi_mask; - hit_qphiFines[i] = GetPhiBinFine(hi.phi) + (hi.qbin << 16); - }//i < size + m_hit_infos[i] = hi; + } } - RadixSort sort; - sort.Sort(&hit_qphiFines[0], size); + operator delete [] (m_hit_ranks); + { + RadixSort sort; + sort.Sort(&m_qphifines[0], size, RADIX_UNSIGNED); + m_hit_ranks = sort.RelinquishRanks(); + } - int curr_qphi = -1; + int curr_qphi = -1; empty_q_bins(0, m_nq, 0); - for (uint16_t i = 0; i < size; ++i) + for (int i = 0; i < size; ++i) { - int j = sort.GetRanks()[i]; - - // Could fix the mis-sorts. Set ha size to size + 1 and fake last entry to avoid ifs. + int j = m_hit_ranks[i]; #ifdef COPY_SORTED_HITS memcpy(&m_hits[i], &hitv[j], sizeof(Hit)); @@ -132,15 +129,17 @@ void LayerOfHits::SuckInHits(const HitVec &hitv) if (Config::usePhiQArrays) { - m_hit_phis[i] = ha[j].phi; - m_hit_qs [i] = ha[j].q; + m_hit_phis[i] = m_hit_infos[j].phi; + m_hit_qs [i] = m_hit_infos[j].q; } - const int phi_bin = ha[j].phibin; - const int q_bin = ha[j].qbin; + // Combined q-phi bin with fine part masked off + const int jqphi = m_qphifines[j] & m_phi_fine_xmask; + + const int phi_bin = (jqphi & m_phi_mask_fine) >> m_phi_bits_shift; + const int q_bin = jqphi >> 16; // Fill the bin info - const int jqphi = hit_qphiFines[j] & m_phi_fine_mask; if (jqphi != curr_qphi) { m_phi_bin_infos[q_bin][phi_bin] = {i, i}; @@ -150,12 +149,6 @@ void LayerOfHits::SuckInHits(const HitVec &hitv) m_phi_bin_infos[q_bin][phi_bin].second++; } -#ifndef COPY_SORTED_HITS - operator delete [] (m_hit_ranks); - m_hit_ranks = sort.RelinquishRanks(); - m_ext_hits = & hitv; -#endif - // Check for mis-sorts due to lost precision (not really important). // float phi_prev = 0; // int bin_prev = -1; @@ -185,6 +178,106 @@ void LayerOfHits::SuckInHits(const HitVec &hitv) // } } +//============================================================================== + + +void LayerOfHits::BeginRegistrationOfHits(const HitVec &hitv) +{ + assert (m_nq > 0 && "SetupLayer() was not called."); + + m_ext_hits = &hitv; + + m_hit_infos.clear(); + m_qphifines.clear(); + m_ext_idcs.clear(); +} + +void LayerOfHits::RegisterHit(int idx) +{ + const Hit &h = (*m_ext_hits)[idx]; + + m_ext_idcs.push_back(idx); + + HitInfo hi = { h.phi(), m_is_barrel ? h.z() : h.r() }; + + m_qphifines.push_back( GetPhiBinFine(hi.phi) + (GetQBinChecked(hi.q) << 16) ); + + if (Config::usePhiQArrays) + { + m_hit_infos.emplace_back(hi); + } +} + +void LayerOfHits::EndRegistrationOfHits() +{ + const int size = m_ext_idcs.size(); + + // radix + operator delete [] (m_hit_ranks); + { + RadixSort sort; + sort.Sort(&m_qphifines[0], size, RADIX_UNSIGNED); + m_hit_ranks = sort.RelinquishRanks(); + } + + // copy q/phi + +#ifdef COPY_SORTED_HITS + if (m_capacity < size) + { + free_hits(); + alloc_hits(1.02 * size); + } +#endif + + if (Config::usePhiQArrays) + { + m_hit_phis.resize(size); + m_hit_qs.resize(size); + } + + int curr_qphi = -1; + empty_q_bins(0, m_nq, 0); + + for (int i = 0; i < size; ++i) + { + int j = m_hit_ranks[i]; // index in intermediate + int k = m_ext_idcs[j]; // index in external hit_vec + +#ifdef COPY_SORTED_HITS + memcpy(&m_hits[i], &hitv[k], sizeof(Hit)); +#endif + + if (Config::usePhiQArrays) + { + m_hit_phis[i] = m_hit_infos[j].phi; + m_hit_qs [i] = m_hit_infos[j].q; + } + + // Combined q-phi bin with fine part masked off + const int jqphi = m_qphifines[j] & m_phi_fine_xmask; + + const int phi_bin = (jqphi & m_phi_mask_fine) >> m_phi_bits_shift; + const int q_bin = jqphi >> 16; + + // Fill the bin info + if (jqphi != curr_qphi) + { + m_phi_bin_infos[q_bin][phi_bin] = {i, i}; + curr_qphi = jqphi; + } + + m_phi_bin_infos[q_bin][phi_bin].second++; + + // Ranks[i] will never be used again ... use it to point to external index. + m_hit_ranks[i] = k; + } + + // Matti: We can release m_hit_infos, m_ext_idcs, m_qphifines -- and realloc on next BeginInput. +} + +//============================================================================== + void LayerOfHits::SelectHitIndices(float q, float phi, float dq, float dphi, std::vector& idcs, bool isForSeeding, bool dump) { // Sanitizes q, dq and dphi. phi is expected to be in -pi, pi. diff --git a/mkFit/HitStructures.h b/mkFit/HitStructures.h index f3cc1c7d..2730b3a8 100644 --- a/mkFit/HitStructures.h +++ b/mkFit/HitStructures.h @@ -76,6 +76,16 @@ class LayerOfHits const HitVec *m_ext_hits; #endif + // Stuff needed during setup + struct HitInfo + { + float phi; + float q; + }; + std::vector m_hit_infos; + std::vector m_qphifines; + std::vector m_ext_idcs; + public: const LayerInfo *m_layer_info = 0; vecvecPhiBinInfo_t m_phi_bin_infos; @@ -84,14 +94,15 @@ class LayerOfHits float m_qmin, m_qmax, m_fq; int m_nq = 0; + bool m_is_barrel; - int layer_id() const { return m_layer_info->m_layer_id; } - bool is_barrel() const { return m_layer_info->is_barrel(); } - bool is_endcap() const { return ! m_layer_info->is_barrel(); } + int layer_id() const { return m_layer_info->m_layer_id; } + bool is_barrel() const { return m_is_barrel; } + bool is_endcap() const { return ! m_is_barrel; } int bin_index(int q, int p) const { return q*Config::m_nphi + p; } PhiBinInfo_t operator[](int i) const { - int q = i/Config::m_nphi; + int q = i / Config::m_nphi; int p = i % Config::m_nphi; return m_phi_bin_infos[q][p]; } @@ -141,11 +152,11 @@ class LayerOfHits static constexpr float m_fphi = Config::m_nphi / Config::TwoPI; static constexpr int m_phi_mask = 0x7f; static constexpr int m_phi_bits = 7; - static constexpr float m_fphi_fine = 1024 / Config::TwoPI; - static constexpr int m_phi_mask_fine = 0x3ff; - static constexpr int m_phi_bits_fine = 10;//can't be more than 16 + static constexpr float m_fphi_fine = 1024 / Config::TwoPI; + static constexpr int m_phi_mask_fine = 0x3ff; + static constexpr int m_phi_bits_fine = 10; //can't be more than 16 static constexpr int m_phi_bits_shift = m_phi_bits_fine - m_phi_bits; - static constexpr int m_phi_fine_mask = ~((1 << m_phi_bits_shift) - 1); + static constexpr int m_phi_fine_xmask = ~((1 << m_phi_bits_shift) - 1); protected: @@ -165,12 +176,14 @@ class LayerOfHits void setup_bins(float qmin, float qmax, float dq); - void set_phi_bin(int q_bin, int phi_bin, uint16_t &hit_count, uint16_t &hits_in_bin) - { - m_phi_bin_infos[q_bin][phi_bin] = { hit_count, hit_count + hits_in_bin }; - hit_count += hits_in_bin; - hits_in_bin = 0; - } + + // Not used. + // void set_phi_bin(int q_bin, int phi_bin, uint16_t &hit_count, uint16_t &hits_in_bin) + // { + // m_phi_bin_infos[q_bin][phi_bin] = { hit_count, hit_count + hits_in_bin }; + // hit_count += hits_in_bin; + // hits_in_bin = 0; + // } void empty_phi_bins(int q_bin, int phi_bin_1, int phi_bin_2, uint16_t hit_count) { @@ -195,37 +208,43 @@ class LayerOfHits { #ifdef COPY_SORTED_HITS free_hits(); -#else - operator delete [] (m_hit_ranks); #endif + operator delete [] (m_hit_ranks); } void SetupLayer(const LayerInfo &li); void Reset() {} - float NormalizeQ(float q) const { if (q < m_qmin) return m_qmin; if (q > m_qmax) return m_qmax; return q; } + float NormalizeQ(float q) const { return std::clamp(q, m_qmin, m_qmax); } int GetQBin(float q) const { return (q - m_qmin) * m_fq; } - int GetQBinChecked(float q) const { int qb = GetQBin(q); if (qb < 0) qb = 0; else if (qb >= m_nq) qb = m_nq - 1; return qb; } + int GetQBinChecked(float q) const { return std::clamp(GetQBin(q), 0, m_nq - 1); } // if you don't pass phi in (-pi, +pi), mask away the upper bits using m_phi_mask or use the Checked version. int GetPhiBinFine(float phi) const { return std::floor(m_fphi_fine * (phi + Config::PI)); } - int GetPhiBin(float phi) const { return GetPhiBinFine(phi)>>m_phi_bits_shift; } + int GetPhiBin (float phi) const { return GetPhiBinFine(phi)>>m_phi_bits_shift; } int GetPhiBinChecked(float phi) const { return GetPhiBin(phi) & m_phi_mask; } const vecPhiBinInfo_t& GetVecPhiBinInfo(float q) const { return m_phi_bin_infos[GetQBin(q)]; } + // Get in all hits from given hit-vec void SuckInHits(const HitVec &hitv); + // Use external hit-vec and only use hits that are passed to me. + void BeginRegistrationOfHits(const HitVec &hitv); + void RegisterHit(int idx); + void EndRegistrationOfHits(); + + // Use this to remap internal hit index to external one. + int GetOriginalHitIndex(int i) const { return m_hit_ranks[i]; } + #ifdef COPY_SORTED_HITS - int GetHitIndex(int i) const { return i; } const Hit& GetHit(int i) const { return m_hits[i]; } const Hit* GetHitArray() const { return m_hits; } #else - int GetHitIndex(int i) const { return m_hit_ranks[i]; } const Hit& GetHit(int i) const { return (*m_ext_hits)[m_hit_ranks[i]]; } const Hit* GetHitArray() const { return & (*m_ext_hits)[0]; } #endif @@ -257,6 +276,13 @@ class EventOfHits void SuckInHits(int layer, const HitVec &hitv) { m_layers_of_hits[layer].SuckInHits(hitv); + /* + int nh = hitv.size(); + auto &loh = m_layers_of_hits[layer]; + loh.BeginRegistrationOfHits(hitv); + for (int i = 0; i < nh; ++i) loh.RegisterHit(i); + loh.EndRegistrationOfHits(); + */ } }; diff --git a/mkFit/MkBuilder.cc b/mkFit/MkBuilder.cc index f6e285b5..c5fffac9 100644 --- a/mkFit/MkBuilder.cc +++ b/mkFit/MkBuilder.cc @@ -911,6 +911,32 @@ void MkBuilder::remap_track_hits(TrackVec & tracks) } } } + + // Matti ... this can now be just something like this (not tested): + /* + for (auto&& track : tracks) + { + for (int i = 0; i < track.nTotalHits(); ++i) + { + int hitidx = track.getHitIdx(i); + int hitlyr = track.getHitLyr(i); + if (hitidx >= 0) + { + const auto & loh = m_event_of_hits.m_layers_of_hits[hitlyr]; + track.setHitIdx(i, loh.GetOriginalHitIndex(hitidx)); + } + } + } + */ + // Note that the indices after this are correct CMSSW "large-vector" indices. + // So no hit/layer mapping code in CMSSW producer is needed. + // + // We could really store original indices into HitOnTrack.index + // from the start. One just needs to be careful when getting hits in backwards fit, + // to use (a currently non-existent) LayerOfHits::GetHitWithOriginalIndex(hot.index) + // + // Further, somebody should walk over quite a bit of validation code that is + // rather involved doing such remappings. } //------------------------------------------------------------------------------