diff --git a/RecoTracker/MkFit/plugins/MkFitGeometryESProducer.cc b/RecoTracker/MkFit/plugins/MkFitGeometryESProducer.cc index 24ef8b2bd7adf..b62e9aeb1f7b0 100644 --- a/RecoTracker/MkFit/plugins/MkFitGeometryESProducer.cc +++ b/RecoTracker/MkFit/plugins/MkFitGeometryESProducer.cc @@ -198,7 +198,7 @@ void MkFitGeometryESProducer::fillShapeAndPlacement(const GeomDet *det, DetId detid = det->geographicalId(); float xy[4][2]; - float dz; + float half_length, dz; const Bounds *b = &((det->surface()).bounds()); if (const TrapezoidalPlaneBounds *b2 = dynamic_cast(b)) { @@ -212,6 +212,7 @@ void MkFitGeometryESProducer::fillShapeAndPlacement(const GeomDet *det, xy[2][1] = par[3]; xy[3][0] = par[0]; xy[3][1] = -par[3]; + half_length = par[3]; dz = par[2]; #ifdef DUMP_MKF_GEO @@ -229,6 +230,7 @@ void MkFitGeometryESProducer::fillShapeAndPlacement(const GeomDet *det, xy[2][1] = dy; xy[3][0] = dx; xy[3][1] = -dy; + half_length = dy; dz = b2->thickness() * 0.5; // half thickness #ifdef DUMP_MKF_GEO @@ -280,7 +282,8 @@ void MkFitGeometryESProducer::fillShapeAndPlacement(const GeomDet *det, const auto &p = det->position(); auto z = det->rotation().z(); auto x = det->rotation().x(); - layer_info.register_module({{p.x(), p.y(), p.z()}, {z.x(), z.y(), z.z()}, {x.x(), x.y(), x.z()}, detid.rawId()}); + layer_info.register_module( + {{p.x(), p.y(), p.z()}, {z.x(), z.y(), z.z()}, {x.x(), x.y(), x.z()}, half_length, detid.rawId()}); // Set some layer parameters (repeatedly, would require hard-coding otherwise) layer_info.set_subdet(detid.subdetId()); layer_info.set_is_pixel(detid.subdetId() <= 2); diff --git a/RecoTracker/MkFitCMS/standalone/Makefile b/RecoTracker/MkFitCMS/standalone/Makefile index 3cf3b83cec91d..fdb15758af987 100644 --- a/RecoTracker/MkFitCMS/standalone/Makefile +++ b/RecoTracker/MkFitCMS/standalone/Makefile @@ -6,13 +6,18 @@ CMS_DIR := ${SRCDIR}/RecoTracker/MkFitCMS LIB_CMS := ../libMicCMS.so MAIN := ../mkFit +MAIN_DEPS := ${LIB_CMS} +MAIN_LIBS := -lMicCore -lMicCMS WRMEMF := ../writeMemoryFile WMF_DICT_PCM := ../WriteMemFileDict_rdict.pcm SHELL_DICT_PCM := ../ShellDict_rdict.pcm ROOTOUT := WriteMemFileDict.cc ShellDict.cc TGTS := ${LIB_CMS} ${MAIN} + ifdef WITH_ROOT +MAIN_DEPS += ../libMicRntDump.so +MAIN_LIBS += -lMicRntDump TGTS += ${WRMEMF} ${WMF_DICT_PCM} ${SHELL_DICT_PCM} endif @@ -21,18 +26,19 @@ endif all: ${TGTS} SRCS := $(wildcard ${CMS_DIR}/src/*.cc) $(wildcard ${SACMS}/*.cc) + ifdef WITH_ROOT SRCS += ${SACMS}/tkNtuple/WriteMemoryFile.cc -WriteMemFileDict.cc: ${SACMS}/tkNtuple/DictsLinkDef.h +WriteMemFileDict.cc ${WMF_DICT_PCM} &: ${SACMS}/tkNtuple/DictsLinkDef.h rootcling -I=${SRCDIR} -f WriteMemFileDict.cc $< -${WMF_DICT_PCM}: WriteMemFileDict.cc mv WriteMemFileDict_rdict.pcm ${WMF_DICT_PCM} + SRCS += ShellDict.cc -ShellDict.cc: ${SACMS}/Shell.h ${SACMS}/ShellLinkDef.h +ShellDict.cc ${SHELL_DICT_PCM} &: ${SACMS}/Shell.h ${SACMS}/ShellLinkDef.h rootcling -I=${SRCDIR} -f ShellDict.cc ${SACMS}/Shell.h ${SACMS}/ShellLinkDef.h -${SHELL_DICT_PCM}: ShellDict.cc mv ShellDict_rdict.pcm ${SHELL_DICT_PCM} endif + SRCB := $(notdir ${SRCS}) DEPS := $(SRCB:.cc=.d) OBJS := $(SRCB:.cc=.o) @@ -60,11 +66,11 @@ ${LIB_CMS}: ${CMS_OBJS} @mkdir -p $(@D) ${CXX} ${CXXFLAGS} ${VEC_HOST} ${CMS_OBJS} -shared -o $@ ${LDFLAGS_HOST} ${LDFLAGS} -${MAIN}: ${LIB_CMS} mkFit.o - ${CXX} ${CXXFLAGS} ${VEC_HOST} ${LDFLAGS} mkFit.o -o $@ ${LDFLAGS_HOST} -ltbb -L.. -lMicCore -lMicCMS -Wl,-rpath=. +${MAIN}: ${MAIN_DEPS} mkFit.o + ${CXX} ${CXXFLAGS} ${VEC_HOST} ${LDFLAGS} mkFit.o -o $@ ${LDFLAGS_HOST} -ltbb -L.. ${MAIN_LIBS} -Wl,-rpath=. -${WRMEMF}: WriteMemoryFile.o WriteMemFileDict.o - ${CXX} ${CXXFLAGS} ${LDFLAGS} $^ -o $@ ${LDFLAGS_HOST} -ltbb -L.. -lMicCore -Wl,-rpath=. +${WRMEMF}: ${MAIN_DEPS} WriteMemoryFile.o WriteMemFileDict.o + ${CXX} ${CXXFLAGS} ${LDFLAGS} $^ -o $@ ${LDFLAGS_HOST} -ltbb -L.. ${MAIN_LIBS} -Wl,-rpath=. ${OBJS}: %.o: %.cc %.d ${CXX} ${CPPFLAGS} ${CXXFLAGS} ${VEC_HOST} -c -o $@ $< @@ -77,6 +83,7 @@ echo: @echo SRCS = ${SRCS} @echo DEPS = ${DEPS} @echo OBJS = ${OBJS} + @echo MAIN_LIBS = ${MAIN_LIBS} echo_cc_defs: ${CXX} -dM -E -mavx2 - < /dev/null diff --git a/RecoTracker/MkFitCMS/standalone/MkStandaloneSeqs.cc b/RecoTracker/MkFitCMS/standalone/MkStandaloneSeqs.cc index 4f9ae7d3c490c..8c015cf91a23c 100644 --- a/RecoTracker/MkFitCMS/standalone/MkStandaloneSeqs.cc +++ b/RecoTracker/MkFitCMS/standalone/MkStandaloneSeqs.cc @@ -159,7 +159,7 @@ namespace mkfit { if (mctrk < 0 || static_cast(mctrk) >= event->simTracks_.size()) { ++m_cnt_nomc; - dprint("XX bad track idx " << mctrk << ", orig label was " << label); + dprintf("XX bad track idx %d, orig label was %d\n", mctrk, label); } else { auto &simtrack = event->simTracks_[mctrk]; pTmc = simtrack.pT(); diff --git a/RecoTracker/MkFitCMS/standalone/Shell.cc b/RecoTracker/MkFitCMS/standalone/Shell.cc index 2f4c1aae539af..ad4fa12c00827 100644 --- a/RecoTracker/MkFitCMS/standalone/Shell.cc +++ b/RecoTracker/MkFitCMS/standalone/Shell.cc @@ -200,6 +200,8 @@ namespace mkfit { s.sortHitsByLayer(); // sort seed hits for the matched hits (I hope it works here) } + m_event->setCurrentSeedTracks(seeds); + builder.find_tracks_load_seeds(seeds, do_seed_clean); builder.findTracksCloneEngine(); @@ -265,6 +267,8 @@ namespace mkfit { printf("Shell::ProcessEvent post remove-duplicates: %d comb-cands\n", (int) out_tracks.size()); + m_event->resetCurrentSeedTracks(); + builder.end_event(); } diff --git a/RecoTracker/MkFitCMS/standalone/buildtestMPlex.cc b/RecoTracker/MkFitCMS/standalone/buildtestMPlex.cc index 6872274f24436..01995e17aed47 100644 --- a/RecoTracker/MkFitCMS/standalone/buildtestMPlex.cc +++ b/RecoTracker/MkFitCMS/standalone/buildtestMPlex.cc @@ -110,9 +110,8 @@ namespace mkfit { const bool validation_on = (Config::sim_val || Config::quality_val); + TrackVec seeds1; if (validation_on) { - TrackVec seeds1; - unsigned int algorithms[] = {4}; //only initialStep for (auto const &s : ev.seedTracks_) { @@ -187,6 +186,10 @@ namespace mkfit { // ev.print_tracks(ev.candidateTracks_, true); + if (validation_on) { + ev.seedTracks_.swap(seeds1); + } + return time; } @@ -199,9 +202,8 @@ namespace mkfit { const bool validation_on = (Config::sim_val || Config::quality_val); + TrackVec seeds1; if (validation_on) { - TrackVec seeds1; - unsigned int algorithms[] = {4}; //only initialStep for (auto const &s : ev.seedTracks_) { @@ -226,6 +228,8 @@ namespace mkfit { bool seeds_sorted = false; // CCCC builder.PrepareSeeds(); + ev.setCurrentSeedTracks(ev.seedTracks_); + ev.simLabelForCurrentSeed(0); builder.find_tracks_load_seeds(ev.seedTracks_, seeds_sorted); @@ -248,6 +252,7 @@ namespace mkfit { check_nan_n_silly_candidates(ev); // first store candidate tracks + ev.candidateTracks_.clear(); builder.export_best_comb_cands(ev.candidateTracks_); job.switch_to_backward(); @@ -272,10 +277,16 @@ namespace mkfit { StdSeq::root_val(&ev); } + ev.resetCurrentSeedTracks(); + builder.end_event(); // ev.print_tracks(ev.candidateTracks_, true); + if (validation_on) { + ev.seedTracks_.swap(seeds1); + } + return time; } @@ -288,9 +299,8 @@ namespace mkfit { const bool validation_on = (Config::sim_val || Config::quality_val); + TrackVec seeds1; if (validation_on) { - TrackVec seeds1; - unsigned int algorithms[] = {4}; //only initialStep for (auto const &s : ev.seedTracks_) { @@ -315,6 +325,7 @@ namespace mkfit { bool seeds_sorted = false; // CCCC builder.PrepareSeeds(); + ev.setCurrentSeedTracks(ev.seedTracks_); builder.find_tracks_load_seeds(ev.seedTracks_, seeds_sorted); @@ -337,6 +348,7 @@ namespace mkfit { check_nan_n_silly_candidates(ev); // first store candidate tracks - needed for BH backward fit and root_validation + ev.candidateTracks_.clear(); builder.export_best_comb_cands(ev.candidateTracks_); job.switch_to_backward(); @@ -365,10 +377,16 @@ namespace mkfit { StdSeq::root_val(&ev); } + ev.resetCurrentSeedTracks(); + builder.end_event(); // ev.print_tracks(ev.candidateTracks_, true); + if (validation_on) { + ev.seedTracks_.swap(seeds1); + } + return time; } @@ -458,6 +476,8 @@ namespace mkfit { // Add protection in case no seeds are found for iteration if (seeds.size() <= 0) continue; + ev.setCurrentSeedTracks(seeds); + ev.simLabelForCurrentSeed(0); builder.find_tracks_load_seeds(seeds, do_seed_clean); @@ -546,6 +566,7 @@ namespace mkfit { builder.export_tracks(ev.fitTracks_); } + ev.resetCurrentSeedTracks(); builder.end_event(); } diff --git a/RecoTracker/MkFitCMS/standalone/mkFit.cc b/RecoTracker/MkFitCMS/standalone/mkFit.cc index 889a7d5d63d69..7c5370b538134 100644 --- a/RecoTracker/MkFitCMS/standalone/mkFit.cc +++ b/RecoTracker/MkFitCMS/standalone/mkFit.cc @@ -16,6 +16,7 @@ #ifndef NO_ROOT #include "RecoTracker/MkFitCore/standalone/Validation.h" +#include "RecoTracker/MkFitCore/standalone/RntDumper/RntDumper.h" #endif //#define DEBUG @@ -127,7 +128,7 @@ namespace { bool g_run_fit_std = false; - bool g_run_build_all = true; + bool g_run_build_default = true; bool g_run_build_cmssw = false; bool g_run_build_bh = false; bool g_run_build_std = false; @@ -333,14 +334,14 @@ void test_standard() { int maxLayer_thisthread = 0; for (int b = 0; b < Config::finderReportBestOutOfN; ++b) { t_cur[0] = 0; // t_cur[0] = (g_run_fit_std) ? runFittingTestPlex(ev, plex_tracks) : 0; - t_cur[1] = (g_run_build_all || g_run_build_bh) ? runBuildingTestPlexBestHit(ev, eoh, mkb) : 0; - t_cur[3] = (g_run_build_all || g_run_build_ce) ? runBuildingTestPlexCloneEngine(ev, eoh, mkb) : 0; - if (g_run_build_all || g_run_build_mimi) + t_cur[1] = (g_run_build_bh) ? runBuildingTestPlexBestHit(ev, eoh, mkb) : 0; + t_cur[2] = (g_run_build_default || g_run_build_std) ? runBuildingTestPlexStandard(ev, eoh, mkb) : 0; + t_cur[3] = (g_run_build_default || g_run_build_ce) ? runBuildingTestPlexCloneEngine(ev, eoh, mkb) : 0; + if (g_run_build_mimi) t_cur_iter = runBtpCe_MultiIter(ev, eoh, mkb, Config::nItersCMSSW); - t_cur[4] = (g_run_build_all || g_run_build_mimi) ? t_cur_iter[Config::nItersCMSSW] : 0; - if (g_run_build_all || g_run_build_cmssw) + t_cur[4] = (g_run_build_mimi) ? t_cur_iter[Config::nItersCMSSW] : 0; + if (g_run_build_cmssw) runBuildingTestPlexDumbCMSSW(ev, eoh, mkb); - t_cur[2] = (g_run_build_all || g_run_build_std) ? runBuildingTestPlexStandard(ev, eoh, mkb) : 0; if (g_run_build_ce || g_run_build_mimi) { ncands_thisthread = mkb.total_cands(); auto const& ln = mkb.max_hits_layer(eoh); @@ -387,7 +388,7 @@ void test_standard() { if (evt > 0) for (int i = 0; i < NT; ++i) t_skip[i] += t_best[i]; - if (g_run_build_all || g_run_build_mimi) { + if (g_run_build_mimi) { for (int i = 0; i < Config::nItersCMSSW; ++i) t_sum_iter[i] += t_cur_iter[i]; if (evt > 0) @@ -427,7 +428,7 @@ void test_standard() { maxHits_all.load(), maxLayer_all.load()); //fflush(stdout); - if (g_run_build_all || g_run_build_mimi) { + if (g_run_build_mimi) { printf("================================================================\n"); for (int i = 0; i < Config::nItersCMSSW; ++i) std::cout << " Iteration " << i << " build time = " << t_sum_iter[i] << " \n"; @@ -537,7 +538,6 @@ int main(int argc, const char* argv[]) { "\n\n" "FittingTestMPlex options\n\n" " --fit-std run standard fitting test (def: %s)\n" - " --fit-std-only run only standard fitting test (def: %s)\n" " --cf-fitting enable conformal fit before fitting tracks to get initial estimate of track " "parameters and errors (def: %s)\n" " --fit-val enable ROOT based validation for fittingMPlex (def: %s)\n" @@ -697,16 +697,14 @@ int main(int argc, const char* argv[]) { Config::numHitsPerTask, b2a(g_run_fit_std), - b2a(g_run_fit_std && - !(g_run_build_all || g_run_build_cmssw || g_run_build_bh || g_run_build_std || g_run_build_ce)), b2a(Config::cf_fitting), b2a(Config::fit_val), - b2a(g_run_build_all || g_run_build_cmssw), - b2a(g_run_build_all || g_run_build_bh), - b2a(g_run_build_all || g_run_build_std), - b2a(g_run_build_all || g_run_build_ce), - b2a(g_run_build_all || g_run_build_mimi), + b2a(g_run_build_cmssw), + b2a(g_run_build_bh), + b2a(g_run_build_default || g_run_build_std), + b2a(g_run_build_default || g_run_build_ce), + b2a(g_run_build_mimi), getOpt(Config::seedInput, g_seed_opts).c_str(), getOpt(Config::seedCleaning, g_clean_opts).c_str(), @@ -835,47 +833,26 @@ int main(int argc, const char* argv[]) { next_arg_or_die(mArgs, i); Config::numHitsPerTask = atoi(i->c_str()); } else if (*i == "--fit-std") { + g_run_build_default = false; g_run_fit_std = true; - } else if (*i == "--fit-std-only") { - g_run_fit_std = true; - g_run_build_all = false; - g_run_build_bh = false; - g_run_build_std = false; - g_run_build_ce = false; } else if (*i == "--cf-fitting") { Config::cf_fitting = true; } else if (*i == "--fit-val") { Config::fit_val = true; } else if (*i == "--build-cmssw") { - g_run_build_all = false; + g_run_build_default = false; g_run_build_cmssw = true; - g_run_build_bh = false; - g_run_build_std = false; - g_run_build_ce = false; } else if (*i == "--build-bh") { - g_run_build_all = false; - g_run_build_cmssw = false; + g_run_build_default = false; g_run_build_bh = true; - g_run_build_std = false; - g_run_build_ce = false; } else if (*i == "--build-std") { - g_run_build_all = false; - g_run_build_cmssw = false; - g_run_build_bh = false; + g_run_build_default = false; g_run_build_std = true; - g_run_build_ce = false; } else if (*i == "--build-ce") { - g_run_build_all = false; - g_run_build_cmssw = false; - g_run_build_bh = false; - g_run_build_std = false; + g_run_build_default = false; g_run_build_ce = true; } else if (*i == "--build-mimi") { - g_run_build_all = false; - g_run_build_cmssw = false; - g_run_build_bh = false; - g_run_build_std = false; - g_run_build_ce = false; + g_run_build_default = false; g_run_build_mimi = true; if (Config::nItersCMSSW == 0) Config::nItersCMSSW = 3; @@ -1047,6 +1024,10 @@ int main(int argc, const char* argv[]) { test_standard(); } +#ifndef NO_ROOT + RntDumper::FinalizeAll(); +#endif + // clang-format on return 0; diff --git a/RecoTracker/MkFitCore/BuildFile.xml b/RecoTracker/MkFitCore/BuildFile.xml index 1c172b635b419..5417d592d4bd1 100644 --- a/RecoTracker/MkFitCore/BuildFile.xml +++ b/RecoTracker/MkFitCore/BuildFile.xml @@ -1,6 +1,7 @@ + diff --git a/RecoTracker/MkFitCore/interface/Config.h b/RecoTracker/MkFitCore/interface/Config.h index b5e3695f90523..ac4dfe1277b0a 100644 --- a/RecoTracker/MkFitCore/interface/Config.h +++ b/RecoTracker/MkFitCore/interface/Config.h @@ -10,7 +10,7 @@ namespace mkfit { constexpr float PIOver4 = Const::PI / 4.0f; constexpr float PI3Over4 = 3.0f * Const::PI / 4.0f; constexpr float InvPI = 1.0f / Const::PI; - constexpr float sol = 0.299792458; // speed of light in nm/s + constexpr float sol = 0.299792458; // speed of light in m/ns // NAN and silly track parameter tracking options constexpr bool nan_etc_sigs_enable = false; diff --git a/RecoTracker/MkFitCore/interface/HitStructures.h b/RecoTracker/MkFitCore/interface/HitStructures.h index 13299358a9a33..f08d942c18cd5 100644 --- a/RecoTracker/MkFitCore/interface/HitStructures.h +++ b/RecoTracker/MkFitCore/interface/HitStructures.h @@ -52,7 +52,7 @@ namespace mkfit { // Setup and filling //------------------- - void reset() {} + void reset(); // Get in all hits from given hit-vec void suckInHits(const HitVec& hitv); @@ -82,8 +82,17 @@ namespace mkfit { bool isBinDead(bin_index_t pi, bin_index_t qi) const { return m_dead_bins[qi * m_ax_phi.size_of_N() + pi]; } - float hit_q(unsigned int i) const { return m_hit_qs[i]; } - float hit_phi(unsigned int i) const { return m_hit_phis[i]; } + struct HitInfo { + float phi; + float q; + float q_half_length; + float qbar; + }; + const HitInfo& hit_info(unsigned int i) const { return m_hit_infos[i]; } + float hit_phi(unsigned int i) const { return m_hit_infos[i].phi; } + float hit_q(unsigned int i) const { return m_hit_infos[i].q; } + float hit_q_half_length(unsigned int i) const { return m_hit_infos[i].q_half_length; } + float hit_qbar(unsigned int i) const { return m_hit_infos[i].qbar; } // Use this to map original indices to sorted internal ones. m_ext_idcs needs to be initialized. unsigned int getHitIndexFromOriginal(unsigned int i) const { return m_ext_idcs[i - m_min_ext_idx]; } @@ -146,19 +155,11 @@ namespace mkfit { // Bin information for dead regions std::vector m_dead_bins; - // Cached hit phi and q values to minimize Hit memory access - std::vector m_hit_phis; - std::vector m_hit_qs; - // Geometry / q-binning constants - initialized in setupLayer() const LayerInfo* m_layer_info = nullptr; bool m_is_barrel; - // Data needed during setup - struct HitInfo { - float phi; - float q; - }; + // Cached hit phi, q and qbar values to minimize Hit memory access std::vector m_hit_infos; }; diff --git a/RecoTracker/MkFitCore/interface/MkBuilder.h b/RecoTracker/MkFitCore/interface/MkBuilder.h index 9e68ac703665a..c68b3e7bf80e6 100644 --- a/RecoTracker/MkFitCore/interface/MkBuilder.h +++ b/RecoTracker/MkFitCore/interface/MkBuilder.h @@ -27,7 +27,7 @@ namespace mkfit { class MkBuilder { public: - using insert_seed_foo = void(const Track &, int, int); + using insert_seed_foo = void(const Track &, int, int, int); typedef std::vector> CandIdx_t; diff --git a/RecoTracker/MkFitCore/interface/Track.h b/RecoTracker/MkFitCore/interface/Track.h index d941b615663f0..a0aaa772fd728 100644 --- a/RecoTracker/MkFitCore/interface/Track.h +++ b/RecoTracker/MkFitCore/interface/Track.h @@ -193,11 +193,20 @@ namespace mkfit { void setLabel(int lbl) { label_ = lbl; } bool hasSillyValues(bool dump, bool fix, const char* pref = ""); - bool hasNanNSillyValues() const; + // ------------------------------------------------------------------------ + float d0BeamSpot(const float x_bs, const float y_bs, bool linearize = false) const; + // used for swimming cmssw rec tracks to mkFit position + float swimPhiToR(const float x, const float y) const; + + bool canReachRadius(float R) const; + float maxReachRadius() const; + float zAtR(float R, float* r_reached = nullptr) const; + float rAtZ(float Z) const; + // ------------------------------------------------------------------------ struct Status { @@ -355,8 +364,7 @@ namespace mkfit { // TrackCand //============================================================================== - // TrackCand depends on stuff in mkFit/HitStructures, CombCand in particular, - // so it is declared / implemented there. + // TrackCand defined in TrackStructures.h along with CombCandidate. // class TrackCand : public TrackBase { ... }; @@ -387,15 +395,7 @@ namespace mkfit { Track(const Track& t) : TrackBase(t), hitsOnTrk_(t.hitsOnTrk_) {} - // used for swimming cmssw rec tracks to mkFit position - float swimPhiToR(const float x, const float y) const; - - bool canReachRadius(float R) const; - float maxReachRadius() const; - float zAtR(float R, float* r_reached = nullptr) const; - float rAtZ(float Z) const; - - //this function is very inefficient, use only for debug and validation! + // This function is very inefficient, use only for debug and validation! HitVec hitsVector(const std::vector& globalHitVec) const { HitVec hitsVec; for (int ihit = 0; ihit < Config::nMaxTrkHits; ++ihit) { @@ -420,6 +420,15 @@ namespace mkfit { } } + int mcHitIDofFirstHit(const std::vector& globalHitVec, const MCHitInfoVec& globalMCHitInfo) const { + const HitOnTrack& hot = hitsOnTrk_[0]; + if ((hot.index >= 0) && (static_cast(hot.index) < globalHitVec[hot.layer].size())) { + return globalHitVec[hot.layer][hot.index].mcTrackID(globalMCHitInfo); + } else { + return hot.index; + } + } + // The following 2 (well, 3) funcs to be fixed once we move lastHitIdx_ and nFoundHits_ // out of TrackBase. If we do it. void reserveHits(int nHits) { hitsOnTrk_.reserve(nHits); } diff --git a/RecoTracker/MkFitCore/interface/TrackStructures.h b/RecoTracker/MkFitCore/interface/TrackStructures.h index d4cc07f26a20d..736fbe8847e55 100644 --- a/RecoTracker/MkFitCore/interface/TrackStructures.h +++ b/RecoTracker/MkFitCore/interface/TrackStructures.h @@ -290,13 +290,9 @@ namespace mkfit { m_lastHitIdx_before_bkwsearch(o.m_lastHitIdx_before_bkwsearch), m_nInsideMinusOneHits_before_bkwsearch(o.m_nInsideMinusOneHits_before_bkwsearch), m_nTailMinusOneHits_before_bkwsearch(o.m_nTailMinusOneHits_before_bkwsearch), -#ifdef DUMPHITWINDOW - m_seed_algo(o.m_seed_algo), - m_seed_label(o.m_seed_label), -#endif + m_seed_origin_index(o.m_seed_origin_index), m_hots_size(o.m_hots_size), - m_hots(o.m_hots) { - } + m_hots(o.m_hots) {} // Required for std::swap(). CombCandidate(CombCandidate&& o) @@ -307,10 +303,7 @@ namespace mkfit { m_lastHitIdx_before_bkwsearch(o.m_lastHitIdx_before_bkwsearch), m_nInsideMinusOneHits_before_bkwsearch(o.m_nInsideMinusOneHits_before_bkwsearch), m_nTailMinusOneHits_before_bkwsearch(o.m_nTailMinusOneHits_before_bkwsearch), -#ifdef DUMPHITWINDOW - m_seed_algo(o.m_seed_algo), - m_seed_label(o.m_seed_label), -#endif + m_seed_origin_index(o.m_seed_origin_index), m_hots_size(o.m_hots_size), m_hots(std::move(o.m_hots)) { // This is not needed as we do EOCC::reset() after EOCCS::resize which @@ -332,10 +325,7 @@ namespace mkfit { m_lastHitIdx_before_bkwsearch = o.m_lastHitIdx_before_bkwsearch; m_nInsideMinusOneHits_before_bkwsearch = o.m_nInsideMinusOneHits_before_bkwsearch; m_nTailMinusOneHits_before_bkwsearch = o.m_nTailMinusOneHits_before_bkwsearch; -#ifdef DUMPHITWINDOW - m_seed_algo = o.m_seed_algo; - m_seed_label = o.m_seed_label; -#endif + m_seed_origin_index = o.m_seed_origin_index; m_hots_size = o.m_hots_size; m_hots = std::move(o.m_hots); @@ -378,7 +368,7 @@ namespace mkfit { m_nTailMinusOneHits_before_bkwsearch = -1; } - void importSeed(const Track& seed, const track_score_func& score_func, int region); + void importSeed(const Track& seed, int seed_idx, const track_score_func& score_func, int region); int addHit(const HitOnTrack& hot, float chi2, int prev_idx) { m_hots.push_back({hot, chi2, prev_idx}); @@ -412,10 +402,7 @@ namespace mkfit { int pickupLayer() const { return m_pickup_layer; } -#ifdef DUMPHITWINDOW - int seed_algo() const { return m_seed_algo; } - int seed_label() const { return m_seed_label; } -#endif + int seed_origin_index() const { return m_seed_origin_index; } private: trk_cand_vec_type m_trk_cands; @@ -425,11 +412,7 @@ namespace mkfit { short int m_lastHitIdx_before_bkwsearch = -1; short int m_nInsideMinusOneHits_before_bkwsearch = -1; short int m_nTailMinusOneHits_before_bkwsearch = -1; - -#ifdef DUMPHITWINDOW - int m_seed_algo = 0; - int m_seed_label = 0; -#endif + int m_seed_origin_index = -1; // seed index in the passed-in seed vector int m_hots_size = 0; std::vector m_hots; }; @@ -623,10 +606,10 @@ namespace mkfit { m_n_seeds_inserted -= n_removed; } - void insertSeed(const Track& seed, const track_score_func& score_func, int region, int pos) { + void insertSeed(const Track& seed, int seed_idx, const track_score_func& score_func, int region, int pos) { assert(pos < m_size); - m_candidates[pos].importSeed(seed, score_func, region); + m_candidates[pos].importSeed(seed, seed_idx, score_func, region); ++m_n_seeds_inserted; } diff --git a/RecoTracker/MkFitCore/interface/TrackerInfo.h b/RecoTracker/MkFitCore/interface/TrackerInfo.h index 940212283114e..965f47976b4f3 100644 --- a/RecoTracker/MkFitCore/interface/TrackerInfo.h +++ b/RecoTracker/MkFitCore/interface/TrackerInfo.h @@ -6,6 +6,7 @@ #include "RecoTracker/MkFitCore/interface/Config.h" #include #include +#include namespace mkfit { @@ -30,10 +31,12 @@ namespace mkfit { SVector3 pos; SVector3 zdir; SVector3 xdir; + float half_length; unsigned int detid; ModuleInfo() = default; - ModuleInfo(SVector3 p, SVector3 zd, SVector3 xd, unsigned int id) : pos(p), zdir(zd), xdir(xd), detid(id) {} + ModuleInfo(SVector3 p, SVector3 zd, SVector3 xd, float hl, unsigned int id) + : pos(p), zdir(zd), xdir(xd), half_length(hl), detid(id) {} }; //============================================================================== diff --git a/RecoTracker/MkFitCore/src/HitStructures.cc b/RecoTracker/MkFitCore/src/HitStructures.cc index b4201f6475757..da2afc2f60b0f 100644 --- a/RecoTracker/MkFitCore/src/HitStructures.cc +++ b/RecoTracker/MkFitCore/src/HitStructures.cc @@ -61,9 +61,22 @@ namespace mkfit { //============================================================================== + void LayerOfHits::reset() { + m_hit_infos.clear(); + m_ext_idcs.clear(); + m_min_ext_idx = std::numeric_limits::max(); + m_max_ext_idx = std::numeric_limits::min(); + m_n_hits = 0; + m_binnor.reset_contents(); + } + + //============================================================================== + void LayerOfHits::suckInHits(const HitVec &hitv) { - m_n_hits = hitv.size(); m_ext_hits = &hitv; + m_n_hits = hitv.size(); + + m_binnor.begin_registration(m_n_hits); #ifdef COPY_SORTED_HITS if (m_capacity < m_n_hits) { @@ -72,24 +85,31 @@ namespace mkfit { } #endif + std::vector hinfos; if (Config::usePhiQArrays) { - m_hit_phis.resize(m_n_hits); - m_hit_qs.resize(m_n_hits); - m_hit_infos.resize(m_n_hits); + hinfos.reserve(m_n_hits); + m_hit_infos.reserve(m_n_hits); } - m_binnor.reset_contents(); - m_binnor.begin_registration(m_n_hits); - for (unsigned int i = 0; i < m_n_hits; ++i) { const Hit &h = hitv[i]; - HitInfo hi = {h.phi(), m_is_barrel ? h.z() : h.r()}; + float phi = h.phi(); + float q = m_is_barrel ? h.z() : h.r(); - m_binnor.register_entry_safe(hi.phi, hi.q); + m_binnor.register_entry_safe(phi, q); if (Config::usePhiQArrays) { - m_hit_infos[i] = hi; + const float sqrt3 = std::sqrt(3); + float half_length, qbar; + if (m_is_barrel) { + half_length = sqrt3 * std::sqrt(h.ezz()); + qbar = h.r(); + } else { + half_length = sqrt3 * std::sqrt(h.exx() + h.eyy()); + qbar = h.z(); + } + hinfos.emplace_back(HitInfo({phi, q, half_length, qbar})); } } @@ -101,8 +121,7 @@ namespace mkfit { memcpy(&m_hits[i], &hitv[j], sizeof(Hit)); #endif if (Config::usePhiQArrays) { - m_hit_phis[i] = m_hit_infos[j].phi; - m_hit_qs[i] = m_hit_infos[j].q; + m_hit_infos.emplace_back(hinfos[j]); } } } @@ -131,14 +150,8 @@ namespace mkfit { void LayerOfHits::beginRegistrationOfHits(const HitVec &hitv) { m_ext_hits = &hitv; - m_n_hits = 0; - m_hit_infos.clear(); - m_ext_idcs.clear(); - m_min_ext_idx = std::numeric_limits::max(); - m_max_ext_idx = std::numeric_limits::min(); - m_binnor.reset_contents(); m_binnor.begin_registration(128); // initial reserve for cons vectors } @@ -149,12 +162,22 @@ namespace mkfit { m_min_ext_idx = std::min(m_min_ext_idx, idx); m_max_ext_idx = std::max(m_max_ext_idx, idx); - HitInfo hi = {h.phi(), m_is_barrel ? h.z() : h.r()}; + float phi = h.phi(); + float q = m_is_barrel ? h.z() : h.r(); - m_binnor.register_entry_safe(hi.phi, hi.q); + m_binnor.register_entry_safe(phi, q); if (Config::usePhiQArrays) { - m_hit_infos.emplace_back(hi); + const float sqrt3 = std::sqrt(3); + float half_length, qbar; + if (m_is_barrel) { + half_length = sqrt3 * std::sqrt(h.ezz()); + qbar = h.r(); + } else { + half_length = sqrt3 * std::sqrt(h.exx() + h.eyy()); + qbar = h.z(); + } + m_hit_infos.emplace_back(HitInfo({phi, q, half_length, qbar})); } } @@ -174,9 +197,10 @@ namespace mkfit { } #endif + std::vector hinfos; if (Config::usePhiQArrays) { - m_hit_phis.resize(m_n_hits); - m_hit_qs.resize(m_n_hits); + hinfos.swap(m_hit_infos); + m_hit_infos.reserve(m_n_hits); } for (unsigned int i = 0; i < m_n_hits; ++i) { @@ -188,8 +212,7 @@ namespace mkfit { #endif if (Config::usePhiQArrays) { - m_hit_phis[i] = m_hit_infos[j].phi; - m_hit_qs[i] = m_hit_infos[j].q; + m_hit_infos.emplace_back(hinfos[j]); } // Redirect m_binnor.m_ranks[i] to point to external/original index. diff --git a/RecoTracker/MkFitCore/src/Matriplex/Matriplex.h b/RecoTracker/MkFitCore/src/Matriplex/Matriplex.h index d7fea243db2ce..712afd240789b 100644 --- a/RecoTracker/MkFitCore/src/Matriplex/Matriplex.h +++ b/RecoTracker/MkFitCore/src/Matriplex/Matriplex.h @@ -8,7 +8,7 @@ namespace Matriplex { //------------------------------------------------------------------------------ template - class Matriplex { + class __attribute__((aligned(MPLEX_ALIGN))) Matriplex { public: typedef T value_type; @@ -21,7 +21,7 @@ namespace Matriplex { /// size of the whole matriplex static constexpr int kTotSize = N * kSize; - T fArray[kTotSize] __attribute__((aligned(64))); + T fArray[kTotSize]; Matriplex() {} Matriplex(T v) { setVal(v); } @@ -61,6 +61,124 @@ namespace Matriplex { return *this; } + Matriplex& operator=(T t) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] = t; + return *this; + } + + Matriplex& operator+=(T t) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] += t; + return *this; + } + + Matriplex& operator-=(T t) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] -= t; + return *this; + } + + Matriplex& operator*=(T t) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] *= t; + return *this; + } + + Matriplex& operator/=(T t) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] /= t; + return *this; + } + + Matriplex& operator+=(const Matriplex& a) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] += a.fArray[i]; + return *this; + } + + Matriplex& operator-=(const Matriplex& a) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] -= a.fArray[i]; + return *this; + } + + Matriplex& operator*=(const Matriplex& a) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] *= a.fArray[i]; + return *this; + } + + Matriplex& operator/=(const Matriplex& a) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] /= a.fArray[i]; + return *this; + } + + Matriplex& sqrt(const Matriplex& a) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] = std::sqrt(a.fArray[i]); + return *this; + } + Matriplex& sqrt() { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] = std::sqrt(fArray[i]); + return *this; + } + + Matriplex& sqr(const Matriplex& a) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] = a.fArray[i] * a.fArray[i]; + return *this; + } + Matriplex& sqr() { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] = fArray[i] * fArray[i]; + return *this; + } + + Matriplex& hypot(const Matriplex& a, const Matriplex& b) { + for (idx_t i = 0; i < kTotSize; ++i) { + fArray[i] = a.fArray[i] * a.fArray[i] + b.fArray[i] * b.fArray[i]; + } + return sqrt(); + } + + Matriplex& sin(const Matriplex& a) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] = std::sin(a.fArray[i]); + return *this; + } + Matriplex& sin() { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] = std::sin(fArray[i]); + return *this; + } + + Matriplex& cos(const Matriplex& a) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] = std::cos(a.fArray[i]); + return *this; + } + Matriplex& cos() { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] = std::cos(fArray[i]); + return *this; + } + + Matriplex& tan(const Matriplex& a) { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] = std::tan(a.fArray[i]); + return *this; + } + Matriplex& tan() { + for (idx_t i = 0; i < kTotSize; ++i) + fArray[i] = std::tan(fArray[i]); + return *this; + } + + //--------------------------------------------------------- + void copySlot(idx_t n, const Matriplex& m) { for (idx_t i = n; i < kTotSize; i += N) { fArray[i] = m.fArray[i]; @@ -182,11 +300,180 @@ namespace Matriplex { *(arr++) = fArray[i]; } } + + Matriplex ReduceFixedIJ(idx_t i, idx_t j) const { + Matriplex t; + for (idx_t n = 0; n < N; ++n) { + t[n] = constAt(n, i, j); + } + return t; + } }; template using MPlex = Matriplex; + //============================================================================== + // Operators + //============================================================================== + + template + MPlex operator+(const MPlex& a, const MPlex& b) { + MPlex t = a; + t += b; + return t; + } + + template + MPlex operator-(const MPlex& a, const MPlex& b) { + MPlex t = a; + t -= b; + return t; + } + + template + MPlex operator*(const MPlex& a, const MPlex& b) { + MPlex t = a; + t *= b; + return t; + } + + template + MPlex operator/(const MPlex& a, const MPlex& b) { + MPlex t = a; + t /= b; + return t; + } + + template + MPlex operator+(const MPlex& a, T b) { + MPlex t = a; + t += b; + return t; + } + + template + MPlex operator-(const MPlex& a, T b) { + MPlex t = a; + t -= b; + return t; + } + + template + MPlex operator*(const MPlex& a, T b) { + MPlex t = a; + t *= b; + return t; + } + + template + MPlex operator/(const MPlex& a, T b) { + MPlex t = a; + t /= b; + return t; + } + + template + MPlex operator+(T a, const MPlex& b) { + MPlex t = a; + t += b; + return t; + } + + template + MPlex operator-(T a, const MPlex& b) { + MPlex t = a; + t -= b; + return t; + } + + template + MPlex operator*(T a, const MPlex& b) { + MPlex t = a; + t *= b; + return t; + } + + template + MPlex operator/(T a, const MPlex& b) { + MPlex t = a; + t /= b; + return t; + } + + template + MPlex sqrt(const MPlex& a) { + MPlex t; + return t.sqrt(a); + } + + template + MPlex sqr(const MPlex& a) { + MPlex t; + return t.sqrt(a); + } + + template + MPlex hypot(const MPlex& a, const MPlex& b) { + MPlex t; + return t.hypot(a, b); + } + + template + MPlex sin(const MPlex& a) { + MPlex t; + return t.sin(a); + } + + template + MPlex cos(const MPlex& a) { + MPlex t; + return t.cos(a); + } + + template + void sincos(const MPlex& a, MPlex& s, MPlex& c) { + for (idx_t i = 0; i < a.kTotSize; ++i) { + s.fArray[i] = std::sin(a.fArray[i]); + c.fArray[i] = std::cos(a.fArray[i]); + } + } + + template + MPlex tan(const MPlex& a) { + MPlex t; + return t.tan(a); + } + + template + void min_max(const MPlex& a, + const MPlex& b, + MPlex& min, + MPlex& max) { + for (idx_t i = 0; i < a.kTotSize; ++i) { + min.fArray[i] = std::min(a.fArray[i], b.fArray[i]); + max.fArray[i] = std::max(a.fArray[i], b.fArray[i]); + } + } + + template + MPlex min(const MPlex& a, const MPlex& b) { + MPlex t; + for (idx_t i = 0; i < a.kTotSize; ++i) { + t.fArray[i] = std::min(a.fArray[i], b.fArray[i]); + } + return t; + } + + template + MPlex max(const MPlex& a, const MPlex& b) { + MPlex t; + for (idx_t i = 0; i < a.kTotSize; ++i) { + t.fArray[i] = std::max(a.fArray[i], b.fArray[i]); + } + return t; + } + //============================================================================== // Multiplications //============================================================================== diff --git a/RecoTracker/MkFitCore/src/Matriplex/MatriplexCommon.h b/RecoTracker/MkFitCore/src/Matriplex/MatriplexCommon.h index 653f5601d597d..d1a62ce4bdd07 100644 --- a/RecoTracker/MkFitCore/src/Matriplex/MatriplexCommon.h +++ b/RecoTracker/MkFitCore/src/Matriplex/MatriplexCommon.h @@ -16,6 +16,18 @@ #include #endif +#ifndef MPLEX_ALIGN +#if defined(__AVX512F__) +#define MPLEX_ALIGN 64 +#elif defined(__AVX__) || defined(__AVX2__) +#define MPLEX_ALIGN 32 +#elif defined(__SSE3__) +#define MPLEX_ALIGN 16 +#else +#define MPLEX_ALIGN 32 +#endif +#endif + #if defined(MPLEX_USE_INTRINSICS) // This seems unnecessary: __AVX__ is usually defined for all higher ISA extensions #if defined(__AVX__) || defined(__AVX512F__) diff --git a/RecoTracker/MkFitCore/src/Matriplex/MatriplexSym.h b/RecoTracker/MkFitCore/src/Matriplex/MatriplexSym.h index 6b8cf44222adc..151f616402440 100644 --- a/RecoTracker/MkFitCore/src/Matriplex/MatriplexSym.h +++ b/RecoTracker/MkFitCore/src/Matriplex/MatriplexSym.h @@ -22,7 +22,7 @@ namespace Matriplex { //------------------------------------------------------------------------------ template - class MatriplexSym { + class __attribute__((aligned(MPLEX_ALIGN))) MatriplexSym { public: typedef T value_type; @@ -35,7 +35,7 @@ namespace Matriplex { /// size of the whole matriplex static constexpr int kTotSize = N * kSize; - T fArray[kTotSize] __attribute__((aligned(64))); + T fArray[kTotSize]; MatriplexSym() {} MatriplexSym(T v) { setVal(v); } diff --git a/RecoTracker/MkFitCore/src/Matrix.h b/RecoTracker/MkFitCore/src/Matrix.h index ed5e6939b8692..809b438964e5e 100644 --- a/RecoTracker/MkFitCore/src/Matrix.h +++ b/RecoTracker/MkFitCore/src/Matrix.h @@ -65,6 +65,8 @@ namespace mkfit { typedef Matriplex::Matriplex MPlexQF; typedef Matriplex::Matriplex MPlexQI; typedef Matriplex::Matriplex MPlexQUI; + typedef Matriplex::Matriplex MPlexQH; + typedef Matriplex::Matriplex MPlexQUH; typedef Matriplex::Matriplex MPlexQB; diff --git a/RecoTracker/MkFitCore/src/MiniPropagators.cc b/RecoTracker/MkFitCore/src/MiniPropagators.cc new file mode 100644 index 0000000000000..78b99cab76c0f --- /dev/null +++ b/RecoTracker/MkFitCore/src/MiniPropagators.cc @@ -0,0 +1,241 @@ +#include "RecoTracker/MkFitCore/src/MiniPropagators.h" +#include "vdt/atan2.h" +#include "vdt/tan.h" +#include "vdt/sincos.h" +#include "vdt/sqrt.h" + +namespace mkfit::mini_propagators { + + State::State(const MPlexLV& par, int ti) { + x = par.constAt(ti, 0, 0); + y = par.constAt(ti, 1, 0); + z = par.constAt(ti, 2, 0); + const float pt = 1.0f / par.constAt(ti, 3, 0); + px = pt * std::cos(par.constAt(ti, 4, 0)); + py = pt * std::sin(par.constAt(ti, 4, 0)); + pz = pt / std::tan(par.constAt(ti, 5, 0)); + } + + bool InitialState::propagate_to_r(PropAlgo_e algo, float R, State& c, bool update_momentum) const { + switch (algo) { + case PA_Line: { + } + case PA_Quadratic: { + } + + case PA_Exact: { + // Momentum is always updated -- used as temporary for stepping. + const float k = 1.0f / inv_k; + + const float curv = 0.5f * inv_k * inv_pt; + const float oo_curv = 1.0f / curv; // 2 * radius of curvature + const float lambda = pz * inv_pt; + + float D = 0; + + c = *this; + c.dalpha = 0; + for (int i = 0; i < Config::Niter; ++i) { + // compute tangental and ideal distance for the current iteration. + // 3-rd order asin for symmetric incidence (shortest arc lenght). + float r0 = std::hypot(c.x, c.y); + float td = (R - r0) * curv; + float id = oo_curv * td * (1.0f + 0.16666666f * td * td); + // This would be for line approximation: + // float id = R - r0; + D += id; + + //printf("%-3d r0=%f R-r0=%f td=%f id=%f id_line=%f delta_id=%g\n", + // i, r0, R-r0, td, id, R - r0, id - (R-r0)); + + float alpha = id * inv_pt * inv_k; + float sina, cosa; + vdt::fast_sincosf(alpha, sina, cosa); + + // update parameters + c.dalpha += alpha; + c.x += k * (c.px * sina - c.py * (1.0f - cosa)); + c.y += k * (c.py * sina + c.px * (1.0f - cosa)); + + const float o_px = c.px; // copy before overwriting + c.px = c.px * cosa - c.py * sina; + c.py = c.py * cosa + o_px * sina; + } + + c.z += lambda * D; + } + } + // should have some epsilon constant / member? relative vs. abs? + c.fail_flag = std::abs(std::hypot(c.x, c.y) - R) < 0.1f ? 0 : 1; + return c.fail_flag; + } + + bool InitialState::propagate_to_z(PropAlgo_e algo, float Z, State& c, bool update_momentum) const { + switch (algo) { + case PA_Line: { + } + case PA_Quadratic: { + } + + case PA_Exact: { + const float k = 1.0f / inv_k; + + const float dz = Z - z; + const float alpha = dz * inv_k / pz; + + float sina, cosa; + vdt::fast_sincosf(alpha, sina, cosa); + + c.dalpha = alpha; + c.x = x + k * (px * sina - py * (1.0f - cosa)); + c.y = y + k * (py * sina + px * (1.0f - cosa)); + c.z = Z; + + if (update_momentum) { + c.px = px * cosa - py * sina; + c.py = py * cosa + px * sina; + c.pz = pz; + } + } break; + } + c.fail_flag = 0; + return c.fail_flag; + } + + //=========================================================================== + // Vectorized version + //=========================================================================== + + MPF fast_atan2(const MPF& y, const MPF& x) { + MPF t; + for (int i = 0; i < y.kTotSize; ++i) { + t[i] = vdt::fast_atan2f(y[i], x[i]); + } + return t; + } + + MPF fast_tan(const MPF& a) { + MPF t; + for (int i = 0; i < a.kTotSize; ++i) { + t[i] = vdt::fast_tanf(a[i]); + } + return t; + } + + void fast_sincos(const MPF& a, MPF& s, MPF& c) { + for (int i = 0; i < a.kTotSize; ++i) { + vdt::fast_sincosf(a[i], s[i], c[i]); + } + } + + StatePlex::StatePlex(const MPlexLV& par) { + x = par.ReduceFixedIJ(0, 0); + y = par.ReduceFixedIJ(1, 0); + z = par.ReduceFixedIJ(2, 0); + const MPF pt = 1.0f / par.ReduceFixedIJ(3, 0); + fast_sincos(par.ReduceFixedIJ(4, 0), py, px); + px *= pt; + py *= pt; + pz = pt / fast_tan(par.ReduceFixedIJ(5, 0)); + } + + // propagate to radius; returns number of failed propagations + int InitialStatePlex::propagate_to_r( + PropAlgo_e algo, const MPF& R, StatePlex& c, bool update_momentum, int N_proc) const { + switch (algo) { + case PA_Line: { + } + case PA_Quadratic: { + } + + case PA_Exact: { + // Momentum is always updated -- used as temporary for stepping. + const MPF k = 1.0f / inv_k; + + const MPF curv = 0.5f * inv_k * inv_pt; + const MPF oo_curv = 1.0f / curv; // 2 * radius of curvature + const MPF lambda = pz * inv_pt; + + MPF D = 0; + + c = *this; + c.dalpha = 0; + for (int i = 0; i < Config::Niter; ++i) { + // compute tangental and ideal distance for the current iteration. + // 3-rd order asin for symmetric incidence (shortest arc lenght). + MPF r0 = Matriplex::hypot(c.x, c.y); + MPF td = (R - r0) * curv; + MPF id = oo_curv * td * (1.0f + 0.16666666f * td * td); + // This would be for line approximation: + // float id = R - r0; + D += id; + + //printf("%-3d r0=%f R-r0=%f td=%f id=%f id_line=%f delta_id=%g\n", + // i, r0, R-r0, td, id, R - r0, id - (R-r0)); + + MPF alpha = id * inv_pt * inv_k; + + MPF sina, cosa; + fast_sincos(alpha, sina, cosa); + + // update parameters + c.dalpha += alpha; + c.x += k * (c.px * sina - c.py * (1.0f - cosa)); + c.y += k * (c.py * sina + c.px * (1.0f - cosa)); + + MPF o_px = c.px; // copy before overwriting + c.px = c.px * cosa - c.py * sina; + c.py = c.py * cosa + o_px * sina; + } + + c.z += lambda * D; + } + } + + // should have some epsilon constant / member? relative vs. abs? + MPF r = Matriplex::hypot(c.x, c.y); + c.fail_flag = 0; + int n_fail = 0; + for (int i = 0; i < N_proc; ++i) { + if (std::abs(R[i] - r[i]) > 0.1f) { + c.fail_flag[i] = 1; + ++n_fail; + } + } + return n_fail; + } + + int InitialStatePlex::propagate_to_z( + PropAlgo_e algo, const MPF& Z, StatePlex& c, bool update_momentum, int N_proc) const { + switch (algo) { + case PA_Line: { + } + case PA_Quadratic: { + } + + case PA_Exact: { + MPF k = 1.0f / inv_k; + + MPF dz = Z - z; + MPF alpha = dz * inv_k / pz; + + MPF sina, cosa; + fast_sincos(alpha, sina, cosa); + + c.dalpha = alpha; + c.x = x + k * (px * sina - py * (1.0f - cosa)); + c.y = y + k * (py * sina + px * (1.0f - cosa)); + c.z = Z; + + if (update_momentum) { + c.px = px * cosa - py * sina; + c.py = py * cosa + px * sina; + c.pz = pz; + } + } break; + } + c.fail_flag = 0; + return 0; + } + +} // namespace mkfit::mini_propagators diff --git a/RecoTracker/MkFitCore/src/MiniPropagators.h b/RecoTracker/MkFitCore/src/MiniPropagators.h new file mode 100644 index 0000000000000..6a8e07ac37fa5 --- /dev/null +++ b/RecoTracker/MkFitCore/src/MiniPropagators.h @@ -0,0 +1,79 @@ +#ifndef RecoTracker_MkFitCore_src_MiniPropagators_h +#define RecoTracker_MkFitCore_src_MiniPropagators_h + +#include "RecoTracker/MkFitCore/interface/Config.h" +#include "RecoTracker/MkFitCore/interface/MatrixSTypes.h" +#include "Matrix.h" + +namespace mkfit::mini_propagators { + + enum PropAlgo_e { PA_Line, PA_Quadratic, PA_Exact }; + + struct State { + float x, y, z; + float px, py, pz; + float dalpha; + int fail_flag; + + State() = default; + State(const MPlexLV& par, int ti); + }; + + struct InitialState : public State { + float inv_pt, inv_k; + float theta; + + InitialState(const MPlexLV& par, const MPlexQI& chg, int ti) + : InitialState(State(par, ti), chg.constAt(ti, 0, 0), par.constAt(ti, 3, 0), par.constAt(ti, 5, 0)) {} + + InitialState(State s, short charge, float ipt, float tht, float bf = Config::Bfield) + : State(s), inv_pt(ipt), theta(tht) { + inv_k = ((charge < 0) ? 0.01f : -0.01f) * Const::sol * bf; + } + + bool propagate_to_r(PropAlgo_e algo, float R, State& c, bool update_momentum) const; + bool propagate_to_z(PropAlgo_e algo, float Z, State& c, bool update_momentum) const; + }; + + //----------------------------------------------------------- + // Vectorized version + //----------------------------------------------------------- + + using MPF = MPlexQF; + using MPI = MPlexQI; + + MPF fast_atan2(const MPF& y, const MPF& x); + MPF fast_tan(const MPF& a); + void fast_sincos(const MPF& a, MPF& s, MPF& c); + + struct StatePlex { + MPF x, y, z; + MPF px, py, pz; + MPF dalpha; + MPI fail_flag{0}; + + StatePlex() = default; + StatePlex(const MPlexLV& par); + }; + + struct InitialStatePlex : public StatePlex { + MPF inv_pt, inv_k; + MPF theta; + + InitialStatePlex(const MPlexLV& par, const MPI& chg) + : InitialStatePlex(StatePlex(par), chg, par.ReduceFixedIJ(3, 0), par.ReduceFixedIJ(5, 0)) {} + + InitialStatePlex(StatePlex s, MPI charge, MPF ipt, MPF tht, float bf = Config::Bfield) + : StatePlex(s), inv_pt(ipt), theta(tht) { + for (int i = 0; i < inv_k.kTotSize; ++i) { + inv_k[i] = ((charge[i] < 0) ? 0.01f : -0.01f) * Const::sol * bf; + } + } + + int propagate_to_r(PropAlgo_e algo, const MPF& R, StatePlex& c, bool update_momentum, int N_proc = NN) const; + int propagate_to_z(PropAlgo_e algo, const MPF& Z, StatePlex& c, bool update_momentum, int N_proc = NN) const; + }; + +}; // namespace mkfit::mini_propagators + +#endif diff --git a/RecoTracker/MkFitCore/src/MkBase.h b/RecoTracker/MkFitCore/src/MkBase.h index e7ddfd91faf35..2244ea1f8cccb 100644 --- a/RecoTracker/MkFitCore/src/MkBase.h +++ b/RecoTracker/MkFitCore/src/MkBase.h @@ -21,6 +21,7 @@ namespace mkfit { float getPar(int itrack, int i, int par) const { return m_Par[i].constAt(itrack, par, 0); } float radiusSqr(int itrack, int i) const { return hipo_sqr(getPar(itrack, i, 0), getPar(itrack, i, 1)); } + float radius(int itrack, int i) const { return hipo(getPar(itrack, i, 0), getPar(itrack, i, 1)); } //---------------------------------------------------------------------------- diff --git a/RecoTracker/MkFitCore/src/MkBuilder.cc b/RecoTracker/MkFitCore/src/MkBuilder.cc index 7df60a442ee03..97c0c6dd06f01 100644 --- a/RecoTracker/MkFitCore/src/MkBuilder.cc +++ b/RecoTracker/MkFitCore/src/MkBuilder.cc @@ -268,7 +268,7 @@ namespace mkfit { int j = seeds_sorted ? i : ranks[i]; int reg = part.m_region[j]; const Track &seed = in_seeds[j]; - insert_seed(seed, reg, seed_cursors[reg]++); + insert_seed(seed, j, reg, seed_cursors[reg]++); HitOnTrack hot = seed.getLastHitOnTrack(); m_seedMinLastLayer[reg] = std::min(m_seedMinLastLayer[reg], hot.layer); @@ -441,7 +441,7 @@ namespace mkfit { assert(!in_seeds.empty()); m_tracks.resize(in_seeds.size()); - import_seeds(in_seeds, seeds_sorted, [&](const Track &seed, int region, int pos) { + import_seeds(in_seeds, seeds_sorted, [&](const Track &seed, int seed_idx, int region, int pos) { m_tracks[pos] = seed; m_tracks[pos].setNSeedHits(seed.nTotalHits()); m_tracks[pos].setEtaRegion(region); @@ -515,10 +515,13 @@ namespace mkfit { prev_layer = curr_layer; curr_layer = layer_plan_it.layer(); mkfndr->setup(prop_config, + m_job->m_iter_config, m_job->m_iter_config.m_params, m_job->m_iter_config.m_layer_configs[curr_layer], st_par, m_job->get_mask_for_layer(curr_layer), + m_event, + region, m_job->m_in_fwd); const LayerOfHits &layer_of_hits = m_job->m_event_of_hits[curr_layer]; @@ -610,8 +613,8 @@ namespace mkfit { m_event_of_comb_cands.reset((int)in_seeds.size(), m_job->max_max_cands()); - import_seeds(in_seeds, seeds_sorted, [&](const Track &seed, int region, int pos) { - m_event_of_comb_cands.insertSeed(seed, m_job->steering_params(region).m_track_scorer, region, pos); + import_seeds(in_seeds, seeds_sorted, [&](const Track &seed, int seed_idx, int region, int pos) { + m_event_of_comb_cands.insertSeed(seed, seed_idx, m_job->steering_params(region).m_track_scorer, region, pos); }); } @@ -699,10 +702,7 @@ namespace mkfit { TrackCand &cand = m_event_of_comb_cands[seed_cand_idx[ti].first][seed_cand_idx[ti].second]; WSR_Result &w = mkfndr->m_XWsrResult[ti - itrack]; - // XXXX-4 Low pT tracks can miss a barrel layer ... and should be stopped - const float cand_r = - std::hypot(mkfndr->getPar(ti - itrack, MkBase::iP, 0), mkfndr->getPar(ti - itrack, MkBase::iP, 1)); - + // Low pT tracks can miss a barrel layer ... and should be stopped dprintf("WSR Check label %d, seed %d, cand %d score %f -> wsr %d, in_gap %d\n", cand.label(), seed_cand_idx[ti].first, @@ -711,20 +711,25 @@ namespace mkfit { w.m_wsr, w.m_in_gap); - if (layer_info.is_barrel() && cand_r < layer_info.rin()) { - // Fake outside so it does not get processed in FindTracks Std/CE... and - // create a stopped replica in barrel and original copy if there is - // still chance to hit endcaps. - dprintf("Barrel cand propagated to r=%f ... layer is %f - %f\n", cand_r, layer_info.rin(), layer_info.rout()); - - mkfndr->m_XHitSize[ti - itrack] = 0; + if (w.m_wsr == WSR_Failed) { + // Fake outside so it does not get processed in FindTracks BH/Std/CE. + // [ Should add handling of WSR_Failed there, perhaps. ] w.m_wsr = WSR_Outside; - tmp_cands[seed_cand_idx[ti].first - start_seed].push_back(cand); - if (region == TrackerInfo::Reg_Barrel) { - dprintf(" creating extra stopped held back candidate\n"); - tmp_cands[seed_cand_idx[ti].first - start_seed].back().addHitIdx(-2, layer_info.layer_id(), 0); + if (layer_info.is_barrel()) { + dprintf("Barrel cand propagation failed, got to r=%f ... layer is %f - %f\n", + mkfndr->radius(ti - itrack, MkBase::iP), + layer_info.rin(), + layer_info.rout()); + // In barrel region, create a stopped replica. In transition region keep the original copy + // as there is still a chance to hit endcaps. + tmp_cands[seed_cand_idx[ti].first - start_seed].push_back(cand); + if (region == TrackerInfo::Reg_Barrel) { + dprintf(" creating extra stopped held back candidate\n"); + tmp_cands[seed_cand_idx[ti].first - start_seed].back().addHitIdx(-2, layer_info.layer_id(), 0); + } } + // Never happens for endcap / propToZ } else if (w.m_wsr == WSR_Outside) { dprintf(" creating extra held back candidate\n"); tmp_cands[seed_cand_idx[ti].first - start_seed].push_back(cand); @@ -801,18 +806,22 @@ namespace mkfit { prev_layer = curr_layer; curr_layer = layer_plan_it.layer(); mkfndr->setup(prop_config, + m_job->m_iter_config, iter_params, m_job->m_iter_config.m_layer_configs[curr_layer], st_par, m_job->get_mask_for_layer(curr_layer), + m_event, + region, m_job->m_in_fwd); - dprintf("\n* Processing layer %d\n", curr_layer); - const LayerOfHits &layer_of_hits = m_job->m_event_of_hits[curr_layer]; const LayerInfo &layer_info = trk_info.layer(curr_layer); const FindingFoos &fnd_foos = FindingFoos::get_finding_foos(layer_info.is_barrel()); + dprintf("\n* Processing layer %d\n", curr_layer); + mkfndr->begin_layer(layer_of_hits); + int theEndCand = find_tracks_unroll_candidates(seed_cand_idx, start_seed, end_seed, @@ -902,7 +911,7 @@ namespace mkfit { tmp_cands[is].clear(); } } - + mkfndr->end_layer(); } // end of layer loop mkfndr->release(); @@ -1010,20 +1019,24 @@ namespace mkfit { prev_layer = curr_layer; curr_layer = layer_plan_it.layer(); mkfndr->setup(prop_config, + m_job->m_iter_config, iter_params, m_job->m_iter_config.m_layer_configs[curr_layer], st_par, m_job->get_mask_for_layer(curr_layer), + m_event, + region, m_job->m_in_fwd); const bool pickup_only = layer_plan_it.is_pickup_only(); - dprintf("\n\n* Processing layer %d, %s\n\n", curr_layer, pickup_only ? "pickup only" : "full finding"); - const LayerInfo &layer_info = trk_info.layer(curr_layer); const LayerOfHits &layer_of_hits = m_job->m_event_of_hits[curr_layer]; const FindingFoos &fnd_foos = FindingFoos::get_finding_foos(layer_info.is_barrel()); + dprintf("\n\n* Processing layer %d, %s\n\n", curr_layer, pickup_only ? "pickup only" : "full finding"); + mkfndr->begin_layer(layer_of_hits); + const int theEndCand = find_tracks_unroll_candidates( seed_cand_idx, start_seed, end_seed, curr_layer, prev_layer, pickup_only, iteration_dir); @@ -1074,10 +1087,6 @@ namespace mkfit { dprint("now get hit range"); -#ifdef DUMPHITWINDOW - mkfndr->m_event = m_event; -#endif - mkfndr->selectHitIndices(layer_of_hits, end - itrack); find_tracks_handle_missed_layers( @@ -1134,7 +1143,7 @@ namespace mkfit { } } #endif - + mkfndr->end_layer(); } // end of layer loop cloner.end_eta_bin(); @@ -1184,6 +1193,7 @@ namespace mkfit { void MkBuilder::fit_cands_BH(MkFinder *mkfndr, int start_cand, int end_cand, int region) { const SteeringParams &st_par = m_job->steering_params(region); const PropagationConfig &prop_config = m_job->m_trk_info.prop_config(); + mkfndr->setup_bkfit(prop_config, st_par, m_event); #ifdef DEBUG_FINAL_FIT EventOfCombCandidates &eoccs = m_event_of_comb_cands; bool debug = true; @@ -1280,7 +1290,7 @@ namespace mkfit { EventOfCombCandidates &eoccs = m_event_of_comb_cands; const SteeringParams &st_par = m_job->steering_params(region); const PropagationConfig &prop_config = m_job->m_trk_info.prop_config(); - mkfndr->setup_bkfit(prop_config, st_par); + mkfndr->setup_bkfit(prop_config, st_par, m_event); int step = NN; for (int icand = start_cand; icand < end_cand; icand += step) { @@ -1304,7 +1314,6 @@ namespace mkfit { "sx_t/F:sy_t/F:sz_t/F:d_xy/F:d_z/F\n"); first = false; } - mkfndr->m_event = m_event; #endif // input tracks diff --git a/RecoTracker/MkFitCore/src/MkFinder.cc b/RecoTracker/MkFitCore/src/MkFinder.cc index 7347a31a04525..c6eeb64f0b119 100644 --- a/RecoTracker/MkFitCore/src/MkFinder.cc +++ b/RecoTracker/MkFitCore/src/MkFinder.cc @@ -6,46 +6,92 @@ #include "FindingFoos.h" #include "KalmanUtilsMPlex.h" #include "MatriplexPackers.h" +#include "MiniPropagators.h" //#define DEBUG #include "Debug.h" -#if (defined(DUMPHITWINDOW) || defined(DEBUG_BACKWARD_FIT)) && defined(MKFIT_STANDALONE) +#if defined(MKFIT_STANDALONE) #include "RecoTracker/MkFitCore/standalone/Event.h" #endif +#ifdef RNT_DUMP_MkF_SelHitIdcs +// declares struct RntIfc_selectHitIndices rnt_shi in unnamed namespace; +#include "RecoTracker/MkFitCore/standalone/RntDumper/MkFinder_selectHitIndices.icc" +#endif + +#include "vdt/atan2.h" + #include +#include namespace mkfit { void MkFinder::setup(const PropagationConfig &pc, + const IterationConfig &ic, const IterationParams &ip, const IterationLayerConfig &ilc, const SteeringParams &sp, const std::vector *ihm, + const Event *ev, + int region, bool infwd) { m_prop_config = &pc; + m_iteration_config = ⁣ m_iteration_params = &ip; m_iteration_layer_config = &ilc; m_steering_params = &sp; m_iteration_hit_mask = ihm; + m_event = ev; + m_current_region = region; m_in_fwd = infwd; } - void MkFinder::setup_bkfit(const PropagationConfig &pc, const SteeringParams &sp) { + void MkFinder::setup_bkfit(const PropagationConfig &pc, const SteeringParams &sp, const Event *ev) { m_prop_config = &pc; m_steering_params = &sp; + m_event = ev; } void MkFinder::release() { m_prop_config = nullptr; + m_iteration_config = nullptr; m_iteration_params = nullptr; m_iteration_layer_config = nullptr; m_steering_params = nullptr; m_iteration_hit_mask = nullptr; + m_event = nullptr; + m_current_region = -1; m_in_fwd = true; } + void MkFinder::begin_layer(const LayerOfHits &layer_of_hits) { +#ifdef RNT_DUMP_MkF_SelHitIdcs + const LayerOfHits &L = layer_of_hits; + const LayerInfo &LI = *L.layer_info(); + rnt_shi.ResetH(); + rnt_shi.ResetF(); + *rnt_shi.h = {m_event->evtID(), + m_iteration_config->m_iteration_index, + m_iteration_config->m_track_algorithm, + m_current_region, + L.layer_id(), + L.is_barrel() ? LI.rin() : LI.zmin(), + LI.is_barrel() ? LI.rout() : LI.zmax(), + L.is_barrel(), + L.is_pixel(), + L.is_stereo()}; + *rnt_shi.f = *rnt_shi.h; +#endif + } + + void MkFinder::end_layer() { +#ifdef RNT_DUMP_MkF_SelHitIdcs + rnt_shi.FillH(); + rnt_shi.FillF(); +#endif + } + //============================================================================== // Input / Output TracksAndHitIdx //============================================================================== @@ -94,13 +140,9 @@ namespace mkfit { copy_in(trk, imp, iI); -#ifdef DUMPHITWINDOW - m_SeedAlgo(imp, 0, 0) = tracks[idxs[i].first].seed_algo(); - m_SeedLabel(imp, 0, 0) = tracks[idxs[i].first].seed_label(); -#endif - m_SeedIdx(imp, 0, 0) = idxs[i].first; m_CandIdx(imp, 0, 0) = idxs[i].second; + m_SeedOriginIdx[imp] = tracks[idxs[i].first].seed_origin_index(); } } @@ -122,13 +164,9 @@ namespace mkfit { copy_in(trk, imp, iI); -#ifdef DUMPHITWINDOW - m_SeedAlgo(imp, 0, 0) = tracks[idxs[i].seed_idx].seed_algo(); - m_SeedLabel(imp, 0, 0) = tracks[idxs[i].seed_idx.seed_label(); -#endif - m_SeedIdx(imp, 0, 0) = idxs[i].seed_idx; m_CandIdx(imp, 0, 0) = idxs[i].cand_idx; + m_SeedOriginIdx[imp] = tracks[idxs[i].seed_idx].seed_origin_index(); const Hit &hit = layer_of_hits.refHit(idxs[i].hit_idx); m_msErr.copyIn(imp, hit.errArray()); @@ -153,13 +191,9 @@ namespace mkfit { copy_in(trk, imp, iI); -#ifdef DUMPHITWINDOW - m_SeedAlgo(imp, 0, 0) = tracks[idxs[i].first].seed_algo(); - m_SeedLabel(imp, 0, 0) = tracks[idxs[i].first].seed_label(); -#endif - m_SeedIdx(imp, 0, 0) = idxs[i].first; m_CandIdx(imp, 0, 0) = idxs[i].second.trkIdx; + m_SeedOriginIdx[imp] = tracks[idxs[i].first].seed_origin_index(); } } @@ -245,7 +279,7 @@ namespace mkfit { // SelectHitIndices //============================================================================== - void MkFinder::selectHitIndices(const LayerOfHits &layer_of_hits, const int N_proc) { + void MkFinder::selectHitIndices(const LayerOfHits &layer_of_hits, const int N_proc, bool fill_binsearch_only) { // bool debug = true; using bidx_t = LayerOfHits::bin_index_t; using bcnt_t = LayerOfHits::bin_content_t; @@ -391,19 +425,38 @@ namespace mkfit { } } +#ifdef RNT_DUMP_MkF_SelHitIdcs + if (fill_binsearch_only) { + // XXX loop over good indices (prepared in V2) and put in V1 BinSearch results + for (auto i : rnt_shi.f_h_idcs) { + CandInfo &ci = (*rnt_shi.ci)[rnt_shi.f_h_remap[i]]; + ci.bso = BinSearch({phiv[i], + dphiv[i], + qv[i], + dqv[i], + pb1v[i], + pb2v[i], + qb1v[i], + qb2v[i], + m_XWsrResult[i].m_wsr, + m_XWsrResult[i].m_in_gap, + false}); + } + return; + } +#endif + // Vectorizing this makes it run slower! //#pragma omp simd for (int itrack = 0; itrack < N_proc; ++itrack) { // PROP-FAIL-ENABLE The following to be enabled when propagation failure // detection is properly implemented in propagate-to-R/Z. - // if (m_FailFlag[itrack]) { - // m_XWsrResult[itrack].m_wsr = WSR_Failed; - // m_XHitSize[itrack] = -1; - // continue; - // } + if (m_FailFlag[itrack]) { + m_XWsrResult[itrack].m_wsr = WSR_Failed; + continue; + } if (m_XWsrResult[itrack].m_wsr == WSR_Outside) { - m_XHitSize[itrack] = -1; continue; } @@ -413,44 +466,23 @@ namespace mkfit { const bidx_t pb1 = pb1v[itrack]; const bidx_t pb2 = pb2v[itrack]; - // Used only by usePhiQArrays const float q = qv[itrack]; const float phi = phiv[itrack]; const float dphi = dphiv[itrack]; const float dq = dqv[itrack]; + // clang-format off dprintf(" %2d/%2d: %6.3f %6.3f %6.6f %7.5f %3u %3u %4u %4u\n", L.layer_id(), itrack, q, phi, dq, dphi, qb1, qb2, pb1, pb2); // clang-format on - // MT: One could iterate in "spiral" order, to pick hits close to the center. - // http://stackoverflow.com/questions/398299/looping-in-a-spiral - // This would then work best with relatively small bin sizes. - // Or, set them up so I can always take 3x3 array around the intersection. #if defined(DUMPHITWINDOW) && defined(MKFIT_STANDALONE) - int thisseedmcid = -999999; - { - int seedlabel = m_SeedLabel(itrack, 0, 0); - TrackVec &seedtracks = m_event->seedTracks_; - int thisidx = -999999; - for (int i = 0; i < int(seedtracks.size()); ++i) { - auto &thisseed = seedtracks[i]; - if (thisseed.label() == seedlabel) { - thisidx = i; - break; - } - } - if (thisidx > -999999) { - auto &seedtrack = m_event->seedTracks_[thisidx]; - std::vector thismcHitIDs; - seedtrack.mcHitIDsVec(m_event->layerHits_, m_event->simHitsInfo_, thismcHitIDs); - if (std::adjacent_find(thismcHitIDs.begin(), thismcHitIDs.end(), std::not_equal_to<>()) == - thismcHitIDs.end()) { - thisseedmcid = thismcHitIDs.at(0); - } - } - } + const auto ngr = [](float f) { return isFinite(f) ? f : -999.0f; }; + + const int seed_lbl = m_event->currentSeed(m_SeedOriginIdx[itrack]).label(); + Event::SimLabelFromHits slfh = m_event->simLabelForCurrentSeed(m_SeedOriginIdx[itrack]); + const int seed_mcid = (slfh.is_set() && slfh.good_frac() > 0.7f) ? slfh.label : -999999; #endif for (bidx_t qi = qb1; qi != qb2; ++qi) { @@ -474,7 +506,7 @@ namespace mkfit { for (bcnt_t hi = pbi.begin(); hi < pbi.end(); ++hi) { // MT: Access into m_hit_zs and m_hit_phis is 1% run-time each. - unsigned int hi_orig = L.getOriginalHitIndex(hi); + const unsigned int hi_orig = L.getOriginalHitIndex(hi); if (m_iteration_hit_mask && (*m_iteration_hit_mask)[hi_orig]) { dprintf( @@ -489,7 +521,15 @@ namespace mkfit { const float ddq = std::abs(q - L.hit_q(hi)); const float ddphi = cdist(std::abs(phi - L.hit_phi(hi))); + // clang-format off + dprintf(" SHI %3u %4u %5u %6.3f %6.3f %6.4f %7.5f %s\n", + qi, pi, hi, L.hit_q(hi), L.hit_phi(hi), + ddq, ddphi, (ddq < dq && ddphi < dphi) ? "PASS" : "FAIL"); + // clang-format on + #if defined(DUMPHITWINDOW) && defined(MKFIT_STANDALONE) + // clang-format off + MPlexQF thisOutChi2; { const MCHitInfo &mchinfo = m_event->simHitsInfo_[L.refHit(hi).mcHitID()]; int mchid = mchinfo.mcTrackID(); @@ -517,11 +557,11 @@ namespace mkfit { st_z = simtrack.z(); } - const Hit &thishit = L.refHit(hi); + const Hit &thishit = L.refHit(hi_orig); m_msErr.copyIn(itrack, thishit.errArray()); m_msPar.copyIn(itrack, thishit.posArray()); - MPlexQF thisOutChi2; + MPlexQI propFail; MPlexLV tmpPropPar; const FindingFoos &fnd_foos = FindingFoos::get_finding_foos(L.is_barrel()); (*fnd_foos.m_compute_chi2_foo)(m_Err[iI], @@ -531,6 +571,7 @@ namespace mkfit { m_msPar, thisOutChi2, tmpPropPar, + propFail, N_proc, m_prop_config->finding_intra_layer_pflags, m_prop_config->finding_requires_propagation_to_hit_pos); @@ -540,53 +581,29 @@ namespace mkfit { float hz = thishit.z(); float hr = std::hypot(hx, hy); float hphi = std::atan2(hy, hx); - float hex = std::sqrt(thishit.exx()); - if (std::isnan(hex)) - hex = -999.; - float hey = std::sqrt(thishit.eyy()); - if (std::isnan(hey)) - hey = -999.; - float hez = std::sqrt(thishit.ezz()); - if (std::isnan(hez)) - hez = -999.; - float her = std::sqrt( + float hex = ngr( std::sqrt(thishit.exx()) ); + float hey = ngr( std::sqrt(thishit.eyy()) ); + float hez = ngr( std::sqrt(thishit.ezz()) ); + float her = ngr( std::sqrt( (hx * hx * thishit.exx() + hy * hy * thishit.eyy() + 2.0f * hx * hy * m_msErr.At(itrack, 0, 1)) / - (hr * hr)); - if (std::isnan(her)) - her = -999.; - float hephi = std::sqrt(thishit.ephi()); - if (std::isnan(hephi)) - hephi = -999.; - float hchi2 = thisOutChi2[itrack]; - if (std::isnan(hchi2)) - hchi2 = -999.; + (hr * hr)) ); + float hephi = ngr( std::sqrt(thishit.ephi()) ); + float hchi2 = ngr( thisOutChi2[itrack] ); float tx = m_Par[iI].At(itrack, 0, 0); float ty = m_Par[iI].At(itrack, 1, 0); float tz = m_Par[iI].At(itrack, 2, 0); float tr = std::hypot(tx, ty); float tphi = std::atan2(ty, tx); - float tchi2 = m_Chi2(itrack, 0, 0); - if (std::isnan(tchi2)) - tchi2 = -999.; - float tex = std::sqrt(m_Err[iI].At(itrack, 0, 0)); - if (std::isnan(tex)) - tex = -999.; - float tey = std::sqrt(m_Err[iI].At(itrack, 1, 1)); - if (std::isnan(tey)) - tey = -999.; - float tez = std::sqrt(m_Err[iI].At(itrack, 2, 2)); - if (std::isnan(tez)) - tez = -999.; - float ter = std::sqrt( + // float tchi2 = ngr( m_Chi2(itrack, 0, 0) ); // unused + float tex = ngr( std::sqrt(m_Err[iI].At(itrack, 0, 0)) ); + float tey = ngr( std::sqrt(m_Err[iI].At(itrack, 1, 1)) ); + float tez = ngr( std::sqrt(m_Err[iI].At(itrack, 2, 2)) ); + float ter = ngr( std::sqrt( (tx * tx * tex * tex + ty * ty * tey * tey + 2.0f * tx * ty * m_Err[iI].At(itrack, 0, 1)) / - (tr * tr)); - if (std::isnan(ter)) - ter = -999.; - float tephi = std::sqrt( + (tr * tr)) ); + float tephi = ngr( std::sqrt( (ty * ty * tex * tex + tx * tx * tey * tey - 2.0f * tx * ty * m_Err[iI].At(itrack, 0, 1)) / - (tr * tr * tr * tr)); - if (std::isnan(tephi)) - tephi = -999.; + (tr * tr * tr * tr)) ); float ht_dxy = std::hypot(hx - tx, hy - ty); float ht_dz = hz - tz; float ht_dphi = cdist(std::abs(hphi - tphi)); @@ -595,12 +612,12 @@ namespace mkfit { if (first) { printf( "HITWINDOWSEL " - "evt_id/I:" + "evt_id/I:track_algo/I:" "lyr_id/I:lyr_isbrl/I:hit_idx/I:" "trk_cnt/I:trk_idx/I:trk_label/I:" "trk_pt/F:trk_eta/F:trk_mphi/F:trk_chi2/F:" "nhits/I:" - "seed_idx/I:seed_label/I:seed_algo/I:seed_mcid/I:" + "seed_idx/I:seed_label/I:seed_mcid/I:" "hit_mcid/I:" "st_isfindable/I:st_prodtype/I:st_label/I:" "st_pt/F:st_eta/F:st_phi/F:" @@ -618,14 +635,13 @@ namespace mkfit { if (!(std::isnan(phi)) && !(std::isnan(getEta(m_Par[iI].At(itrack, 5, 0))))) { //|| std::isnan(ter) || std::isnan(her) || std::isnan(m_Chi2(itrack, 0, 0)) || std::isnan(hchi2))) - // clang-format off printf("HITWINDOWSEL " - "%d " + "%d %d" "%d %d %d " "%d %d %d " "%6.3f %6.3f %6.3f %6.3f " "%d " - "%d %d %d %d " + "%d %d %d " "%d " "%d %d %d " "%6.3f %6.3f %6.3f " @@ -638,12 +654,12 @@ namespace mkfit { "%6.3f %6.3f %6.3f " "%6.3f" "\n", - m_event->evtID(), - L.layer_id(), L.is_barrel(), L.getOriginalHitIndex(hi), + m_event->evtID(), m_iteration_config->m_track_algorithm, + L.layer_id(), L.is_barrel(), hi_orig, itrack, m_CandIdx(itrack, 0, 0), m_Label(itrack, 0, 0), 1.0f / m_Par[iI].At(itrack, 3, 0), getEta(m_Par[iI].At(itrack, 5, 0)), m_Par[iI].At(itrack, 4, 0), m_Chi2(itrack, 0, 0), m_NFoundHits(itrack, 0, 0), - m_SeedIdx(itrack, 0, 0), m_SeedLabel(itrack, 0, 0), m_SeedAlgo(itrack, 0, 0), thisseedmcid, + m_SeedOriginIdx(itrack, 0, 0), seed_lbl, seed_mcid, mchid, st_isfindable, st_prodtype, st_label, st_pt, st_eta, st_phi, @@ -655,20 +671,15 @@ namespace mkfit { hex, hey, her, hephi, hez, ht_dxy, ht_dz, ht_dphi, hchi2); - // clang-format on } } + // clang-format on #endif if (ddq >= dq) continue; if (ddphi >= dphi) continue; - // clang-format off - dprintf(" SHI %3u %4u %5u %6.3f %6.3f %6.4f %7.5f %s\n", - qi, pi, hi, L.hit_q(hi), L.hit_phi(hi), - ddq, ddphi, (ddq < dq && ddphi < dphi) ? "PASS" : "FAIL"); - // clang-format on // MT: Removing extra check gives full efficiency ... // and means our error estimations are wrong! @@ -694,6 +705,280 @@ namespace mkfit { } //itrack } + //============================================================================== + // SelectHitIndicesV2 + //============================================================================== + + void MkFinder::selectHitIndicesV2(const LayerOfHits &layer_of_hits, const int N_proc) { + // bool debug = true; + using bidx_t = LayerOfHits::bin_index_t; + using bcnt_t = LayerOfHits::bin_content_t; + const LayerOfHits &L = layer_of_hits; + const LayerInfo &LI = *L.layer_info(); + + const int iI = iP; + + dprintf("LayerOfHits::SelectHitIndicesV2 %s layer=%d N_proc=%d\n", + L.is_barrel() ? "barrel" : "endcap", + L.layer_id(), + N_proc); + +#ifdef RNT_DUMP_MkF_SelHitIdcs + rnt_shi.InnerIdcsReset(N_proc); + for (int i = 0; i < N_proc; ++i) { + auto slfh = m_event->simLabelForCurrentSeed(m_SeedOriginIdx[i]); + if (m_FailFlag[i]) { + rnt_shi.RegisterFailedProp(i, m_Par[1 - iI], m_Par[iI], m_event, m_SeedOriginIdx[i]); + } else if (slfh.is_set()) { + rnt_shi.RegisterGoodProp(i, m_Par[iI], m_event, m_SeedOriginIdx[i]); + // get BinSearch result from V1. + selectHitIndices(layer_of_hits, N_proc, true); + } // else ... could do something about the bad seeds ... probably better to collect elsewhere. + } +#endif + + constexpr int NEW_MAX_HIT = 6; // 4 - 6 give about the same # of tracks in quality-val + constexpr float DDPHI_PRESEL_FAC = 2.0f; + constexpr float DDQ_PRESEL_FAC = 1.2f; + constexpr float PHI_BIN_EXTRA_FAC = 2.75f; + constexpr float Q_BIN_EXTRA_FAC = 1.6f; + + namespace mp = mini_propagators; + struct Bins { + MPlexQUH q0, q1, q2, p1, p2; + mp::InitialStatePlex isp; + mp::StatePlex sp1, sp2; + int n_proc; + + // debug & ntuple dump -- to be local in functions + MPlexQF phi_c, dphi; + MPlexQF q_c, qmin, qmax; + + Bins(const MPlexLV &par, const MPlexQI &chg, int np = NN) : isp(par, chg), n_proc(np) {} + + void prop_to_limits(const LayerInfo &li) { + // Positions 1 and 2 should really be by "propagation order", 1 is the closest/ + // This should also work for backward propagation so not exactly trivial. + // Also, do not really need propagation to center. + if (li.is_barrel()) { + isp.propagate_to_r(mp::PA_Exact, li.rin(), sp1, true, n_proc); + isp.propagate_to_r(mp::PA_Exact, li.rout(), sp2, true, n_proc); + } else { + isp.propagate_to_z(mp::PA_Exact, li.zmin(), sp1, true, n_proc); + isp.propagate_to_z(mp::PA_Exact, li.zmax(), sp2, true, n_proc); + } + } + + void find_bin_ranges(const LayerInfo &li, const LayerOfHits &loh) { + // Below made members for debugging + // MPlexQF phi_c, dphi_min, dphi_max; + phi_c = mp::fast_atan2(isp.y, isp.x); + + // Matriplex::min_max(sp1.dphi, sp2.dphi, dphi_min, dphi_max); + // the above is wrong: dalpha is not dphi --> renamed variable in State + MPlexQF xp1, xp2, pmin, pmax; + xp1 = mp::fast_atan2(sp1.y, sp1.x); + xp2 = mp::fast_atan2(sp2.y, sp2.x); + Matriplex::min_max(xp1, xp2, pmin, pmax); + // Matriplex::min_max(mp::fast_atan2(sp1.y, sp1.x), smp::fast_atan2(sp2.y, sp2.x), pmin, pmax); + MPlexQF dp = pmax - pmin; + phi_c = 0.5f * (pmax + pmin); + for (int ii = 0; ii < n_proc; ++ii) { + if (dp[ii] > Const::PI) { + std::swap(pmax[ii], pmin[ii]); + dp[ii] = Const::TwoPI - dp[ii]; + phi_c[ii] = Const::PI - phi_c[ii]; + } + dphi[ii] = 0.5f * dp[ii]; + // printf("phic: %f p1: %f p2: %f pmin: %f pmax: %f dphi: %f\n", + // phi_c[ii], xp1[ii], xp2[ii], pmin[ii], pmax[ii], dphi[ii]); + } + + // MPlexQF qmin, qmax; + if (li.is_barrel()) { + Matriplex::min_max(sp1.z, sp2.z, qmin, qmax); + q_c = isp.z; + } else { + Matriplex::min_max(Matriplex::hypot(sp1.x, sp1.y), Matriplex::hypot(sp2.x, sp2.y), qmin, qmax); + q_c = Matriplex::hypot(isp.x, isp.y); + } + + for (int i = 0; i < p1.kTotSize; ++i) { + // Clamp crazy sizes. This actually only happens when prop-fail flag is set. + // const float dphi_clamp = 0.1; + // if (dphi_min[i] > 0.0f || dphi_min[i] < -dphi_clamp) dphi_min[i] = -dphi_clamp; + // if (dphi_max[i] < 0.0f || dphi_max[i] > dphi_clampf) dphi_max[i] = dphi_clamp; + p1[i] = loh.phiBinChecked(pmin[i] - PHI_BIN_EXTRA_FAC * 0.0123f); + p2[i] = loh.phiBinChecked(pmax[i] + PHI_BIN_EXTRA_FAC * 0.0123f); + + q0[i] = loh.qBinChecked(q_c[i]); + q1[i] = loh.qBinChecked(qmin[i] - Q_BIN_EXTRA_FAC * 0.5f * li.q_bin()); + q2[i] = loh.qBinChecked(qmax[i] + Q_BIN_EXTRA_FAC * 0.5f * li.q_bin()) + 1; + } + } + }; + + Bins B(m_Par[iI], m_Chg, N_proc); + B.prop_to_limits(LI); + B.find_bin_ranges(LI, L); + + for (int i = 0; i < N_proc; ++i) { + m_XHitSize[i] = 0; + // Notify failure. Ideally should be detected before selectHitIndices(). + if (m_FailFlag[i]) { + m_XWsrResult[i].m_wsr = WSR_Failed; + } else { + if (LI.is_barrel()) { + m_XWsrResult[i] = L.is_within_z_sensitive_region(B.q_c[i], 0.5f * (B.q2[i] - B.q1[i])); + } else { + m_XWsrResult[i] = L.is_within_r_sensitive_region(B.q_c[i], 0.5f * (B.q2[i] - B.q1[i])); + } + } + } + + // for (int i = 0; i < N_proc; ++i) { + // printf("BinCheck %c %+8.6f %+8.6f | %3d %3d - %3d %3d || | %2d %2d - %2d %2d\n", LI.is_barrel() ? 'B' : 'E', + // B.phi[i], B.dphi[i], B.p1[i], B.p2[i], pb1v[i], pb2v[i], + // B.q[i], B.dq[i], B.q1[i], B.q2[i], qb1v[i], qb2v[i]); + // } + +#ifdef RNT_DUMP_MkF_SelHitIdcs + for (auto i : rnt_shi.f_h_idcs) { + CandInfo &ci = (*rnt_shi.ci)[rnt_shi.f_h_remap[i]]; + ci.bsn = BinSearch({B.phi_c[i], + B.dphi[i], + B.q_c[i], + 0.5f * (B.q2[i] - B.q1[i]), + B.p1[i], + B.p2[i], + B.q1[i], + B.q2[i], + m_XWsrResult[i].m_wsr, + m_XWsrResult[i].m_in_gap, + false}); + ci.ps_min = statep2propstate(B.sp1, i); + ci.ps_max = statep2propstate(B.sp2, i); + } +#endif + + struct PQE { + float score; + unsigned int hit_index; + }; + auto pqe_cmp = [](const PQE &a, const PQE &b) { return a.score < b.score; }; + std::priority_queue, decltype(pqe_cmp)> pqueue(pqe_cmp); + int pqueue_size = 0; + + // Vectorizing this makes it run slower! + //#pragma omp simd + for (int itrack = 0; itrack < N_proc; ++itrack) { + if (m_FailFlag[itrack]) { + m_XWsrResult[itrack].m_wsr = WSR_Failed; + continue; + } + + if (m_XWsrResult[itrack].m_wsr == WSR_Outside) { + continue; + } + + // New binning -- known to be too restrictive, so scaled up. Probably esp. in stereo layers. + // Also, could take track covariance dphi / dq extras + known tilt stuff. + const bidx_t qb = B.q0[itrack]; + const bidx_t qb1 = B.q1[itrack]; + const bidx_t qb2 = B.q2[itrack]; + const bidx_t pb1 = B.p1[itrack]; + const bidx_t pb2 = B.p2[itrack]; + + // clang-format off + dprintf(" %2d/%2d: %6.3f %6.3f %6.6f %7.5f %3u %3u %4u %4u\n", + L.layer_id(), itrack, qv[itrack], phi[itrack], dqv[itrack], dphiv[itrack], + qb1, qb2, pb1, pb2); + // clang-format on + + mp::InitialState mp_is(m_Par[iI], m_Chg, itrack); + mp::State mp_s; + + for (bidx_t qi = qb1; qi != qb2; ++qi) { + for (bidx_t pi = pb1; pi != pb2; pi = L.phiMaskApply(pi + 1)) { + // Limit to central Q-bin + if (qi == qb && L.isBinDead(pi, qi) == true) { + dprint("dead module for track in layer=" << L.layer_id() << " qb=" << qi << " pi=" << pi << " q=" << q + << " phi=" << phi); + m_XWsrResult[itrack].m_in_gap = true; + } + + // It might make sense to make first loop to extract bin indices + // and issue prefetches at the same time. + // Then enter vectorized loop to actually collect the hits in proper order. + + //SK: ~20x1024 bin sizes give mostly 1 hit per bin. Commented out for 128 bins or less + // #pragma nounroll + auto pbi = L.phiQBinContent(pi, qi); + for (bcnt_t hi = pbi.begin(); hi < pbi.end(); ++hi) { + // MT: Access into m_hit_zs and m_hit_phis is 1% run-time each. + + const unsigned int hi_orig = L.getOriginalHitIndex(hi); + + if (m_iteration_hit_mask && (*m_iteration_hit_mask)[hi_orig]) { + dprintf( + "Yay, denying masked hit on layer %u, hi %u, orig idx %u\n", L.layer_info()->layer_id(), hi, hi_orig); + continue; + } + + if (m_XHitSize[itrack] >= MPlexHitIdxMax) + break; + + float new_q, new_phi, new_ddphi, new_ddq; + bool prop_fail; + + if (L.is_barrel()) { + prop_fail = mp_is.propagate_to_r(mp::PA_Exact, L.hit_qbar(hi), mp_s, true); + new_q = mp_s.z; + } else { + prop_fail = mp_is.propagate_to_z(mp::PA_Exact, L.hit_qbar(hi), mp_s, true); + new_q = std::hypot(mp_s.x, mp_s.y); + } + + new_phi = vdt::fast_atan2f(mp_s.y, mp_s.x); + new_ddphi = cdist(std::abs(new_phi - L.hit_phi(hi))); + new_ddq = std::abs(new_q - L.hit_q(hi)); + + bool dqdphi_presel = + new_ddq < DDQ_PRESEL_FAC * L.hit_q_half_length(hi) && new_ddphi < DDPHI_PRESEL_FAC * 0.0123f; + + // clang-format off + dprintf(" SHI %3u %4u %5u %6.3f %6.3f %6.4f %7.5f PROP-%s %s\n", + qi, pi, hi, L.hit_q(hi), L.hit_phi(hi), + ddq, ddphi, prop_fail ? "FAIL" : "OK", dqdphi_presel ? "PASS" : "REJECT"); + // clang-format on + + if (prop_fail || !dqdphi_presel) + continue; + if (pqueue_size < NEW_MAX_HIT) { + pqueue.push({new_ddphi, hi_orig}); + ++pqueue_size; + } else if (new_ddphi < pqueue.top().score) { + pqueue.pop(); + pqueue.push({new_ddphi, hi_orig}); + } + } //hi + } //pi + } //qi + + dprintf(" PQUEUE (%d)", pqueue_size); + // Reverse hits so best dphis/scores come first in the hit-index list. + m_XHitSize[itrack] = pqueue_size; + while (pqueue_size) { + --pqueue_size; + m_XHitArr.At(itrack, pqueue_size, 0) = pqueue.top().hit_index; + dprintf(" %d: %f %d", pqueue_size, pqueue.top().score, pqueue.top().hit_index); + pqueue.pop(); + } + dprintf("\n"); + + } //itrack + } + //============================================================================== // AddBestHit - Best Hit Track Finding //============================================================================== diff --git a/RecoTracker/MkFitCore/src/MkFinder.h b/RecoTracker/MkFitCore/src/MkFinder.h index e4514dc10d7fa..17da11ad8aa9d 100644 --- a/RecoTracker/MkFitCore/src/MkFinder.h +++ b/RecoTracker/MkFitCore/src/MkFinder.h @@ -24,9 +24,7 @@ namespace mkfit { class IterationLayerConfig; class SteeringParams; -#if defined(DUMPHITWINDOW) or defined(DEBUG_BACKWARD_FIT) class Event; -#endif struct UpdateIndices { int seed_idx; @@ -50,14 +48,20 @@ namespace mkfit { MkFinder() {} void setup(const PropagationConfig &pc, + const IterationConfig &ic, const IterationParams &ip, const IterationLayerConfig &ilc, const SteeringParams &sp, const std::vector *ihm, + const Event *ev, + int region, bool infwd); - void setup_bkfit(const PropagationConfig &pc, const SteeringParams &sp); + void setup_bkfit(const PropagationConfig &pc, const SteeringParams &sp, const Event *ev); void release(); + void begin_layer(const LayerOfHits &layer_of_hits); + void end_layer(); + //---------------------------------------------------------------------------- void inputTracksAndHitIdx(const std::vector &tracks, int beg, int end, bool inputProp); @@ -116,7 +120,8 @@ namespace mkfit { float getHitSelDynamicChi2Cut(const int itrk, const int ipar); - void selectHitIndices(const LayerOfHits &layer_of_hits, const int N_proc); + void selectHitIndices(const LayerOfHits &layer_of_hits, const int N_proc, bool fill_binsearch_only = false); + void selectHitIndicesV2(const LayerOfHits &layer_of_hits, const int N_proc); void addBestHit(const LayerOfHits &layer_of_hits, const int N_proc, const FindingFoos &fnd_foos); @@ -285,14 +290,9 @@ namespace mkfit { HitOnTrack m_HoTArrs[NN][Config::nMaxTrkHits]; -#if defined(DUMPHITWINDOW) or defined(DEBUG_BACKWARD_FIT) - MPlexQI m_SeedAlgo; // seed algorithm - MPlexQI m_SeedLabel; // seed label - Event *m_event; -#endif - - MPlexQI m_SeedIdx; // seed index in local thread (for bookkeeping at thread level) - MPlexQI m_CandIdx; // candidate index for the given seed (for bookkeeping of clone engine) + MPlexQI m_SeedIdx; // seed index in local thread (for bookkeeping at thread level) + MPlexQI m_CandIdx; // candidate index for the given seed (for bookkeeping of clone engine) + MPlexQI m_SeedOriginIdx; // seed index in MkBuilder seed input vector MPlexQI m_Stopped; // Flag for BestHit that a track has been stopped (and copied out already) @@ -328,10 +328,13 @@ namespace mkfit { // MPlexLV candParAtCurrHit; const PropagationConfig *m_prop_config = nullptr; + const IterationConfig *m_iteration_config = nullptr; const IterationParams *m_iteration_params = nullptr; const IterationLayerConfig *m_iteration_layer_config = nullptr; const SteeringParams *m_steering_params = nullptr; const std::vector *m_iteration_hit_mask = nullptr; + const Event *m_event = nullptr; + int m_current_region = -1; bool m_in_fwd = true; // Backward fit diff --git a/RecoTracker/MkFitCore/src/Track.cc b/RecoTracker/MkFitCore/src/Track.cc index 24ebb7d530417..4fdd34437b01f 100644 --- a/RecoTracker/MkFitCore/src/Track.cc +++ b/RecoTracker/MkFitCore/src/Track.cc @@ -200,77 +200,7 @@ namespace mkfit { return std::hypot(x_center - x_bs, y_center - y_bs) - abs_ooc_half; } } - - const char* TrackBase::algoint_to_cstr(int algo) { - static const char* const names[] = {"undefAlgorithm", - "ctf", - "duplicateMerge", - "cosmics", - "initialStep", - "lowPtTripletStep", - "pixelPairStep", - "detachedTripletStep", - "mixedTripletStep", - "pixelLessStep", - "tobTecStep", - "jetCoreRegionalStep", - "conversionStep", - "muonSeededStepInOut", - "muonSeededStepOutIn", - "outInEcalSeededConv", - "inOutEcalSeededConv", - "nuclInter", - "standAloneMuon", - "globalMuon", - "cosmicStandAloneMuon", - "cosmicGlobalMuon", - "highPtTripletStep", - "lowPtQuadStep", - "detachedQuadStep", - "reservedForUpgrades1", - "reservedForUpgrades2", - "bTagGhostTracks", - "beamhalo", - "gsf", - "hltPixel", - "hltIter0", - "hltIter1", - "hltIter2", - "hltIter3", - "hltIter4", - "hltIterX", - "hiRegitMuInitialStep", - "hiRegitMuLowPtTripletStep", - "hiRegitMuPixelPairStep", - "hiRegitMuDetachedTripletStep", - "hiRegitMuMixedTripletStep", - "hiRegitMuPixelLessStep", - "hiRegitMuTobTecStep", - "hiRegitMuMuonSeededStepInOut", - "hiRegitMuMuonSeededStepOutIn", - "algoSize"}; - - if (algo < 0 || algo >= (int)TrackAlgorithm::algoSize) - return names[0]; - return names[algo]; - } - - //============================================================================== - // Track - //============================================================================== - - void Track::resizeHitsForInput() { - bzero(&hitsOnTrk_, sizeof(hitsOnTrk_)); - hitsOnTrk_.resize(lastHitIdx_ + 1); - } - - void Track::sortHitsByLayer() { - std::stable_sort(&hitsOnTrk_[0], &hitsOnTrk_[lastHitIdx_ + 1], [](const auto& h1, const auto& h2) { - return h1.layer < h2.layer; - }); - } - - float Track::swimPhiToR(const float x0, const float y0) const { + float TrackBase::swimPhiToR(const float x0, const float y0) const { const float dR = getHypot(x() - x0, y() - y0); // XXX-ASSUMPTION-ERROR can not always reach R, should see what callers expect. // For now return PI to signal apex on the ohter side of the helix. @@ -280,13 +210,13 @@ namespace mkfit { return squashPhiGeneral(momPhi() - dPhi); } - bool Track::canReachRadius(float R) const { + bool TrackBase::canReachRadius(float R) const { const float k = ((charge() < 0) ? 100.0f : -100.0f) / (Const::sol * Config::Bfield); const float ooc = 2.0f * k * pT(); return std::abs(ooc) > R - std::hypot(x(), y()); } - float Track::maxReachRadius() const { + float TrackBase::maxReachRadius() const { const float k = ((charge() < 0) ? 100.0f : -100.0f) / (Const::sol * Config::Bfield); const float abs_ooc_half = std::abs(k * pT()); // center of helix in x,y plane @@ -295,7 +225,7 @@ namespace mkfit { return std::hypot(x_center, y_center) + abs_ooc_half; } - float Track::zAtR(float R, float* r_reached) const { + float TrackBase::zAtR(float R, float* r_reached) const { float xc = x(); float yc = y(); float pxc = px(); @@ -372,7 +302,7 @@ namespace mkfit { // ---------------------------------------------------------------- } - float Track::rAtZ(float Z) const { + float TrackBase::rAtZ(float Z) const { float xc = x(); float yc = y(); float pxc = px(); @@ -398,6 +328,75 @@ namespace mkfit { return std::hypot(xc, yc); } + const char* TrackBase::algoint_to_cstr(int algo) { + static const char* const names[] = {"undefAlgorithm", + "ctf", + "duplicateMerge", + "cosmics", + "initialStep", + "lowPtTripletStep", + "pixelPairStep", + "detachedTripletStep", + "mixedTripletStep", + "pixelLessStep", + "tobTecStep", + "jetCoreRegionalStep", + "conversionStep", + "muonSeededStepInOut", + "muonSeededStepOutIn", + "outInEcalSeededConv", + "inOutEcalSeededConv", + "nuclInter", + "standAloneMuon", + "globalMuon", + "cosmicStandAloneMuon", + "cosmicGlobalMuon", + "highPtTripletStep", + "lowPtQuadStep", + "detachedQuadStep", + "reservedForUpgrades1", + "reservedForUpgrades2", + "bTagGhostTracks", + "beamhalo", + "gsf", + "hltPixel", + "hltIter0", + "hltIter1", + "hltIter2", + "hltIter3", + "hltIter4", + "hltIterX", + "hiRegitMuInitialStep", + "hiRegitMuLowPtTripletStep", + "hiRegitMuPixelPairStep", + "hiRegitMuDetachedTripletStep", + "hiRegitMuMixedTripletStep", + "hiRegitMuPixelLessStep", + "hiRegitMuTobTecStep", + "hiRegitMuMuonSeededStepInOut", + "hiRegitMuMuonSeededStepOutIn", + "algoSize"}; + + if (algo < 0 || algo >= (int)TrackAlgorithm::algoSize) + return names[0]; + return names[algo]; + } + + //============================================================================== + // Track + //============================================================================== + + void Track::resizeHitsForInput() { + bzero(&hitsOnTrk_, sizeof(hitsOnTrk_)); + hitsOnTrk_.resize(lastHitIdx_ + 1); + } + + void Track::sortHitsByLayer() { + std::stable_sort(&hitsOnTrk_[0], &hitsOnTrk_[lastHitIdx_ + 1], [](const auto& h1, const auto& h2) { + return h1.layer < h2.layer; + }); + } + //============================================================================== void print(const TrackState& s) { diff --git a/RecoTracker/MkFitCore/src/TrackStructures.cc b/RecoTracker/MkFitCore/src/TrackStructures.cc index c5b7ffd03fc08..e26d21130979a 100644 --- a/RecoTracker/MkFitCore/src/TrackStructures.cc +++ b/RecoTracker/MkFitCore/src/TrackStructures.cc @@ -50,15 +50,12 @@ namespace mkfit { // CombCandidate //============================================================================== - void CombCandidate::importSeed(const Track &seed, const track_score_func &score_func, int region) { + void CombCandidate::importSeed(const Track &seed, int seed_idx, const track_score_func &score_func, int region) { m_trk_cands.emplace_back(TrackCand(seed, this)); m_state = CombCandidate::Dormant; m_pickup_layer = seed.getLastHitLyr(); -#ifdef DUMPHITWINDOW - m_seed_algo = seed.algoint(); - m_seed_label = seed.label(); -#endif + m_seed_origin_index = seed_idx; TrackCand &cand = m_trk_cands.back(); cand.setNSeedHits(seed.nTotalHits()); diff --git a/RecoTracker/MkFitCore/standalone/ConfigStandalone.h b/RecoTracker/MkFitCore/standalone/ConfigStandalone.h index 1d49a8bb8c897..82e79998568fc 100644 --- a/RecoTracker/MkFitCore/standalone/ConfigStandalone.h +++ b/RecoTracker/MkFitCore/standalone/ConfigStandalone.h @@ -81,13 +81,6 @@ namespace mkfit { constexpr float varZ = Config::hitposerrZ * Config::hitposerrZ; constexpr float varR = Config::hitposerrR * Config::hitposerrR; - // scattering simulation - constexpr float X0 = - 9.370; // cm, from http://pdg.lbl.gov/2014/AtomicNuclearProperties/HTML/silicon_Si.html // Pb = 0.5612 cm - constexpr float xr = - 0.1; // -assumes radial impact. This is bigger than what we have in main --> shouldn't it be the parameter below??? if radial impact?? - //const float xr = std::sqrt(Config::beamspotX*Config::beamspotX + Config::beamspotY*Config::beamspotY); - // Config for seeding constexpr int nlayers_per_seed_max = 4; // Needed for allocation of arrays on stack. constexpr float chi2seedcut = 9.0; @@ -115,7 +108,7 @@ namespace mkfit { 0.0078; // 0.0075; // errors used for MC only fit, straight from sim tracks, outward with simple geometry constexpr float phierr049 = 0.0017; // 0.0017; constexpr float thetaerr049 = 0.0033; // 0.0031; - // parameters for layers 0,1,2 // --> ENDTOEND with "real seeding", fit is outward by definition, with poly geo + // parameters for layers 0,1,2 --> with "real seeding", fit is outward by definition, with poly geo constexpr float ptinverr012 = 0.12007; // 0.1789; -->old values from only MC seeds constexpr float phierr012 = 1.0; // found empirically 0.00646; // 0.0071 constexpr float thetaerr012 = 0.2; // also found empirically 0.01366; // 0.0130; diff --git a/RecoTracker/MkFitCore/standalone/Event.cc b/RecoTracker/MkFitCore/standalone/Event.cc index 474e4ccc388dd..825a15f5e0d7a 100644 --- a/RecoTracker/MkFitCore/standalone/Event.cc +++ b/RecoTracker/MkFitCore/standalone/Event.cc @@ -333,8 +333,10 @@ namespace mkfit { } #endif #ifdef DUMP_LAYER_HITS + // clang-format off printf("Read %i layers\n", nl); int total_hits = 0; + float min_delta = 100, max_delta = -100; for (int il = 0; il < nl; il++) { if (layerHits_[il].empty()) continue; @@ -343,18 +345,29 @@ namespace mkfit { total_hits += layerHits_[il].size(); for (int ih = 0; ih < (int)layerHits_[il].size(); ih++) { const Hit &hit = layerHits_[il][ih]; - printf(" mcHitID=%5d r=%10g x=%10g y=%10g z=%10g sx=%10.4g sy=%10.4e sz=%10.4e\n", - hit.mcHitID(), - hit.r(), - hit.x(), - hit.y(), - hit.z(), - std::sqrt(hit.exx()), - std::sqrt(hit.eyy()), - std::sqrt(hit.ezz())); + const LayerInfo &linfo = Config::TrkInfo[il]; + unsigned int mid = hit.detIDinLayer(); + const ModuleInfo &mi = linfo.module_info(mid); + + if ( ! linfo.is_pixel() && ! linfo.is_barrel()) { + float delta = mi.half_length - std::sqrt(3)*std::sqrt(hit.exx() + hit.eyy()); + min_delta = std::min(min_delta, delta); + max_delta = std::max(max_delta, delta); + } + // continue; + + printf(" mcHitID=%5d r=%10g x=%10g y=%10g z=%10g" + " sx=%10.4g sy=%10.4e sz=%10.4e sxy=%10.4e, mhl=%10.4e, delta=%10.4e\n", + hit.mcHitID(), hit.r(), hit.x(), hit.y(), hit.z(), + std::sqrt(hit.exx()), std::sqrt(hit.eyy()), std::sqrt(hit.ezz()), + std::sqrt(hit.exx() + hit.eyy()), + mi.half_length, + mi.half_length - std::sqrt(3)*std::sqrt(hit.exx() + hit.eyy())); } } - printf("Total hits in all layers = %d\n", total_hits); + printf("Total hits in all layers = %d; endcap strips: min_delta=%.5f max_delta=%.5f\n", + total_hits, min_delta, max_delta); + // clang-format on #endif #ifdef DUMP_REC_TRACKS printf("Read %i rectracks\n", nert); @@ -845,6 +858,60 @@ namespace mkfit { } } + //============================================================================== + // Handling of current seed vectors and MC label reconstruction from hit data + //============================================================================== + + void Event::setCurrentSeedTracks(const TrackVec &seeds) { currentSeedTracks_ = &seeds; } + const Track &Event::currentSeed(int i) const { return (*currentSeedTracks_)[i]; } + + Event::SimLabelFromHits Event::simLabelForCurrentSeed(int i) const { + assert(currentSeedTracks_ != nullptr); + + if (currentSeedSimFromHits_.empty()) { + currentSeedSimFromHits_.resize(currentSeedTracks_->size()); + + for (int si = 0; si < (int)currentSeedTracks_->size(); ++si) { + const Track &s = currentSeed(si); + // printf("%3d (%d): [", si, s.label()); + std::map lab_cnt; + for (int hi = 0; hi < s.nTotalHits(); ++hi) { + auto hot = s.getHitOnTrack(hi); + // printf(" %d", hot.index); + if (hot.index < 0) + continue; + const Hit &h = layerHits_[hot.layer][hot.index]; + int hl = simHitsInfo_[h.mcHitID()].mcTrackID_; + // printf(" (%d)", hl); + if (hl >= 0) + ++lab_cnt[hl]; + } + int max_c = -1, max_l = -1; + for (auto &x : lab_cnt) { + if (x.second > max_c) { + max_l = x.first; + max_c = x.second; + } else if (x.second == max_c) { + max_l = -1; + } + } + if (max_c < 0) { + max_c = 0; + max_l = -1; + } + // printf(" ] -> %d %d => %d\n", s.nTotalHits(), max_c, max_l); + currentSeedSimFromHits_[si] = {s.nTotalHits(), max_c, max_l}; + } + } + + return currentSeedSimFromHits_[i]; + } + + void Event::resetCurrentSeedTracks() { + currentSeedTracks_ = nullptr; + currentSeedSimFromHits_.clear(); + } + //============================================================================== // DataFile //============================================================================== diff --git a/RecoTracker/MkFitCore/standalone/Event.h b/RecoTracker/MkFitCore/standalone/Event.h index f8ec0970100c0..18419be8d64ad 100644 --- a/RecoTracker/MkFitCore/standalone/Event.h +++ b/RecoTracker/MkFitCore/standalone/Event.h @@ -49,8 +49,21 @@ namespace mkfit { Validation &validation_; + // For seed access in deep data dumpers. + struct SimLabelFromHits { + int n_hits = 0, n_match = 0, label = -1; + float good_frac() const { return (float)n_match / n_hits; } + bool is_set() const { return label >= 0; } + }; + void setCurrentSeedTracks(const TrackVec &seeds); + const Track ¤tSeed(int i) const; + SimLabelFromHits simLabelForCurrentSeed(int i) const; + void resetCurrentSeedTracks(); + private: int evtID_; + const TrackVec *currentSeedTracks_ = nullptr; + mutable std::vector currentSeedSimFromHits_; public: BeamSpot beamSpot_; // XXXX Read/Write of BeamSpot + file-version bump or extra-section to be added. diff --git a/RecoTracker/MkFitCore/standalone/Geoms/CylCowWLids.cc b/RecoTracker/MkFitCore/standalone/Geoms/CylCowWLids.cc index 735057b755c77..0f03b103f687a 100644 --- a/RecoTracker/MkFitCore/standalone/Geoms/CylCowWLids.cc +++ b/RecoTracker/MkFitCore/standalone/Geoms/CylCowWLids.cc @@ -108,7 +108,7 @@ namespace { // Actual coverage for tracks with z = 3cm is 2.4 float full_eta = 2.5; float full_eta_pix_0 = 2.55; // To account for BS z-spread - float full_eta_ec_in[] = {0, 2.525, 2.515}; + float full_eta_ec_in[] = {0, 2.525, 2.515, 2.505}; float pix_0 = 4, pix_sep = 6; float pix_z0 = 24, pix_zgrow = 6; diff --git a/RecoTracker/MkFitCore/standalone/Makefile b/RecoTracker/MkFitCore/standalone/Makefile index 40e2ce6f91083..84dfebe0136b2 100644 --- a/RecoTracker/MkFitCore/standalone/Makefile +++ b/RecoTracker/MkFitCore/standalone/Makefile @@ -10,7 +10,6 @@ TGTS := ${LIB_CORE} .PHONY: all clean distclean -all: ${TGTS} SRCS := $(wildcard ${CORE_DIR}/src/*.cc) \ $(wildcard ${CORE_DIR}/src/Matriplex/*.cc) \ @@ -23,9 +22,35 @@ vpath %.cc ${CORE_DIR}/src ${CORE_DIR}/src/Matriplex ${SADIR} AUTO_TGTS := +# Begin ROOT dumpers +# XXX can remove, is var empty if undefined ???? +RNTD_OBJS := + +ifdef WITH_ROOT + +LIB_RNTDUMP := ../libMicRntDump.so +RNTDUMP_DICT_PCM := ../MicRntDumpDict_rdict.pcm + +RNTD_SRCS := $(wildcard ${SADIR}/RntDumper/*.cc) +RNTD_SRCB := $(notdir ${RNTD_SRCS}) +RNTD_OBJS := $(RNTD_SRCB:.cc=.o) +DEPS += $(RNTD_SRCB:.cc=.d) +TGTS += ${LIB_RNTDUMP} + +#LIB_CORE_DEPS := ${LIB_RNTDUMP} +#LIB_CORE_LD_EXTRA := -L.. -lMicRntDump + +vpath %.cc ${SADIR}/RntDumper + +endif +# End ROOT dumpers + +all: ${TGTS} + # Begin Matriplex auto-matriplex: +### Matriplex auto-generated fragments are included in source. ### ${MAKE} -f ${CORE_DIR}/src/Matriplex auto && touch $@ touch $@ @@ -41,7 +66,7 @@ endif clean-local: -rm -f ${TGTS} *.d *.o *.om *.so - -rm -rf main.dSYM + -rm -rf MicRntDumpDict.cc auto-matriplex main.dSYM -rm -rf plotting/*.so plotting/*.d plotting/*.pcm clean: clean-local @@ -52,13 +77,25 @@ distclean: clean-local -rm -f ${LIB_CORE} ### cd Matriplex && ${MAKE} distclean -${LIB_CORE}: ${OBJS} +${LIB_CORE}: ${OBJS} ${LIB_RNTDUMP} @mkdir -p $(@D) - ${CXX} ${CXXFLAGS} ${VEC_HOST} ${OBJS} -shared -o $@ ${LDFLAGS_HOST} ${LDFLAGS} + ${CXX} ${CXXFLAGS} ${VEC_HOST} ${OBJS} -shared -o $@ ${LDFLAGS_HOST} ${LDFLAGS} ${LIB_CORE_LD_EXTRA} + +ifdef WITH_ROOT +${LIB_RNTDUMP}: ${RNTD_OBJS} MicRntDumpDict.o + ${CXX} ${CXXFLAGS} ${VEC_HOST} ${RNTD_OBJS} MicRntDumpDict.o -shared -o $@ ${LDFLAGS_HOST} ${LDFLAGS} -lROOTEve -${OBJS}: %.o: %.cc %.d +MicRntDumpDict.cc ${RNTDUMP_DICT_PCM} &: ${SADIR}/RntDumper/RntStructs.h ${SADIR}/RntDumper/RntStructs_Linkdef.h + rootcling -I. -I${SRCDIR} -f MicRntDumpDict.cc ${SADIR}/RntDumper/RntStructs.h ${SADIR}/RntDumper/RntStructs_Linkdef.h + mv MicRntDumpDict_rdict.pcm ${RNTDUMP_DICT_PCM} +endif + +${OBJS} ${RNTD_OBJS}: %.o: %.cc %.d ${CXX} ${CPPFLAGS} ${CXXFLAGS} ${VEC_HOST} -c -o $@ $< +%.s: %.cc + ${CXX} ${CPPFLAGS} ${CXXFLAGS} ${VEC_HOST} -S -o $@ $< + %.d: %.cc ${MAKEDEPEND} -o $@ $< @@ -67,6 +104,7 @@ echo: @echo SRCS = ${SRCS} @echo DEPS = ${DEPS} @echo OBJS = ${OBJS} + @echo RNTD_OBJS = ${RNTD_OBJS} echo_cc_defs: ${CXX} -dM -E -mavx2 - < /dev/null diff --git a/RecoTracker/MkFitCore/standalone/Makefile.config b/RecoTracker/MkFitCore/standalone/Makefile.config index de0fc4afd5486..67398920fdd6b 100644 --- a/RecoTracker/MkFitCore/standalone/Makefile.config +++ b/RecoTracker/MkFitCore/standalone/Makefile.config @@ -47,7 +47,8 @@ endif # 3. Optimization # -O3 implies vectorization and simd (but not AVX) -OPT := -g -O3 +# But use -Ofast rather than -O3 so gcc can auto-vectorize sin/cos with libmvec (glibc >= 2.22) +OPT := -g -Ofast # 4. Vectorization settings ifdef AVX_512 @@ -71,35 +72,12 @@ USE_INTRINSICS := -DMPLEX_USE_INTRINSICS #USE_INTRINSICS := -DMPT_SIZE=1 USE_VTUNE_NOTIFY := 1 -# 6. MIC stuff - obsolete - -# 7. OSX hack (is there a good way to detect clang?) -# MT needs this on OSX-10.8, c++ -v -# Apple LLVM version 5.1 (clang-503.0.40) (based on LLVM 3.4svn) -# OSX_CXXFLAGS := -stdlib=libc++ -# And with gcc-4.8.1 from cms we need this -# OSX_LDFLAGS := -lstdc++ - -# 9. Check track state propagation for success, turns on simple -# checks of filter convergence: used in SMatrix code mostly, still retain as toyMC propagation still uses this -USE_STATE_VALIDITY_CHECKS := -DCHECKSTATEVALID - -# 10. Turn on multiple scattering: for toyMC SMatrix code. Scattering handled through material map in CMSSW -#USE_SCATTERING := -DSCATTERING - -# 11. In SMatrix track building, use linear interpolation across a -# a volume instead of using the geometry -#USE_LINEAR_INTERPOLATION := -DLINEARINTERP - -# 12. Use built tracks for fitting in SMatrix code, comment out to fit sim tracks -#ENDTOEND := -DENDTOEND - -# 13. Intel Threading Building Blocks. With icc uses system +# 6. Intel Threading Building Blocks. With icc uses system # TBB, otherwise requires a local installation, and paths will # have to be adjusted. WITH_TBB := 1 -# 14. Use inward fit in Conformal fit + final KF Fit: unsed in mkFit, used in SMatrix +# 7. Use inward fit in Conformal fit #INWARD_FIT := -DINWARDFIT ################################################################ @@ -111,7 +89,14 @@ CXXFLAGS := -fPIC ${OPT} ${OSX_CXXFLAGS} LDFLAGS_HOST := -CPPFLAGS += ${USE_STATE_VALIDITY_CHECKS} ${USE_SCATTERING} ${USE_LINEAR_INTERPOLATION} ${ENDTOEND} ${INWARD_FIT} +CPPFLAGS += ${INWARD_FIT} + +# Dump hit selection window tuples +CPPFLAGS += -DDUMPHITWINDOW # -DDUMP_HIT_PRESEL +# The new dumpers +ifdef WITH_ROOT + CPPFLAGS += -DRNT_DUMP_MkF_SelHitIdcs +endif ifdef USE_VTUNE_NOTIFY ifdef VTUNE_AMPLIFIER_XE_2017_DIR @@ -129,11 +114,13 @@ endif ifeq ($(CXX), g++) # For gcc, compile with flags approximating "scram build echo_CXXFLAGS" in a CMSSW env - CXXFLAGS += -pthread -pipe -Werror=main -Werror=pointer-arith -Werror=overlength-strings -Wno-vla -Werror=overflow -std=c++17 -ftree-vectorize -Werror=array-bounds -Werror=format-contains-nul -Werror=type-limits -fvisibility-inlines-hidden -fno-math-errno --param vect-max-version-for-alias-checks=50 -Xassembler --compress-debug-sections -fuse-ld=bfd -msse3 -felide-constructors -fmessage-length=0 -Wall -Wno-non-template-friend -Wno-long-long -Wreturn-type -Wextra -Wpessimizing-move -Wclass-memaccess -Wno-cast-function-type -Wno-unused-but-set-parameter -Wno-ignored-qualifiers -Wno-deprecated-copy -Wno-unused-parameter -Wunused -Wparentheses -Wno-deprecated -Werror=return-type -Werror=missing-braces -Werror=unused-value -Werror=unused-label -Werror=address -Werror=format -Werror=sign-compare -Werror=write-strings -Werror=delete-non-virtual-dtor -Werror=strict-aliasing -Werror=narrowing -Werror=unused-but-set-variable -Werror=reorder -Werror=unused-variable -Werror=conversion-null -Werror=return-local-addr -Wnon-virtual-dtor -Werror=switch -fdiagnostics-show-option -Wno-unused-local-typedefs -Wno-attributes -Wno-psabi + # -msse3 is enabled above, if so configured (i.e., not avx) + # -Ofast is specified above in section 3. + # -std is set to c++1z below (was c++17 here) + CXXFLAGS += -pthread -pipe -Werror=main -Werror=pointer-arith -Werror=overlength-strings -Wno-vla -Werror=overflow -ftree-vectorize -Werror=array-bounds -Werror=format-contains-nul -Werror=type-limits -fvisibility-inlines-hidden -fno-math-errno --param vect-max-version-for-alias-checks=50 -Xassembler --compress-debug-sections -fuse-ld=bfd -felide-constructors -fmessage-length=0 -Wall -Wno-non-template-friend -Wno-long-long -Wreturn-type -Wextra -Wpessimizing-move -Wclass-memaccess -Wno-cast-function-type -Wno-unused-but-set-parameter -Wno-ignored-qualifiers -Wno-deprecated-copy -Wno-unused-parameter -Wunused -Wparentheses -Wno-deprecated -Werror=return-type -Werror=missing-braces -Werror=unused-value -Werror=unused-label -Werror=address -Werror=format -Werror=sign-compare -Werror=write-strings -Werror=delete-non-virtual-dtor -Werror=strict-aliasing -Werror=narrowing -Werror=unused-but-set-variable -Werror=reorder -Werror=unused-variable -Werror=conversion-null -Werror=return-local-addr -Wnon-virtual-dtor -Werror=switch -fdiagnostics-show-option -Wno-unused-local-typedefs -Wno-attributes -Wno-psabi CXXFLAGS += -fdiagnostics-color=auto -fdiagnostics-show-option -fopenmp-simd - # Add flags to enable auto-vectorization of sin/cos with libmvec (glibc >= 2.22) # Disable reciprocal math for Intel/AMD consistency, see cms-sw/cmsdist#8280 - CXXFLAGS += -Ofast -fno-reciprocal-math -mrecip=none + CXXFLAGS += -fno-reciprocal-math -mrecip=none endif # Try to find a new enough TBB @@ -172,8 +159,8 @@ ifdef DEBUG CPPFLAGS += -DDEBUG endif -# Set stdlib at the very end, as other flags (i.e. ROOT) can override our choice for which version of c++ -CPPFLAGS += -std=c++1z +# Set C++ standard version at the very end, as other flags (i.e. ROOT) can override our choice for which version of c++ +CXXFLAGS += -std=c++1z ################################################################ # Dependency generation diff --git a/RecoTracker/MkFitCore/standalone/RntDumper/MkFinder_selectHitIndices.icc b/RecoTracker/MkFitCore/standalone/RntDumper/MkFinder_selectHitIndices.icc new file mode 100644 index 0000000000000..e5af09b9be68b --- /dev/null +++ b/RecoTracker/MkFitCore/standalone/RntDumper/MkFinder_selectHitIndices.icc @@ -0,0 +1,159 @@ +// Include-cc file for ntuple dumpers in MkFinder::selectHitIndices +// To be included in cc file, outside of any namespace. + +#include "RecoTracker/MkFitCore/standalone/RntDumper/RntDumper.h" +#include "RecoTracker/MkFitCore/standalone/RntDumper/RntStructs.h" + +#include "RecoTracker/MkFitCore/src/Matrix.h" +#include "RecoTracker/MkFitCore/interface/Track.h" +#include "RecoTracker/MkFitCore/src/MiniPropagators.h" + +#include +#include + +#include + +namespace { + namespace miprops = mkfit::mini_propagators; + using namespace mkfit; + + RVec state2pos(const miprops::State &s) { return {s.x, s.y, s.z}; } + RVec state2mom(const miprops::State &s) { return {s.px, s.py, s.pz}; } + State state2state(const miprops::State &s) { return {state2pos(s), state2mom(s)}; } + + RVec statep2pos(const miprops::StatePlex &s, int i) { return {s.x[i], s.y[i], s.z[i]}; } + RVec statep2mom(const miprops::StatePlex &s, int i) { return {s.px[i], s.py[i], s.pz[i]}; } + State statep2state(const miprops::StatePlex &s, int i) { return {statep2pos(s, i), statep2mom(s, i)}; } + PropState statep2propstate(const miprops::StatePlex &s, int i) { + return {statep2state(s, i), s.dalpha[i], s.fail_flag[i]}; + } + + RVec track2pos(const TrackBase &s) { return {s.x(), s.y(), s.z()}; } + RVec track2mom(const TrackBase &s) { return {s.px(), s.py(), s.pz()}; } + State track2state(const TrackBase &s) { return {track2pos(s), track2mom(s)}; } + + SimSeedInfo evsi2ssinfo(const Event *ev, int seed_idx) { + SimSeedInfo ssi; + Event::SimLabelFromHits slfh = ev->simLabelForCurrentSeed(seed_idx); + if (slfh.is_set()) { + ssi.s_sim = track2state(ev->simTracks_[slfh.label]); + ssi.sim_lbl = slfh.label; + ssi.n_hits = slfh.n_hits; + ssi.n_match = slfh.n_match; + ssi.has_sim = true; + } + auto seed = ev->currentSeed(seed_idx); + ssi.s_seed = track2state(seed); + ssi.seed_lbl = seed.label(); + ssi.seed_idx = seed_idx; + return ssi; + } + + // --- === --- + + struct RntIfc_selectHitIndices { + using RNTupleWriter = ROOT::Experimental::RNTupleWriter; + + RntDumper *f_dumper; + const bool f_do_rnt; + const bool f_do_tree; + std::vector f_h_idcs; + std::vector f_h_remap; + int f_h_cnt; + + RNTupleWriter *f_H_writer = nullptr; + TTree *f_H_tree = nullptr; + std::shared_ptr h; + std::shared_ptr> ci; + + RNTupleWriter *f_F_writer = nullptr; + TTree *f_F_tree = nullptr; + std::shared_ptr f; + std::shared_ptr> fpi; + + RntIfc_selectHitIndices(bool rntp, bool treep) : f_do_rnt(rntp), f_do_tree(treep) { + f_dumper = RntDumper::Create("SelHitIdcs.root"); + + auto mh = f_dumper->CreateModel(); + // Register branches -- entry objects used for both RNTuple and TTree! + h = mh->MakeField("h"); + ci = mh->MakeField>("ci"); + if (f_do_rnt) { + f_H_writer = f_dumper->WritifyModel(mh, "H_rnt"); // setup for writing + } + if (f_do_tree) { + // printf("Addresses %p %p %p %p\n", h.get(), bso.get(), bsn.get(), lpi.get()); + f_H_tree = new TTree("H_tree", "info from selectHitIndices"); + f_H_tree->Branch("h", h.get()); + f_H_tree->Branch("ci", ci.get()); + f_dumper->RegisterTree(f_H_tree); + } + + auto mf = f_dumper->CreateModel(); + f = mf->MakeField("f"); + fpi = mf->MakeField>("fpi"); + if (f_do_rnt) { + f_F_writer = f_dumper->WritifyModel(mf, "F_rnt"); // setup for writing + } + if (f_do_tree) { + // printf("Addresses %p %p %p %p\n", h.get(), bso.get(), bsn.get(), lpi.get()); + f_F_tree = new TTree("F_tree", "info on failed propagations from selectHitIndices"); + f_F_tree->Branch("f", f.get()); + f_F_tree->Branch("fpi", fpi.get()); + f_dumper->RegisterTree(f_F_tree); + } + } + + ~RntIfc_selectHitIndices() {} + + void ResetH() { + ci->clear(); + f_h_cnt = 0; + } + void ResetF() { fpi->clear(); } + + void InnerIdcsReset(int maxN) { + f_h_idcs.clear(); + std::vector v(maxN, -666666); + f_h_remap.swap(v); + } + CandInfo &RegisterGoodProp(int i, const MPlexLV &ctr, const Event *ev, int seed_idx) { + f_h_idcs.push_back(i); + f_h_remap[i] = f_h_cnt; + ++f_h_cnt; + // for (auto i : f_h_idcs) { + // auto gi = f_h_remap[i]; + // // use i for local access, gi for access into tree vectors + // } + + mini_propagators::State c(ctr, i); + return ci->emplace_back(evsi2ssinfo(ev, seed_idx), state2state(c)); + } + void RegisterFailedProp(int i, const MPlexLV &beg, const MPlexLV &end, const Event *ev, int seed_idx) { + mini_propagators::State b(beg, i), e(end, i); + fpi->emplace_back(evsi2ssinfo(ev, seed_idx), state2state(b), state2state(e)); + } + int MapHIdx(int i) { return f_h_remap[i]; } + + void FillH() { + for (auto &e : *ci) + e.nan_check(); + if (f_do_rnt) + f_H_writer->Fill(); + if (f_do_tree) + f_H_tree->Fill(); + } + void FillF() { + if (fpi->empty()) + return; + for (auto &e : *fpi) + e.nan_check(); + if (f_do_rnt) + f_F_writer->Fill(); + if (f_do_tree) + f_F_tree->Fill(); + } + }; + + static RntIfc_selectHitIndices rnt_shi(false, true); +} // namespace diff --git a/RecoTracker/MkFitCore/standalone/RntDumper/RntDumper.cc b/RecoTracker/MkFitCore/standalone/RntDumper/RntDumper.cc new file mode 100644 index 0000000000000..f4314e670e454 --- /dev/null +++ b/RecoTracker/MkFitCore/standalone/RntDumper/RntDumper.cc @@ -0,0 +1,62 @@ +#include "RntDumper.h" + +#include +#include + +#include "TFile.h" +#include "TTree.h" + +#include + +namespace REX = ROOT::Experimental; + +std::vector RntDumper::s_instances; + +RntDumper::RntDumper(const char *fname) : m_file(TFile::Open(fname, "recreate")) { + if (!m_file || !m_file->IsOpen()) { + printf("RntDumper::RntDumper() failed creeating file '%s'.\n", fname); + throw std::runtime_error("Failed creating file"); + } + printf("RntDumper::RntDumper() succesfully opened file '%s'.\n", fname); +} + +RntDumper::~RntDumper() { + printf("RntDumper::~RntDumper() destroying writers and closing file '%s'.\n", m_file->GetName()); + // Finish up trees first, ntuple-writers seem to write everything reulting + // in two cycles of trees. + for (auto &tp : m_trees) { + tp->Write(); + delete tp; + } + m_trees.clear(); + m_writers.clear(); + if (m_file) { + m_file->Close(); + } +} + +std::unique_ptr RntDumper::CreateModel() { return RNTupleModel::Create(); } + +REX::RNTupleWriter *RntDumper::WritifyModel(std::unique_ptr &model, std::string_view mname) { + auto wup = RNTupleWriter::Append(std::move(model), mname, *m_file); + REX::RNTupleWriter *w = wup.get(); + m_writers.insert({std::string(mname), std::move(wup)}); + return w; +} + +void RntDumper::RegisterTree(TTree *t) { m_trees.push_back(t); } + +// === static === + +RntDumper *RntDumper::Create(const char *fname) { + // Should check fnames ? + RntDumper *d = new RntDumper(fname); + s_instances.push_back(d); + return d; +} + +void RntDumper::FinalizeAll() { + printf("RntDumper::FinalizeAll() shutting down %d instances.\n", (int)s_instances.size()); + for (auto &d : s_instances) + delete d; +} diff --git a/RecoTracker/MkFitCore/standalone/RntDumper/RntDumper.h b/RecoTracker/MkFitCore/standalone/RntDumper/RntDumper.h new file mode 100644 index 0000000000000..c2478f4080664 --- /dev/null +++ b/RecoTracker/MkFitCore/standalone/RntDumper/RntDumper.h @@ -0,0 +1,43 @@ +#ifndef RecoTracker_MkFitCore_standalone_RntDumper_RntDumper_h +#define RecoTracker_MkFitCore_standalone_RntDumper_RntDumper_h + +#include +#include +#include +#include + +class TFile; +class TTree; + +namespace ROOT::Experimental { + class RNTupleModel; + class RNTupleWriter; +} // namespace ROOT::Experimental + +class RntDumper { + using RNTupleModel = ROOT::Experimental::RNTupleModel; + using RNTupleWriter = ROOT::Experimental::RNTupleWriter; + +public: + std::unique_ptr CreateModel(); + RNTupleWriter *WritifyModel(std::unique_ptr &model, std::string_view mname); + + void RegisterTree(TTree *t); + + static RntDumper *Create(const char *fname); + static void FinalizeAll(); + + TFile *file() { return m_file.get(); } + +private: + explicit RntDumper(const char *fname); + ~RntDumper(); + + std::unique_ptr m_file; + std::unordered_map> m_writers; + std::vector m_trees; + + static std::vector s_instances; +}; + +#endif diff --git a/RecoTracker/MkFitCore/standalone/RntDumper/RntStructs.cc b/RecoTracker/MkFitCore/standalone/RntDumper/RntStructs.cc new file mode 100644 index 0000000000000..166f5c5886a43 --- /dev/null +++ b/RecoTracker/MkFitCore/standalone/RntDumper/RntStructs.cc @@ -0,0 +1,53 @@ +#include "RntStructs.h" + +#include + +namespace { + constexpr bool isFinite(float x) { + const unsigned int mask = 0x7f800000; + union { + unsigned int l; + float d; + } v = {.d = x}; + return (v.l & mask) != mask; + } + + // nan-guard, in place, return true if nan detected. + bool ngr(float &f) { + bool is_bad = !isFinite(f); + if (is_bad) + f = -999.0f; + return is_bad; + } + bool ngr(RVec &v) { + bool is_bad = ngr(v.fX); + is_bad |= ngr(v.fY); + is_bad |= ngr(v.fZ); + return is_bad; + } + bool ngr(State &s) { + bool is_bad = ngr(s.pos); + is_bad |= ngr(s.mom); + return is_bad; + } +} // namespace + +bool BinSearch::nan_check() { + has_nans = ngr(phi); + has_nans |= ngr(dphi); + has_nans |= ngr(q); + has_nans |= ngr(dq); + return has_nans; +} + +void CandInfo::nan_check() { + has_nans = ngr(ps_min); + has_nans |= ngr(ps_max); + has_nans |= bso.nan_check(); + has_nans |= bsn.nan_check(); +} + +void FailedPropInfo::nan_check() { + has_nans = ngr(s_prev); + has_nans |= ngr(s_final); +} diff --git a/RecoTracker/MkFitCore/standalone/RntDumper/RntStructs.h b/RecoTracker/MkFitCore/standalone/RntDumper/RntStructs.h new file mode 100644 index 0000000000000..35d6dcacb68ac --- /dev/null +++ b/RecoTracker/MkFitCore/standalone/RntDumper/RntStructs.h @@ -0,0 +1,104 @@ +#ifndef RecoTracker_MkFitCore_standalone_RntDumper_RntStructs_h +#define RecoTracker_MkFitCore_standalone_RntDumper_RntStructs_h + +// Avoid MkFit includes for now to simpligy pure ROOT builds. +// #include "RecoTracker/MkFitCore/interface/ + +#include "ROOT/REveVector.hxx" +#include "Math/Point3D.h" +#include "Math/Vector3D.h" + +// From CMSSW data formats +/// point in space with cartesian internal representation +typedef ROOT::Math::PositionVector3D > XYZPointF; +/// spatial vector with cartesian internal representation +typedef ROOT::Math::DisplacementVector3D > XYZVectorF; +/// spatial vector with cylindrical internal representation using pseudorapidity +typedef ROOT::Math::DisplacementVector3D > RhoEtaPhiVectorF; +/// spatial vector with polar internal representation +/// WARNING: ROOT dictionary not provided for the type below +// typedef ROOT::Math::DisplacementVector3D > RThetaPhiVectorF; + +using RVec = ROOT::Experimental::REveVector; + +struct HeaderLayer { + int event, iter_idx, iter_algo, eta_region, layer; + float qb_min, qb_max; // qbar layer limits, r for barrel, z for endcap + bool is_barrel, is_pix, is_stereo; + + HeaderLayer() = default; + HeaderLayer& operator=(const HeaderLayer&) = default; +}; + +struct State { + RVec pos, mom; + + State() = default; + State& operator=(const State&) = default; +}; + +struct PropState : public State { + float dalpha; // helix angle during propagation + int fail_flag; + + PropState() = default; + PropState& operator=(const PropState&) = default; +}; + +struct SimSeedInfo { + State s_sim; + State s_seed; + int sim_lbl, seed_lbl, seed_idx; + int n_hits, n_match; + bool has_sim = false; + + float good_frac() const { return (float)n_match / n_hits; } + + SimSeedInfo() = default; + SimSeedInfo& operator=(const SimSeedInfo&) = default; +}; + +struct BinSearch { + float phi, dphi, q, dq; + short unsigned int p1, p2, q1, q2; + short int wsr; + bool wsr_in_gap; + bool has_nans = false; + + bool nan_check(); + + BinSearch() = default; + BinSearch& operator=(const BinSearch&) = default; +}; + +struct CandInfo { + SimSeedInfo ssi; + State s_ctr; + PropState ps_min, ps_max; + BinSearch bso; + BinSearch bsn; + bool has_nans = false; + + CandInfo(const SimSeedInfo& s, const State& c) : ssi(s), s_ctr(c) {} + + void nan_check(); + + CandInfo() = default; + CandInfo& operator=(const CandInfo&) = default; +}; + +struct FailedPropInfo { + SimSeedInfo ssi; + State s_prev; + State s_final; + bool has_nans = false; + + FailedPropInfo(const SimSeedInfo& s, const State& p, const State& f) : ssi(s), s_prev(p), s_final(f) {} + + void nan_check(); + + FailedPropInfo() = default; + FailedPropInfo& operator=(const FailedPropInfo&) = default; +}; + +#endif diff --git a/RecoTracker/MkFitCore/standalone/RntDumper/RntStructs_Linkdef.h b/RecoTracker/MkFitCore/standalone/RntDumper/RntStructs_Linkdef.h new file mode 100644 index 0000000000000..da56983314ab3 --- /dev/null +++ b/RecoTracker/MkFitCore/standalone/RntDumper/RntStructs_Linkdef.h @@ -0,0 +1,14 @@ +// #pragma link C++ class ROOT::Experimental::REveVector+; + +#pragma link C++ class HeaderLayer + ; + +#pragma link C++ class State + ; +#pragma link C++ class PropState + ; +#pragma link C++ class SimSeedInfo + ; +#pragma link C++ class BinSearch + ; + +#pragma link C++ class CandInfo + ; +#pragma link C++ class std::vector < CandInfo> + ; + +#pragma link C++ class FailedPropInfo + ; +#pragma link C++ class std::vector < FailedPropInfo> + ;