diff --git a/RecoTracker/MkFit/README.md b/RecoTracker/MkFit/README.md index 831a5bce698c3..b5bf82ce03275 100644 --- a/RecoTracker/MkFit/README.md +++ b/RecoTracker/MkFit/README.md @@ -47,12 +47,41 @@ $ runTheMatrix.py -l --apply 2 --command "--procModifiers tracking * *m_track_algorithm:* CMSSW track algorithm (used internally for reporting and consistency checks) * *m_requires_seed_hit_sorting:* do hits on seed tracks need to be sorted (required for seeds that include strip layers) -* *m_requires_quality_filter:* is additional post-processing required for result tracks -* *m_requires_dupclean_tight:* is tight duplicate removal post-processing required for result tracks * *m_params:* IterationParams structure for this iteration * *m_backward_params:* IterationParams structure for backward search for this iteration * *m_layer_configs:* std::vector of per-layer parameters +#### Seed cleaning params (based on elliptical dR-dz cut) [in class IterationConfig] + +* *m_seed_cleaner_name:* name of standard function to call for seed cleaning; if not set or empty seed cleaning is not performed +* *sc_ptthr_hpt:* pT threshold used to tighten seed cleaning requirements +* *sc_drmax_bh:* dR cut used for seed tracks with std::fabs(eta)<0.9 and pT > c_ptthr_hpt +* *sc_dzmax_bh:* dz cut used for seed tracks with std::fabs(eta)<0.9 and pT > c_ptthr_hpt +* *sc_drmax_eh:* dR cut used for seed tracks with std::fabs(eta)>0.9 and pT > c_ptthr_hpt +* *sc_dzmax_eh:* dz cut used for seed tracks with std::fabs(eta)>0.9 and pT > c_ptthr_hpt +* *sc_drmax_bl:* dR cut used for seed tracks with std::fabs(eta)<0.9 and pT < c_ptthr_hpt +* *sc_dzmax_bl:* dz cut used for seed tracks with std::fabs(eta)<0.9 and pT < c_ptthr_hpt +* *sc_drmax_el:* dR cut used for seed tracks with std::fabs(eta)>0.9 and pT < c_ptthr_hpt +* *sc_dzmax_el:* dz cut used for seed tracks with std::fabs(eta)>0.9 and pT < c_ptthr_hpt + +#### Seed partitioning params [in class IterationConfig] + +* *m_seed_partitioner_name:* name of standard function to call for seed partitioning + +#### Pre / post backward-fit candidate top-level params [in class IterationConfig] + +* *m_pre_bkfit_filter_name:* name of standard function used for candidate filtering after forward +search but before backward fit; if not set or empty no candidate filtering is performed at this stage +* *m_post_bkfit_filter_name:* name of standard function used for candidate filtering after backward fit; if not set or empty no candidate filtering is performed at this stage + +#### Duplicate cleaning parameters [in class IterationConfig] + +* *m_duplicate_cleaner_name:* name of standard function used for duplicate cleaning; if not set or empty duplicate cleaning is not performed +* *dc_fracSharedHits:* min fraction of shared hits to determine duplicate track candidate +* *dc_drth_central:* dR cut used to identify duplicate candidates if std::abs(cotan(theta))<1.99 (abs(eta)<1.44) +* *dc_drth_obarrel:* dR cut used to identify duplicate candidates if 1.996.05 (abs(eta)>2.5) + ### Iteration parameters [class IterationParams] * *nlayers_per_seed:* internal mkFit parameter used for standalone validation @@ -63,25 +92,8 @@ $ runTheMatrix.py -l --apply 2 --command "--procModifiers tracking * *chi2CutOverlap:* chi2 cut for accepting an overlap hit (currently NOT used) * *pTCutOverlap:* pT cut below which the overlap hits are not picked up -#### Seed cleaning params (based on elliptical dR-dz cut) - -* *c_ptthr_hpt:* pT threshold used to tighten seed cleaning requirements -* *c_drmax_bh:* dR cut used for seed tracks with std::fabs(eta)<0.9 and pT > c_ptthr_hpt -* *c_dzmax_bh:* dz cut used for seed tracks with std::fabs(eta)<0.9 and pT > c_ptthr_hpt -* *c_drmax_eh:* dR cut used for seed tracks with std::fabs(eta)>0.9 and pT > c_ptthr_hpt -* *c_dzmax_eh:* dz cut used for seed tracks with std::fabs(eta)>0.9 and pT > c_ptthr_hpt -* *c_drmax_bl:* dR cut used for seed tracks with std::fabs(eta)<0.9 and pT < c_ptthr_hpt -* *c_dzmax_bl:* dz cut used for seed tracks with std::fabs(eta)<0.9 and pT < c_ptthr_hpt -* *c_drmax_el:* dR cut used for seed tracks with std::fabs(eta)>0.9 and pT < c_ptthr_hpt -* *c_dzmax_el:* dz cut used for seed tracks with std::fabs(eta)>0.9 and pT < c_ptthr_hpt - -#### Duplicate cleaning parameters - -* *minHitsQF:* min number of hits on track candidate to apply duplicate cleaning based on fraction of shared hits -* *fracSharedHits:* min fraction of shared hits to determine duplicate track candidate -* *drth_central:* dR cut used to identify duplicate candidates if std::abs(cotan(theta))<1.99 (abs(eta)<1.44) -* *drth_obarrel:* dR cut used to identify duplicate candidates if 1.996.05 (abs(eta)>2.5) +#### Pre / post backward-fit candidate filtering params +* *minHitsQF:* minimum number of hits, interpretation depends on particular filtering function used ### Per-layer parameters [class IterationLayerConfig] diff --git a/RecoTracker/MkFit/interface/MkFitGeometry.h b/RecoTracker/MkFit/interface/MkFitGeometry.h index bd2f73a4f3fb5..29ba8c392b23e 100644 --- a/RecoTracker/MkFit/interface/MkFitGeometry.h +++ b/RecoTracker/MkFit/interface/MkFitGeometry.h @@ -25,7 +25,8 @@ class MkFitGeometry { explicit MkFitGeometry(const TrackerGeometry& geom, const GeometricSearchTracker& tracker, const TrackerTopology& ttopo, - std::unique_ptr trackerInfo); + std::unique_ptr trackerInfo, + const mkfit::LayerNumberConverter& layNConv); ~MkFitGeometry(); int mkFitLayerNumber(DetId detId) const; diff --git a/RecoTracker/MkFit/plugins/MkFitGeometryESProducer.cc b/RecoTracker/MkFit/plugins/MkFitGeometryESProducer.cc index b8c88a713a291..172460b17c7bc 100644 --- a/RecoTracker/MkFit/plugins/MkFitGeometryESProducer.cc +++ b/RecoTracker/MkFit/plugins/MkFitGeometryESProducer.cc @@ -21,6 +21,8 @@ #include +// #define DUMP_MKF_GEO + //------------------------------------------------------------------------------ class MkFitGeometryESProducer : public edm::ESProducer { @@ -192,7 +194,9 @@ void MkFitGeometryESProducer::fillShapeAndPlacement(const GeomDet *det, xy[3][1] = -par[3]; dz = par[2]; - // printf("TRAP 0x%x %f %f %f %f\n", detid.rawId(), par[0], par[1], par[2], par[3]); +#ifdef DUMP_MKF_GEO + printf("TRAP 0x%x %f %f %f %f ", detid.rawId(), par[0], par[1], par[2], par[3]); +#endif } else if (const RectangularPlaneBounds *b2 = dynamic_cast(b)) { // Rectangular float dx = b2->width() * 0.5; // half width @@ -207,7 +211,9 @@ void MkFitGeometryESProducer::fillShapeAndPlacement(const GeomDet *det, xy[3][1] = -dy; dz = b2->thickness() * 0.5; // half thickness - // printf("RECT 0x%x %f %f %f\n", detid.rawId(), dx, dy, dz); +#ifdef DUMP_MKF_GEO + printf("RECT 0x%x %f %f %f ", detid.rawId(), dx, dy, dz); +#endif } else { throw cms::Exception("UnimplementedFeature") << "unsupported Bounds class"; } @@ -219,6 +225,14 @@ void MkFitGeometryESProducer::fillShapeAndPlacement(const GeomDet *det, useMatched, trackerTopo_->isStereo(detid), trackerTopo_->side(detid) == static_cast(TrackerDetSide::PosEndcap)); +#ifdef DUMP_MKF_GEO + printf(" subdet=%d layer=%d side=%d is_stereo=%d --> mkflayer=%d\n", + detid.subdetId(), + trackerTopo_->layer(detid), + trackerTopo_->side(detid), + trackerTopo_->isStereo(detid), + lay); +#endif mkfit::LayerInfo &layer_info = trk_info.layer_nc(lay); if (lgc_map) { @@ -252,51 +266,64 @@ void MkFitGeometryESProducer::fillShapeAndPlacement(const GeomDet *det, //============================================================================== -// Ideally these functions would also: -// 0. Setup LayerInfo data (which is now done in auto-generated code). -// Some data-members are a bit over specific, esp/ bools for CMS sub-detectors. -// 1. Establish short module ids (now done in MkFitGeometry constructor). -// 2. Store module normal and strip direction vectors -// 3. ? Any other information ? +// These functions do the following: +// 0. Detect bounding cylinder of each layer. +// 1. Setup LayerInfo data. +// 2. Establish short module ids. +// 3. Store module normal and strip direction vectors. // 4. Extract stereo coverage holes where they exist (TEC, all but last 3 double-layers). // -// Plugin DumpMkFitGeometry.cc can then be used to export this for stand-alone. -// Would also need to be picked up with tk-ntuple converter (to get module ids as -// they will now be used as indices into module info vectors). -// -// An attempt at export cmsRun config is in python/dumpMkFitGeometry.py +// See python/dumpMkFitGeometry.py and dumpMkFitGeometryPhase2.py void MkFitGeometryESProducer::addPixBGeometry(mkfit::TrackerInfo &trk_info) { +#ifdef DUMP_MKF_GEO + printf("\n*** addPixBGeometry\n\n"); +#endif for (auto &det : trackerGeom_->detsPXB()) { fillShapeAndPlacement(det, trk_info); } } void MkFitGeometryESProducer::addPixEGeometry(mkfit::TrackerInfo &trk_info) { +#ifdef DUMP_MKF_GEO + printf("\n*** addPixEGeometry\n\n"); +#endif for (auto &det : trackerGeom_->detsPXF()) { fillShapeAndPlacement(det, trk_info); } } void MkFitGeometryESProducer::addTIBGeometry(mkfit::TrackerInfo &trk_info) { +#ifdef DUMP_MKF_GEO + printf("\n*** addTIBGeometry\n\n"); +#endif for (auto &det : trackerGeom_->detsTIB()) { fillShapeAndPlacement(det, trk_info); } } void MkFitGeometryESProducer::addTOBGeometry(mkfit::TrackerInfo &trk_info) { +#ifdef DUMP_MKF_GEO + printf("\n*** addTOBGeometry\n\n"); +#endif for (auto &det : trackerGeom_->detsTOB()) { fillShapeAndPlacement(det, trk_info); } } void MkFitGeometryESProducer::addTIDGeometry(mkfit::TrackerInfo &trk_info) { +#ifdef DUMP_MKF_GEO + printf("\n*** addTIDGeometry\n\n"); +#endif for (auto &det : trackerGeom_->detsTID()) { fillShapeAndPlacement(det, trk_info); } } void MkFitGeometryESProducer::addTECGeometry(mkfit::TrackerInfo &trk_info) { +#ifdef DUMP_MKF_GEO + printf("\n*** addTECGeometry\n\n"); +#endif // For TEC we also need to discover hole in radial extents. layer_gap_map_t lgc_map; for (auto &det : trackerGeom_->detsTEC()) { @@ -330,6 +357,15 @@ namespace { // PIXE-, TID-, TEC- 1.0, 1.0, 1.0, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5, 10.25, 7.5 }; + const float phase2QBins[] = { + // TODO: Review these numbers. + // PIXB, TOB + 2.0, 2.0, 2.0, 2.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, + // PIXE+, TEC+ + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, + // PIXE-, TEC- + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6, 5.6 + }; } // clang-format on //------------------------------------------------------------------------------ @@ -340,13 +376,19 @@ std::unique_ptr MkFitGeometryESProducer::produce(const TrackerRec trackerGeom_ = &iRecord.get(geomToken_); trackerTopo_ = &iRecord.get(ttopoToken_); + const float *qBinDefaults = nullptr; + // std::string path = "Geometry/TrackerCommonData/data/"; if (trackerGeom_->isThere(GeomDetEnumerators::P1PXB) || trackerGeom_->isThere(GeomDetEnumerators::P1PXEC)) { edm::LogInfo("MkFitGeometryESProducer") << "Extracting PhaseI geometry"; trackerInfo->create_layers(18, 27, 27); + qBinDefaults = phase1QBins; } else if (trackerGeom_->isThere(GeomDetEnumerators::P2PXB) || trackerGeom_->isThere(GeomDetEnumerators::P2PXEC) || trackerGeom_->isThere(GeomDetEnumerators::P2OTB) || trackerGeom_->isThere(GeomDetEnumerators::P2OTEC)) { - throw cms::Exception("UnimplementedFeature") << "PhaseII geometry extraction"; + edm::LogInfo("MkFitGeometryESProducer") << "Extracting PhaseII geometry"; + layerNrConv_.reset(mkfit::TkLayout::phase2); + trackerInfo->create_layers(16, 22, 22); + qBinDefaults = phase2QBins; } else { throw cms::Exception("UnimplementedFeature") << "unsupported / unknowen geometry version"; } @@ -358,7 +400,7 @@ std::unique_ptr MkFitGeometryESProducer::produce(const TrackerRec std::numeric_limits::max(), 0, std::numeric_limits::max(), -std::numeric_limits::max()); li.reserve_modules(256); } - // This is sort of CMS-2017 specific ... but fireworks code uses it for PhaseII as well. + // This is sort of CMS-phase1 specific ... but fireworks code uses it for PhaseII as well. addPixBGeometry(*trackerInfo); addPixEGeometry(*trackerInfo); addTIBGeometry(*trackerInfo); @@ -366,19 +408,30 @@ std::unique_ptr MkFitGeometryESProducer::produce(const TrackerRec addTOBGeometry(*trackerInfo); addTECGeometry(*trackerInfo); - // r_in/out kept as squres until here, root them + // r_in/out kept as squares until here, root them + unsigned int n_mod = 0; for (int i = 0; i < trackerInfo->n_layers(); ++i) { auto &li = trackerInfo->layer_nc(i); li.set_r_in_out(std::sqrt(li.rin()), std::sqrt(li.rout())); li.set_propagate_to(li.is_barrel() ? li.r_mean() : li.z_mean()); - li.set_q_bin(phase1QBins[i]); + li.set_q_bin(qBinDefaults[i]); unsigned int maxsid = li.shrink_modules(); - // Make sure the short id fits in the 12 bits... - assert(maxsid < 1u << 11); - } - return std::make_unique( - iRecord.get(geomToken_), iRecord.get(trackerToken_), iRecord.get(ttopoToken_), std::move(trackerInfo)); + n_mod += maxsid; + + // Make sure the short id fits in the 14 bits... + assert(maxsid < 1u << 13); + assert(n_mod > 0); + } +#ifdef DUMP_MKF_GEO + printf("Total number of modules %u, 14-bits fit up to %u modules\n", n_mod, 1u << 13); +#endif + + return std::make_unique(iRecord.get(geomToken_), + iRecord.get(trackerToken_), + iRecord.get(ttopoToken_), + std::move(trackerInfo), + layerNrConv_); } DEFINE_FWK_EVENTSETUP_MODULE(MkFitGeometryESProducer); diff --git a/RecoTracker/MkFit/plugins/MkFitIterationConfigESProducer.cc b/RecoTracker/MkFit/plugins/MkFitIterationConfigESProducer.cc index 186349c026d33..05012104d2895 100644 --- a/RecoTracker/MkFit/plugins/MkFitIterationConfigESProducer.cc +++ b/RecoTracker/MkFit/plugins/MkFitIterationConfigESProducer.cc @@ -6,197 +6,8 @@ #include "RecoTracker/MkFit/interface/MkFitGeometry.h" // mkFit includes -#include "RecoTracker/MkFitCore/interface/Track.h" -#include "RecoTracker/MkFitCore/interface/TrackerInfo.h" -#include "RecoTracker/MkFitCore/interface/HitStructures.h" #include "RecoTracker/MkFitCore/interface/IterationConfig.h" -namespace { - using namespace mkfit; - - [[maybe_unused]] void partitionSeeds0(const TrackerInfo &trk_info, - const TrackVec &in_seeds, - const EventOfHits &eoh, - IterationSeedPartition &part) { - const size_t size = in_seeds.size(); - - for (size_t i = 0; i < size; ++i) { - const Track &S = in_seeds[i]; - - const bool z_dir_pos = S.pz() > 0; - - const auto &hot = S.getLastHitOnTrack(); - const float eta = eoh[hot.layer].refHit(hot.index).eta(); - - // Region to be defined by propagation / intersection tests - TrackerInfo::EtaRegion reg; - - const LayerInfo &outer_brl = trk_info.outer_barrel_layer(); - - // Define first (mkFit) layer IDs for each strip subdetector. - constexpr int tib1_id = 4; - constexpr int tob1_id = 10; - constexpr int tecp1_id = 27; - constexpr int tecn1_id = 54; - - const LayerInfo &tib1 = trk_info.layer(tib1_id); - const LayerInfo &tob1 = trk_info.layer(tob1_id); - - const LayerInfo &tecp1 = trk_info.layer(tecp1_id); - const LayerInfo &tecn1 = trk_info.layer(tecn1_id); - - const LayerInfo &tec_first = z_dir_pos ? tecp1 : tecn1; - - const float maxR = S.maxReachRadius(); - float z_at_maxr; - - bool can_reach_outer_brl = S.canReachRadius(outer_brl.rout()); - float z_at_outer_brl; - bool misses_first_tec; - if (can_reach_outer_brl) { - z_at_outer_brl = S.zAtR(outer_brl.rout()); - if (z_dir_pos) - misses_first_tec = z_at_outer_brl < tec_first.zmin(); - else - misses_first_tec = z_at_outer_brl > tec_first.zmax(); - } else { - z_at_maxr = S.zAtR(maxR); - if (z_dir_pos) - misses_first_tec = z_at_maxr < tec_first.zmin(); - else - misses_first_tec = z_at_maxr > tec_first.zmax(); - } - - if (misses_first_tec) { - reg = TrackerInfo::Reg_Barrel; - } else { - if ((S.canReachRadius(tib1.rin()) && tib1.is_within_z_limits(S.zAtR(tib1.rin()))) || - (S.canReachRadius(tob1.rin()) && tob1.is_within_z_limits(S.zAtR(tob1.rin())))) { - reg = z_dir_pos ? TrackerInfo::Reg_Transition_Pos : TrackerInfo::Reg_Transition_Neg; - } else { - reg = z_dir_pos ? TrackerInfo::Reg_Endcap_Pos : TrackerInfo::Reg_Endcap_Neg; - } - } - - part.m_region[i] = reg; - if (part.m_phi_eta_foo) - part.m_phi_eta_foo(eoh[hot.layer].refHit(hot.index).phi(), eta); - } - } - - [[maybe_unused]] void partitionSeeds1(const TrackerInfo &trk_info, - const TrackVec &in_seeds, - const EventOfHits &eoh, - IterationSeedPartition &part) { - // Define first (mkFit) layer IDs for each strip subdetector. - constexpr int tib1_id = 4; - constexpr int tob1_id = 10; - constexpr int tidp1_id = 21; - constexpr int tidn1_id = 48; - constexpr int tecp1_id = 27; - constexpr int tecn1_id = 54; - - const LayerInfo &tib1 = trk_info.layer(tib1_id); - const LayerInfo &tob1 = trk_info.layer(tob1_id); - - const LayerInfo &tidp1 = trk_info.layer(tidp1_id); - const LayerInfo &tidn1 = trk_info.layer(tidn1_id); - - const LayerInfo &tecp1 = trk_info.layer(tecp1_id); - const LayerInfo &tecn1 = trk_info.layer(tecn1_id); - - // Merge first two layers to account for mono/stereo coverage. - // TrackerInfo could hold joint limits for sub-detectors. - const auto &L = trk_info; - const float tidp_rin = std::min(L[tidp1_id].rin(), L[tidp1_id + 1].rin()); - const float tidp_rout = std::max(L[tidp1_id].rout(), L[tidp1_id + 1].rout()); - const float tecp_rin = std::min(L[tecp1_id].rin(), L[tecp1_id + 1].rin()); - const float tecp_rout = std::max(L[tecp1_id].rout(), L[tecp1_id + 1].rout()); - const float tidn_rin = std::min(L[tidn1_id].rin(), L[tidn1_id + 1].rin()); - const float tidn_rout = std::max(L[tidn1_id].rout(), L[tidn1_id + 1].rout()); - const float tecn_rin = std::min(L[tecn1_id].rin(), L[tecn1_id + 1].rin()); - const float tecn_rout = std::max(L[tecn1_id].rout(), L[tecn1_id + 1].rout()); - - // Bias towards more aggressive transition-region assignemnts. - // With current tunning it seems to make things a bit worse. - const float tid_z_extra = 0.0f; // 5.0f; - const float tec_z_extra = 0.0f; // 10.0f; - - const size_t size = in_seeds.size(); - - auto barrel_pos_check = [](const Track &S, float maxR, float rin, float zmax) -> bool { - bool inside = maxR > rin && S.zAtR(rin) < zmax; - return inside; - }; - - auto barrel_neg_check = [](const Track &S, float maxR, float rin, float zmin) -> bool { - bool inside = maxR > rin && S.zAtR(rin) > zmin; - return inside; - }; - - auto endcap_pos_check = [](const Track &S, float maxR, float rout, float rin, float zmin) -> bool { - bool inside = maxR > rout ? S.zAtR(rout) > zmin : (maxR > rin && S.zAtR(maxR) > zmin); - return inside; - }; - - auto endcap_neg_check = [](const Track &S, float maxR, float rout, float rin, float zmax) -> bool { - bool inside = maxR > rout ? S.zAtR(rout) < zmax : (maxR > rin && S.zAtR(maxR) < zmax); - return inside; - }; - - for (size_t i = 0; i < size; ++i) { - const Track &S = in_seeds[i]; - - const auto &hot = S.getLastHitOnTrack(); - const float eta = eoh[hot.layer].refHit(hot.index).eta(); - - // Region to be defined by propagation / intersection tests - TrackerInfo::EtaRegion reg; - - const bool z_dir_pos = S.pz() > 0; - const float maxR = S.maxReachRadius(); - - if (z_dir_pos) { - const bool in_tib = barrel_pos_check(S, maxR, tib1.rin(), tib1.zmax()); - const bool in_tob = barrel_pos_check(S, maxR, tob1.rin(), tob1.zmax()); - - if (!in_tib && !in_tob) { - reg = TrackerInfo::Reg_Endcap_Pos; - } else { - const bool in_tid = endcap_pos_check(S, maxR, tidp_rout, tidp_rin, tidp1.zmin() - tid_z_extra); - const bool in_tec = endcap_pos_check(S, maxR, tecp_rout, tecp_rin, tecp1.zmin() - tec_z_extra); - - if (!in_tid && !in_tec) { - reg = TrackerInfo::Reg_Barrel; - } else { - reg = TrackerInfo::Reg_Transition_Pos; - } - } - } else { - const bool in_tib = barrel_neg_check(S, maxR, tib1.rin(), tib1.zmin()); - const bool in_tob = barrel_neg_check(S, maxR, tob1.rin(), tob1.zmin()); - - if (!in_tib && !in_tob) { - reg = TrackerInfo::Reg_Endcap_Neg; - } else { - const bool in_tid = endcap_neg_check(S, maxR, tidn_rout, tidn_rin, tidn1.zmax() + tid_z_extra); - const bool in_tec = endcap_neg_check(S, maxR, tecn_rout, tecn_rin, tecn1.zmax() + tec_z_extra); - - if (!in_tid && !in_tec) { - reg = TrackerInfo::Reg_Barrel; - } else { - reg = TrackerInfo::Reg_Transition_Neg; - } - } - } - - part.m_region[i] = reg; - if (part.m_phi_eta_foo) - part.m_phi_eta_foo(eoh[hot.layer].refHit(hot.index).phi(), eta); - } - } -} // namespace - class MkFitIterationConfigESProducer : public edm::ESProducer { public: MkFitIterationConfigESProducer(const edm::ParameterSet &iConfig); @@ -233,9 +44,9 @@ std::unique_ptr MkFitIterationConfigESProducer::produce( auto it_conf = cj.load_File(configFile_); it_conf->m_params.minPtCut = minPtCut_; it_conf->m_backward_params.minPtCut = minPtCut_; - it_conf->m_partition_seeds = partitionSeeds1; it_conf->m_params.maxClusterSize = maxClusterSize_; it_conf->m_backward_params.maxClusterSize = maxClusterSize_; + it_conf->setupStandardFunctionsFromNames(); return it_conf; } diff --git a/RecoTracker/MkFit/src/MkFitGeometry.cc b/RecoTracker/MkFit/src/MkFitGeometry.cc index 80d20bbc55b03..3096c20f178c7 100644 --- a/RecoTracker/MkFit/src/MkFitGeometry.cc +++ b/RecoTracker/MkFit/src/MkFitGeometry.cc @@ -17,14 +17,13 @@ namespace { MkFitGeometry::MkFitGeometry(const TrackerGeometry& geom, const GeometricSearchTracker& tracker, const TrackerTopology& ttopo, - std::unique_ptr trackerInfo) + std::unique_ptr trackerInfo, + const mkfit::LayerNumberConverter& layNConv) : ttopo_(&ttopo), - lnc_{std::make_unique(mkfit::TkLayout::phase1)}, + lnc_{std::make_unique(layNConv)}, trackerInfo_(std::move(trackerInfo)) { - if (geom.numberOfLayers(PixelSubdetector::PixelBarrel) != 4 || - geom.numberOfLayers(PixelSubdetector::PixelEndcap) != 3) { - throw cms::Exception("Assert") << "For now this code works only with phase1 tracker, you have something else"; - } + if (lnc_->getEra() != mkfit::TkLayout::phase1 && lnc_->getEra() != mkfit::TkLayout::phase2) + throw cms::Exception("Assert") << "This code works only with phase1 and phase2 tracker, you have something else"; // Create DetLayer structure dets_.resize(lnc_->nLayers(), nullptr); @@ -49,12 +48,9 @@ MkFitGeometry::MkFitGeometry(const TrackerGeometry& geom, const auto subdet = detId.subdetId(); const auto layer = ttopo.layer(detId); - // TODO: mono/stereo structure is still hardcoded for phase0/1 strip tracker setDet(subdet, layer, monoLayer, detId, lay); - if (((subdet == StripSubdetector::TIB or subdet == StripSubdetector::TOB) and (layer == 1 or layer == 2)) or - subdet == StripSubdetector::TID or subdet == StripSubdetector::TEC) { + if (lnc_->doesHaveStereo(subdet, layer)) setDet(subdet, layer, stereoLayer, detId, lay); - } } } diff --git a/RecoTracker/MkFit/test/dumpMkFitGeometry.py b/RecoTracker/MkFit/test/dumpMkFitGeometry.py index f5edbb2aa80f5..de6fbb9c60380 100644 --- a/RecoTracker/MkFit/test/dumpMkFitGeometry.py +++ b/RecoTracker/MkFit/test/dumpMkFitGeometry.py @@ -13,8 +13,8 @@ process.load('Configuration.StandardSequences.Reconstruction_cff') process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') -from Configuration.AlCa.GlobalTag import GlobalTag -process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase1_2022_realistic', '') +from Configuration.AlCa.autoCond import autoCond +process.GlobalTag.globaltag = autoCond['phase1_2022_realistic'] process.MessageLogger.cerr.threshold = "INFO" process.MessageLogger.cerr.MkFitGeometryESProducer = dict(limit=-1) @@ -22,7 +22,6 @@ process.source = cms.Source("EmptySource") process.maxEvents.input = 1 - process.add_(cms.ESProducer("MkFitGeometryESProducer")) defaultOutputFileName="phase1-trackerinfo.bin" @@ -30,7 +29,7 @@ # level: 0 - no printout; 1 - print layers, 2 - print modules # outputFileName: binary dump file; no dump if empty string process.dump = cms.EDAnalyzer("DumpMkFitGeometry", - level = cms.untracked.int32(1), + level = cms.untracked.int32(1), outputFileName = cms.untracked.string(defaultOutputFileName) ) diff --git a/RecoTracker/MkFit/test/dumpMkFitGeometryPhase2.py b/RecoTracker/MkFit/test/dumpMkFitGeometryPhase2.py new file mode 100644 index 0000000000000..223ba2a04ed19 --- /dev/null +++ b/RecoTracker/MkFit/test/dumpMkFitGeometryPhase2.py @@ -0,0 +1,40 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.Eras.Era_Phase2C17I13M9_cff import Phase2C17I13M9 +# from Configuration.ProcessModifiers.trackingMkFit_cff import trackingMkFit +from Configuration.ProcessModifiers.trackingMkFitCommon_cff import trackingMkFitCommon +trackingMkFit = cms.ModifierChain(trackingMkFitCommon) + +process = cms.Process('DUMP',Phase2C17I13M9,trackingMkFit) + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.Geometry.GeometryExtended2026D88Reco_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.Reconstruction_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +from Configuration.AlCa.autoCond import autoCond +process.GlobalTag.globaltag = autoCond['phase2_realistic'] + +process.MessageLogger.cerr.threshold = "INFO" +process.MessageLogger.cerr.MkFitGeometryESProducer = dict(limit=-1) + +process.source = cms.Source("EmptySource") +process.maxEvents.input = 1 + + +process.add_(cms.ESProducer("MkFitGeometryESProducer")) + +defaultOutputFileName="phase2-trackerinfo.bin" + +# level: 0 - no printout; 1 - print layers, 2 - print modules +# outputFileName: binary dump file; no dump if empty string +process.dump = cms.EDAnalyzer("DumpMkFitGeometry", + level = cms.untracked.int32(1), + outputFileName = cms.untracked.string(defaultOutputFileName) + ) + +print("Requesting MkFit geometry dump into file:", defaultOutputFileName, "\n"); +process.p = cms.Path(process.dump) diff --git a/RecoTracker/MkFit/test/testDumpMkFitGeometry.sh b/RecoTracker/MkFit/test/testDumpMkFitGeometry.sh index a9aea9b387084..b946a3caa57d6 100755 --- a/RecoTracker/MkFit/test/testDumpMkFitGeometry.sh +++ b/RecoTracker/MkFit/test/testDumpMkFitGeometry.sh @@ -8,3 +8,4 @@ if [ "${SCRAM_TEST_NAME}" != "" ] ; then fi (cmsRun ${LOCAL_TEST_DIR}/dumpMkFitGeometry.py) || die "failed to run dumpMkFitGeometry.py" $? +(cmsRun ${LOCAL_TEST_DIR}/dumpMkFitGeometryPhase2.py) || die "failed to run dumpMkFitGeometryPhase2.py" $? diff --git a/RecoTracker/MkFitCMS/interface/LayerNumberConverter.h b/RecoTracker/MkFitCMS/interface/LayerNumberConverter.h index 5528318dec468..b3ee1f4dee8f3 100644 --- a/RecoTracker/MkFitCMS/interface/LayerNumberConverter.h +++ b/RecoTracker/MkFitCMS/interface/LayerNumberConverter.h @@ -1,21 +1,39 @@ #ifndef RecoTracker_MkFitCMS_interface_LayerNumberConverter_h #define RecoTracker_MkFitCMS_interface_LayerNumberConverter_h +#include + namespace mkfit { - enum struct TkLayout { phase0 = 0, phase1 = 1 }; + enum struct TkLayout { phase0 = 0, phase1 = 1, phase2 = 2 }; class LayerNumberConverter { public: LayerNumberConverter(TkLayout layout) : lo_(layout) {} + void reset(TkLayout layout) { lo_ = layout; } unsigned int nLayers() const { if (lo_ == TkLayout::phase0) return 69; if (lo_ == TkLayout::phase1) return 72; + if (lo_ == TkLayout::phase2) + return 60; // 4 + 12 + 2*(12 + 10) = 16 + 22 + 22 = 60 return 10; } + TkLayout getEra() const { return lo_; } int convertLayerNumber(int det, int lay, bool useMatched, int isStereo, bool posZ) const { + if (lo_ == TkLayout::phase2) { + if (det == 1) + return lay - 1; + if (det == 2) + return 16 + lay - 1 + (posZ ? 0 : 22); + if (det == 5) + return 4 + (2 * (lay - 1)) + isStereo; + if (det == 4) + return 16 + 12 + (2 * (lay - 1)) + isStereo + (posZ ? 0 : 22); + throw std::runtime_error("bad subDet"); + } + if (det == 1 || det == 3 || det == 5) { return convertBarrelLayerNumber(det, lay, useMatched, isStereo); } else { @@ -33,6 +51,22 @@ namespace mkfit { } return -1; } + bool doesHaveStereo(int det, int lay) const { + if (lo_ == TkLayout::phase2) { + if (det == 1 || det == 2) + return false; + if (det == 4 || det == 5) + return true; + throw std::runtime_error("bad subDet"); + } + if (det == 3 || det == 5) { + return lay == 1 || lay == 2; + } + if (det == 4 || det == 6) { + return true; + } + return false; + } int convertBarrelLayerNumber(int cmsswdet, int cmsswlay, bool useMatched, int isStereo) const { int lOffset = 0; diff --git a/RecoTracker/MkFitCMS/interface/MkStdSeqs.h b/RecoTracker/MkFitCMS/interface/MkStdSeqs.h index da473bcee36da..9053114d345b8 100644 --- a/RecoTracker/MkFitCMS/interface/MkStdSeqs.h +++ b/RecoTracker/MkFitCMS/interface/MkStdSeqs.h @@ -11,6 +11,8 @@ namespace mkfit { class EventOfHits; class IterationConfig; class TrackerInfo; + class MkJob; + class TrackCand; namespace StdSeq { @@ -24,113 +26,21 @@ namespace mkfit { void cmssw_Map_TrackHitIndices(const EventOfHits &eoh, TrackVec &seeds); void cmssw_ReMap_TrackHitIndices(const EventOfHits &eoh, TrackVec &out_tracks); - int clean_cms_seedtracks_iter(TrackVec *seed_ptr, const IterationConfig &itrcfg, const BeamSpot &bspot); + int clean_cms_seedtracks_iter(TrackVec &seeds, const IterationConfig &itrcfg, const BeamSpot &bspot); - void find_duplicates(TrackVec &tracks); void remove_duplicates(TrackVec &tracks); - void find_duplicates_sharedhits(TrackVec &tracks, const float fraction); - void find_duplicates_sharedhits_pixelseed(TrackVec &tracks, - const float fraction, - const float drth_central, - const float drth_obarrel, - const float drth_forward); + void clean_duplicates(TrackVec &tracks, const IterationConfig &itconf); + void clean_duplicates_sharedhits(TrackVec &tracks, const IterationConfig &itconf); + void clean_duplicates_sharedhits_pixelseed(TrackVec &tracks, const IterationConfig &itconf); - // quality filter for n hits with seed hit "penalty" for strip-based seeds - // this implicitly separates triplets and doublet seeds with glued layers - template - bool qfilter_n_hits(const TRACK &t, int nMinHits) { - int seedHits = t.getNSeedHits(); - int seedReduction = (seedHits <= 5) ? 2 : 3; - return t.nFoundHits() - seedReduction >= nMinHits; - } + // Quality filters used directly (not through IterationConfig) - // simple hit-count quality filter; used with pixel-based seeds template - bool qfilter_n_hits_pixseed(const TRACK &t, int nMinHits) { - return t.nFoundHits() >= nMinHits; - } - - // layer-dependent quality filter - template - bool qfilter_n_layers(const TRACK &t, const BeamSpot &bspot, const TrackerInfo &trk_inf) { - int enhits = t.nHitsByTypeEncoded(trk_inf); - int npixhits = t.nPixelDecoded(enhits); - int enlyrs = t.nLayersByTypeEncoded(trk_inf); - int npixlyrs = t.nPixelDecoded(enlyrs); - int nmatlyrs = t.nTotMatchDecoded(enlyrs); - int llyr = t.getLastFoundHitLyr(); - int lplyr = t.getLastFoundPixelHitLyr(); - float invpt = t.invpT(); - float invptmin = 1.43; // min 1/pT (=1/0.7) for full filter on (npixhits<=3 .or. npixlyrs<=3) - float d0BS = t.d0BeamSpot(bspot.x, bspot.y); - float d0_max = 0.1; // 1 mm - - bool endsInsidePix = (llyr == 2 || llyr == 18 || llyr == 45); - bool lastInsidePix = ((0 <= lplyr && lplyr < 3) || (18 <= lplyr && lplyr < 20) || (45 <= lplyr && lplyr < 47)); - return !(((npixhits <= 3 || npixlyrs <= 3) && endsInsidePix && - (invpt < invptmin || (invpt >= invptmin && std::abs(d0BS) > d0_max))) || - ((npixlyrs <= 3 && nmatlyrs <= 6) && lastInsidePix && llyr != lplyr && std::abs(d0BS) > d0_max)); - } - - /// quality filter tuned for pixelLess iteration during forward search - template - bool qfilter_pixelLessFwd(const TRACK &t, const BeamSpot &bspot, const TrackerInfo &tk_info) { - float d0BS = t.d0BeamSpot(bspot.x, bspot.y); - float d0_max = 0.05; // 0.5 mm - - int encoded; - encoded = t.nLayersByTypeEncoded(tk_info); - int nLyrs = t.nTotMatchDecoded(encoded); - encoded = t.nHitsByTypeEncoded(tk_info); - int nHits = t.nTotMatchDecoded(encoded); - - int seedReduction = (t.getNSeedHits() <= 5) ? 2 : 3; - - float invpt = t.invpT(); - float invptmin = 1.11; // =1/0.9 - - float thetasym = std::abs(t.theta() - Const::PIOver2); - float thetasymmin = 1.11; // -> |eta|=1.45 - - return (((t.nFoundHits() - seedReduction >= 4 && invpt < invptmin) || - (t.nFoundHits() - seedReduction >= 3 && invpt > invptmin && thetasym <= thetasymmin) || - (t.nFoundHits() - seedReduction >= 4 && invpt > invptmin && thetasym > thetasymmin)) && - !((nLyrs <= 4 || nHits <= 4) && std::abs(d0BS) > d0_max && invpt < invptmin)); - } - - /// quality filter tuned for pixelLess iteration during backward search - template - bool qfilter_pixelLessBkwd(const TRACK &t, const BeamSpot &bspot, const TrackerInfo &tk_info) { - float d0BS = t.d0BeamSpot(bspot.x, bspot.y); - float d0_max = 0.1; // 1 mm - - int encoded; - encoded = t.nLayersByTypeEncoded(tk_info); - int nLyrs = t.nTotMatchDecoded(encoded); - encoded = t.nHitsByTypeEncoded(tk_info); - int nHits = t.nTotMatchDecoded(encoded); - - float invpt = t.invpT(); - float invptmin = 1.11; // =1/0.9 - - float thetasym = std::abs(t.theta() - Const::PIOver2); - float thetasymmin_l = 0.80; // -> |eta|=0.9 - float thetasymmin_h = 1.11; // -> |eta|=1.45 - - return !( - ((nLyrs <= 3 || nHits <= 3)) || - ((nLyrs <= 4 || nHits <= 4) && (invpt < invptmin || (thetasym > thetasymmin_l && std::abs(d0BS) > d0_max))) || - ((nLyrs <= 5 || nHits <= 5) && (invpt > invptmin && thetasym > thetasymmin_h && std::abs(d0BS) > d0_max))); - } - - template - bool qfilter_nan_n_silly(const TRACK &t) { + bool qfilter_nan_n_silly(const TRACK &t, const MkJob &) { return !(t.hasNanNSillyValues()); } - void find_and_remove_duplicates(TrackVec &tracks, const IterationConfig &itconf); - } // namespace StdSeq } // namespace mkfit diff --git a/RecoTracker/MkFitCMS/src/MkSeedPartitioners-phase1.cc b/RecoTracker/MkFitCMS/src/MkSeedPartitioners-phase1.cc new file mode 100644 index 0000000000000..feee6642878c1 --- /dev/null +++ b/RecoTracker/MkFitCMS/src/MkSeedPartitioners-phase1.cc @@ -0,0 +1,356 @@ +#include "RecoTracker/MkFitCore/interface/cms_common_macros.h" +#include "RecoTracker/MkFitCore/interface/Track.h" +#include "RecoTracker/MkFitCore/interface/TrackerInfo.h" +#include "RecoTracker/MkFitCore/interface/HitStructures.h" +#include "RecoTracker/MkFitCore/interface/IterationConfig.h" + +namespace { + using namespace mkfit; + + [[maybe_unused]] void partitionSeeds0(const TrackerInfo &trk_info, + const TrackVec &in_seeds, + const EventOfHits &eoh, + IterationSeedPartition &part) { + const size_t size = in_seeds.size(); + + for (size_t i = 0; i < size; ++i) { + const Track &S = in_seeds[i]; + + const bool z_dir_pos = S.pz() > 0; + + const auto &hot = S.getLastHitOnTrack(); + const float eta = eoh[hot.layer].refHit(hot.index).eta(); + + // Region to be defined by propagation / intersection tests + TrackerInfo::EtaRegion reg; + + const LayerInfo &outer_brl = trk_info.outer_barrel_layer(); + + // Define first (mkFit) layer IDs for each strip subdetector. + constexpr int tib1_id = 4; + constexpr int tob1_id = 10; + constexpr int tecp1_id = 27; + constexpr int tecn1_id = 54; + + const LayerInfo &tib1 = trk_info.layer(tib1_id); + const LayerInfo &tob1 = trk_info.layer(tob1_id); + + const LayerInfo &tecp1 = trk_info.layer(tecp1_id); + const LayerInfo &tecn1 = trk_info.layer(tecn1_id); + + const LayerInfo &tec_first = z_dir_pos ? tecp1 : tecn1; + + const float maxR = S.maxReachRadius(); + float z_at_maxr; + + bool can_reach_outer_brl = S.canReachRadius(outer_brl.rout()); + float z_at_outer_brl; + bool misses_first_tec; + if (can_reach_outer_brl) { + z_at_outer_brl = S.zAtR(outer_brl.rout()); + if (z_dir_pos) + misses_first_tec = z_at_outer_brl < tec_first.zmin(); + else + misses_first_tec = z_at_outer_brl > tec_first.zmax(); + } else { + z_at_maxr = S.zAtR(maxR); + if (z_dir_pos) + misses_first_tec = z_at_maxr < tec_first.zmin(); + else + misses_first_tec = z_at_maxr > tec_first.zmax(); + } + + if (misses_first_tec) { + reg = TrackerInfo::Reg_Barrel; + } else { + if ((S.canReachRadius(tib1.rin()) && tib1.is_within_z_limits(S.zAtR(tib1.rin()))) || + (S.canReachRadius(tob1.rin()) && tob1.is_within_z_limits(S.zAtR(tob1.rin())))) { + reg = z_dir_pos ? TrackerInfo::Reg_Transition_Pos : TrackerInfo::Reg_Transition_Neg; + } else { + reg = z_dir_pos ? TrackerInfo::Reg_Endcap_Pos : TrackerInfo::Reg_Endcap_Neg; + } + } + + part.m_region[i] = reg; + if (part.m_phi_eta_foo) + part.m_phi_eta_foo(eoh[hot.layer].refHit(hot.index).phi(), eta); + } + } + + [[maybe_unused]] void partitionSeeds1(const TrackerInfo &trk_info, + const TrackVec &in_seeds, + const EventOfHits &eoh, + IterationSeedPartition &part) { + // Define first (mkFit) layer IDs for each strip subdetector. + constexpr int tib1_id = 4; + constexpr int tob1_id = 10; + constexpr int tidp1_id = 21; + constexpr int tidn1_id = 48; + constexpr int tecp1_id = 27; + constexpr int tecn1_id = 54; + + const LayerInfo &tib1 = trk_info.layer(tib1_id); + const LayerInfo &tob1 = trk_info.layer(tob1_id); + + const LayerInfo &tidp1 = trk_info.layer(tidp1_id); + const LayerInfo &tidn1 = trk_info.layer(tidn1_id); + + const LayerInfo &tecp1 = trk_info.layer(tecp1_id); + const LayerInfo &tecn1 = trk_info.layer(tecn1_id); + + // Merge first two layers to account for mono/stereo coverage. + // TrackerInfo could hold joint limits for sub-detectors. + const auto &L = trk_info; + const float tidp_rin = std::min(L[tidp1_id].rin(), L[tidp1_id + 1].rin()); + const float tidp_rout = std::max(L[tidp1_id].rout(), L[tidp1_id + 1].rout()); + const float tecp_rin = std::min(L[tecp1_id].rin(), L[tecp1_id + 1].rin()); + const float tecp_rout = std::max(L[tecp1_id].rout(), L[tecp1_id + 1].rout()); + const float tidn_rin = std::min(L[tidn1_id].rin(), L[tidn1_id + 1].rin()); + const float tidn_rout = std::max(L[tidn1_id].rout(), L[tidn1_id + 1].rout()); + const float tecn_rin = std::min(L[tecn1_id].rin(), L[tecn1_id + 1].rin()); + const float tecn_rout = std::max(L[tecn1_id].rout(), L[tecn1_id + 1].rout()); + + // Bias towards more aggressive transition-region assignemnts. + // With current tunning it seems to make things a bit worse. + const float tid_z_extra = 0.0f; // 5.0f; + const float tec_z_extra = 0.0f; // 10.0f; + + const size_t size = in_seeds.size(); + + auto barrel_pos_check = [](const Track &S, float maxR, float rin, float zmax) -> bool { + bool inside = maxR > rin && S.zAtR(rin) < zmax; + return inside; + }; + + auto barrel_neg_check = [](const Track &S, float maxR, float rin, float zmin) -> bool { + bool inside = maxR > rin && S.zAtR(rin) > zmin; + return inside; + }; + + auto endcap_pos_check = [](const Track &S, float maxR, float rout, float rin, float zmin) -> bool { + bool inside = maxR > rout ? S.zAtR(rout) > zmin : (maxR > rin && S.zAtR(maxR) > zmin); + return inside; + }; + + auto endcap_neg_check = [](const Track &S, float maxR, float rout, float rin, float zmax) -> bool { + bool inside = maxR > rout ? S.zAtR(rout) < zmax : (maxR > rin && S.zAtR(maxR) < zmax); + return inside; + }; + + for (size_t i = 0; i < size; ++i) { + const Track &S = in_seeds[i]; + + const auto &hot = S.getLastHitOnTrack(); + const float eta = eoh[hot.layer].refHit(hot.index).eta(); + + // Region to be defined by propagation / intersection tests + TrackerInfo::EtaRegion reg; + + const bool z_dir_pos = S.pz() > 0; + const float maxR = S.maxReachRadius(); + + if (z_dir_pos) { + const bool in_tib = barrel_pos_check(S, maxR, tib1.rin(), tib1.zmax()); + const bool in_tob = barrel_pos_check(S, maxR, tob1.rin(), tob1.zmax()); + + if (!in_tib && !in_tob) { + reg = TrackerInfo::Reg_Endcap_Pos; + } else { + const bool in_tid = endcap_pos_check(S, maxR, tidp_rout, tidp_rin, tidp1.zmin() - tid_z_extra); + const bool in_tec = endcap_pos_check(S, maxR, tecp_rout, tecp_rin, tecp1.zmin() - tec_z_extra); + + if (!in_tid && !in_tec) { + reg = TrackerInfo::Reg_Barrel; + } else { + reg = TrackerInfo::Reg_Transition_Pos; + } + } + } else { + const bool in_tib = barrel_neg_check(S, maxR, tib1.rin(), tib1.zmin()); + const bool in_tob = barrel_neg_check(S, maxR, tob1.rin(), tob1.zmin()); + + if (!in_tib && !in_tob) { + reg = TrackerInfo::Reg_Endcap_Neg; + } else { + const bool in_tid = endcap_neg_check(S, maxR, tidn_rout, tidn_rin, tidn1.zmax() + tid_z_extra); + const bool in_tec = endcap_neg_check(S, maxR, tecn_rout, tecn_rin, tecn1.zmax() + tec_z_extra); + + if (!in_tid && !in_tec) { + reg = TrackerInfo::Reg_Barrel; + } else { + reg = TrackerInfo::Reg_Transition_Neg; + } + } + } + + part.m_region[i] = reg; + if (part.m_phi_eta_foo) + part.m_phi_eta_foo(eoh[hot.layer].refHit(hot.index).phi(), eta); + } + } + + [[maybe_unused]] void partitionSeeds1debug(const TrackerInfo &trk_info, + const TrackVec &in_seeds, + const EventOfHits &eoh, + IterationSeedPartition &part) { + // Define first (mkFit) layer IDs for each strip subdetector. + constexpr int tib1_id = 4; + constexpr int tob1_id = 10; + constexpr int tidp1_id = 21; + constexpr int tidn1_id = 48; + constexpr int tecp1_id = 27; + constexpr int tecn1_id = 54; + + const LayerInfo &tib1 = trk_info.layer(tib1_id); + const LayerInfo &tob1 = trk_info.layer(tob1_id); + + const LayerInfo &tidp1 = trk_info.layer(tidp1_id); + const LayerInfo &tidn1 = trk_info.layer(tidn1_id); + + const LayerInfo &tecp1 = trk_info.layer(tecp1_id); + const LayerInfo &tecn1 = trk_info.layer(tecn1_id); + + // Merge first two layers to account for mono/stereo coverage. + // TrackerInfo could hold joint limits for sub-detectors. + const auto &L = trk_info; + const float tidp_rin = std::min(L[tidp1_id].rin(), L[tidp1_id + 1].rin()); + const float tidp_rout = std::max(L[tidp1_id].rout(), L[tidp1_id + 1].rout()); + const float tecp_rin = std::min(L[tecp1_id].rin(), L[tecp1_id + 1].rin()); + const float tecp_rout = std::max(L[tecp1_id].rout(), L[tecp1_id + 1].rout()); + const float tidn_rin = std::min(L[tidn1_id].rin(), L[tidn1_id + 1].rin()); + const float tidn_rout = std::max(L[tidn1_id].rout(), L[tidn1_id + 1].rout()); + const float tecn_rin = std::min(L[tecn1_id].rin(), L[tecn1_id + 1].rin()); + const float tecn_rout = std::max(L[tecn1_id].rout(), L[tecn1_id + 1].rout()); + + // Bias towards more aggressive transition-region assignemnts. + // With current tunning it seems to make things a bit worse. + const float tid_z_extra = 0.0f; // 5.0f; + const float tec_z_extra = 0.0f; // 10.0f; + + const int size = in_seeds.size(); + + auto barrel_pos_check = [](const Track &S, float maxR, float rin, float zmax, const char *det) -> bool { + bool inside = maxR > rin && S.zAtR(rin) < zmax; + + printf(" in_%s=%d maxR=%7.3f, rin=%7.3f -- ", det, inside, maxR, rin); + if (maxR > rin) { + printf("maxR > rin: S.zAtR(rin) < zmax -- %.3f bool { + bool inside = maxR > rin && S.zAtR(rin) > zmin; + + printf(" in_%s=%d maxR=%7.3f, rin=%7.3f -- ", det, inside, maxR, rin); + if (maxR > rin) { + printf("maxR > rin: S.zAtR(rin) > zmin -- %.3f >? %.3f\n", S.zAtR(rin), zmin); + } else { + printf("maxR < rin: no pie.\n"); + } + + return inside; + }; + + auto endcap_pos_check = [](const Track &S, float maxR, float rout, float rin, float zmin, const char *det) -> bool { + bool inside = maxR > rout ? S.zAtR(rout) > zmin : (maxR > rin && S.zAtR(maxR) > zmin); + + printf(" in_%s=%d maxR=%7.3f, rout=%7.3f, rin=%7.3f -- ", det, inside, maxR, rout, rin); + if (maxR > rout) { + printf("maxR > rout: S.zAtR(rout) > zmin -- %.3f >? %.3f\n", S.zAtR(rout), zmin); + } else if (maxR > rin) { + printf("maxR > rin: S.zAtR(maxR) > zmin) -- %.3f >? %.3f\n", S.zAtR(maxR), zmin); + } else { + printf("maxR < rin: no pie.\n"); + } + + return inside; + }; + + auto endcap_neg_check = [](const Track &S, float maxR, float rout, float rin, float zmax, const char *det) -> bool { + bool inside = maxR > rout ? S.zAtR(rout) < zmax : (maxR > rin && S.zAtR(maxR) < zmax); + + printf(" in_%s=%d maxR=%7.3f, rout=%7.3f, rin=%7.3f -- ", det, inside, maxR, rout, rin); + if (maxR > rout) { + printf("maxR > rout: S.zAtR(rout) < zmax -- %.3f rin) { + printf("maxR > rin: S.zAtR(maxR) < zmax -- %.3f 0; + const float maxR = S.maxReachRadius(); + + printf("partitionSeeds1debug seed index %d, z_dir_pos=%d (pz=%.3f), maxR=%.3f\n", i, z_dir_pos, S.pz(), maxR); + + if (z_dir_pos) { + bool in_tib = barrel_pos_check(S, maxR, tib1.rin(), tib1.zmax(), "TIBp"); + bool in_tob = barrel_pos_check(S, maxR, tob1.rin(), tob1.zmax(), "TOBp"); + + if (!in_tib && !in_tob) { + reg = TrackerInfo::Reg_Endcap_Pos; + printf(" --> region = %d, endcap pos\n", reg); + } else { + bool in_tid = endcap_pos_check(S, maxR, tidp_rout, tidp_rin, tidp1.zmin() - tid_z_extra, "TIDp"); + bool in_tec = endcap_pos_check(S, maxR, tecp_rout, tecp_rin, tecp1.zmin() - tec_z_extra, "TECp"); + + if (!in_tid && !in_tec) { + reg = TrackerInfo::Reg_Barrel; + printf(" --> region = %d, barrel\n", reg); + } else { + reg = TrackerInfo::Reg_Transition_Pos; + printf(" --> region = %d, transition pos\n", reg); + } + } + } else { + bool in_tib = barrel_neg_check(S, maxR, tib1.rin(), tib1.zmin(), "TIBn"); + bool in_tob = barrel_neg_check(S, maxR, tob1.rin(), tob1.zmin(), "TOBn"); + + if (!in_tib && !in_tob) { + reg = TrackerInfo::Reg_Endcap_Neg; + printf(" --> region = %d, endcap neg\n", reg); + } else { + bool in_tid = endcap_neg_check(S, maxR, tidn_rout, tidn_rin, tidn1.zmax() + tid_z_extra, "TIDn"); + bool in_tec = endcap_neg_check(S, maxR, tecn_rout, tecn_rin, tecn1.zmax() + tec_z_extra, "TECn"); + + if (!in_tid && !in_tec) { + reg = TrackerInfo::Reg_Barrel; + printf(" --> region = %d, barrel\n", reg); + } else { + reg = TrackerInfo::Reg_Transition_Neg; + printf(" --> region = %d, transition neg\n", reg); + } + } + } + + part.m_region[i] = reg; + if (part.m_phi_eta_foo) + part.m_phi_eta_foo(eoh[hot.layer].refHit(hot.index).phi(), eta); + } + } + + CMS_SA_ALLOW struct register_seed_partitioners { + register_seed_partitioners() { + IterationConfig::register_seed_partitioner("phase1:0", partitionSeeds0); + IterationConfig::register_seed_partitioner("phase1:1", partitionSeeds1); + IterationConfig::register_seed_partitioner("phase1:1:debug", partitionSeeds1debug); + } + } rsp_instance; +} // namespace diff --git a/RecoTracker/MkFitCMS/src/MkSeedPartitioners-phase2.cc b/RecoTracker/MkFitCMS/src/MkSeedPartitioners-phase2.cc new file mode 100644 index 0000000000000..a10a7a1c60f87 --- /dev/null +++ b/RecoTracker/MkFitCMS/src/MkSeedPartitioners-phase2.cc @@ -0,0 +1,261 @@ +#include "RecoTracker/MkFitCore/interface/cms_common_macros.h" +#include "RecoTracker/MkFitCore/interface/Track.h" +#include "RecoTracker/MkFitCore/interface/TrackerInfo.h" +#include "RecoTracker/MkFitCore/interface/HitStructures.h" +#include "RecoTracker/MkFitCore/interface/IterationConfig.h" + +namespace { + using namespace mkfit; + + // named constants for useful layers (l/u for lower/upper) + constexpr int tecp1l_id = 28; + constexpr int tecp1u_id = 29; + constexpr int tecp2l_id = 30; + constexpr int tecp2u_id = 31; + constexpr int tecn1l_id = 50; + constexpr int tecn1u_id = 51; + constexpr int tecn2l_id = 52; + constexpr int tecntu_id = 53; + + [[maybe_unused]] void partitionSeeds1(const TrackerInfo &trk_info, + const TrackVec &in_seeds, + const EventOfHits &eoh, + IterationSeedPartition &part) { + // Seeds are placed into eta regions and sorted on region + eta. + + // Merge mono and stereo limits for relevant layers / parameters. + // TrackerInfo could hold joint limits for sub-detectors. + const auto &L = trk_info; + const float tecp1_rin = std::min(L[tecp1l_id].rin(), L[tecp1u_id].rin()); + const float tecp1_rout = std::max(L[tecp1l_id].rout(), L[tecp1u_id].rout()); + const float tecp1_zmin = std::min(L[tecp1l_id].zmin(), L[tecp1u_id].zmin()); + + const float tecp2_rin = std::min(L[tecp2l_id].rin(), L[tecp2u_id].rin()); + const float tecp2_zmax = std::max(L[tecp2l_id].zmax(), L[tecp2u_id].zmax()); + + const float tecn1_rin = std::min(L[tecn1l_id].rin(), L[tecn1u_id].rin()); + const float tecn1_rout = std::max(L[tecn1l_id].rout(), L[tecn1u_id].rout()); + const float tecn1_zmax = std::max(L[tecn1l_id].zmax(), L[tecn1u_id].zmax()); + + const float tecn2_rin = std::min(L[tecn2l_id].rin(), L[tecntu_id].rin()); + const float tecn2_zmin = std::min(L[tecn2l_id].zmin(), L[tecntu_id].zmin()); + + const float tec_z_extra = 0.0f; // 10.0f; + + const int size = in_seeds.size(); + + auto barrel_pos_check = [](const Track &S, float maxR, float rin, float zmax) -> bool { + bool inside = maxR > rin && S.zAtR(rin) < zmax; + return inside; + }; + + auto barrel_neg_check = [](const Track &S, float maxR, float rin, float zmin) -> bool { + bool inside = maxR > rin && S.zAtR(rin) > zmin; + return inside; + }; + + auto endcap_pos_check = [](const Track &S, float maxR, float rout, float rin, float zmin) -> bool { + bool inside = maxR > rout ? S.zAtR(rout) > zmin : (maxR > rin && S.zAtR(maxR) > zmin); + return inside; + }; + + auto endcap_neg_check = [](const Track &S, float maxR, float rout, float rin, float zmax) -> bool { + bool inside = maxR > rout ? S.zAtR(rout) < zmax : (maxR > rin && S.zAtR(maxR) < zmax); + return inside; + }; + + for (int i = 0; i < size; ++i) { + const Track &S = in_seeds[i]; + + HitOnTrack hot = S.getLastHitOnTrack(); + float eta = eoh[hot.layer].refHit(hot.index).eta(); + + // Region to be defined by propagation / intersection tests + TrackerInfo::EtaRegion reg; + + const bool z_dir_pos = S.pz() > 0; + const float maxR = S.maxReachRadius(); + + if (z_dir_pos) { + bool in_tec_as_brl = barrel_pos_check(S, maxR, tecp2_rin, tecp2_zmax); + + if (!in_tec_as_brl) { + reg = TrackerInfo::Reg_Endcap_Pos; + } else { + bool in_tec = endcap_pos_check(S, maxR, tecp1_rout, tecp1_rin, tecp1_zmin - tec_z_extra); + + if (!in_tec) { + reg = TrackerInfo::Reg_Barrel; + } else { + reg = TrackerInfo::Reg_Transition_Pos; + } + } + } else { + bool in_tec_as_brl = barrel_neg_check(S, maxR, tecn2_rin, tecn2_zmin); + + if (!in_tec_as_brl) { + reg = TrackerInfo::Reg_Endcap_Neg; + } else { + bool in_tec = endcap_neg_check(S, maxR, tecn1_rout, tecn1_rin, tecn1_zmax + tec_z_extra); + + if (!in_tec) { + reg = TrackerInfo::Reg_Barrel; + } else { + reg = TrackerInfo::Reg_Transition_Neg; + } + } + } + + part.m_region[i] = reg; + if (part.m_phi_eta_foo) + part.m_phi_eta_foo(eoh[hot.layer].refHit(hot.index).phi(), eta); + } + } + + [[maybe_unused]] void partitionSeeds1debug(const TrackerInfo &trk_info, + const TrackVec &in_seeds, + const EventOfHits &eoh, + IterationSeedPartition &part) { + // Seeds are placed into eta regions and sorted on region + eta. + + // Merge mono and stereo limits for relevant layers / parameters. + // TrackerInfo could hold joint limits for sub-detectors. + const auto &L = trk_info; + const float tecp1_rin = std::min(L[tecp1l_id].rin(), L[tecp1u_id].rin()); + const float tecp1_rout = std::max(L[tecp1l_id].rout(), L[tecp1u_id].rout()); + const float tecp1_zmin = std::min(L[tecp1l_id].zmin(), L[tecp1u_id].zmin()); + + const float tecp2_rin = std::min(L[tecp2l_id].rin(), L[tecp2u_id].rin()); + const float tecp2_zmax = std::max(L[tecp2l_id].zmax(), L[tecp2u_id].zmax()); + + const float tecn1_rin = std::min(L[tecn1l_id].rin(), L[tecn1u_id].rin()); + const float tecn1_rout = std::max(L[tecn1l_id].rout(), L[tecn1u_id].rout()); + const float tecn1_zmax = std::max(L[tecn1l_id].zmax(), L[tecn1u_id].zmax()); + + const float tecn2_rin = std::min(L[tecn2l_id].rin(), L[tecntu_id].rin()); + const float tecn2_zmin = std::min(L[tecn2l_id].zmin(), L[tecntu_id].zmin()); + + const float tec_z_extra = 0.0f; // 10.0f; + + const int size = in_seeds.size(); + + auto barrel_pos_check = [](const Track &S, float maxR, float rin, float zmax, const char *det) -> bool { + bool inside = maxR > rin && S.zAtR(rin) < zmax; + + printf(" in_%s=%d maxR=%7.3f, rin=%7.3f -- ", det, inside, maxR, rin); + if (maxR > rin) { + printf("maxR > rin: S.zAtR(rin) < zmax -- %.3f bool { + bool inside = maxR > rin && S.zAtR(rin) > zmin; + + printf(" in_%s=%d maxR=%7.3f, rin=%7.3f -- ", det, inside, maxR, rin); + if (maxR > rin) { + printf("maxR > rin: S.zAtR(rin) > zmin -- %.3f >? %.3f\n", S.zAtR(rin), zmin); + } else { + printf("maxR < rin: no pie.\n"); + } + + return inside; + }; + + auto endcap_pos_check = [](const Track &S, float maxR, float rout, float rin, float zmin, const char *det) -> bool { + bool inside = maxR > rout ? S.zAtR(rout) > zmin : (maxR > rin && S.zAtR(maxR) > zmin); + + printf(" in_%s=%d maxR=%7.3f, rout=%7.3f, rin=%7.3f -- ", det, inside, maxR, rout, rin); + if (maxR > rout) { + printf("maxR > rout: S.zAtR(rout) > zmin -- %.3f >? %.3f\n", S.zAtR(rout), zmin); + } else if (maxR > rin) { + printf("maxR > rin: S.zAtR(maxR) > zmin) -- %.3f >? %.3f\n", S.zAtR(maxR), zmin); + } else { + printf("maxR < rin: no pie.\n"); + } + + return inside; + }; + + auto endcap_neg_check = [](const Track &S, float maxR, float rout, float rin, float zmax, const char *det) -> bool { + bool inside = maxR > rout ? S.zAtR(rout) < zmax : (maxR > rin && S.zAtR(maxR) < zmax); + + printf(" in_%s=%d maxR=%7.3f, rout=%7.3f, rin=%7.3f -- ", det, inside, maxR, rout, rin); + if (maxR > rout) { + printf("maxR > rout: S.zAtR(rout) < zmax -- %.3f rin) { + printf("maxR > rin: S.zAtR(maxR) < zmax -- %.3f 0; + const float maxR = S.maxReachRadius(); + + printf("partitionSeeds1debug seed index %d, z_dir_pos=%d (pz=%.3f), maxR=%.3f\n", i, z_dir_pos, S.pz(), maxR); + + if (z_dir_pos) { + bool in_tec_as_brl = barrel_pos_check(S, maxR, tecp2_rin, tecp2_zmax, "TECasBarrelp"); + + if (!in_tec_as_brl) { + reg = TrackerInfo::Reg_Endcap_Pos; + printf(" --> region = %d, endcap pos\n", reg); + } else { + bool in_tec = endcap_pos_check(S, maxR, tecp1_rout, tecp1_rin, tecp1_zmin - tec_z_extra, "TECp"); + + if (!in_tec) { + reg = TrackerInfo::Reg_Barrel; + printf(" --> region = %d, barrel\n", reg); + } else { + reg = TrackerInfo::Reg_Transition_Pos; + printf(" --> region = %d, transition pos\n", reg); + } + } + } else { + bool in_tec_as_brl = barrel_neg_check(S, maxR, tecn2_rin, tecn2_zmin, "TECasBarreln"); + + if (!in_tec_as_brl) { + reg = TrackerInfo::Reg_Endcap_Neg; + printf(" --> region = %d, endcap neg\n", reg); + } else { + bool in_tec = endcap_neg_check(S, maxR, tecn1_rout, tecn1_rin, tecn1_zmax + tec_z_extra, "TECn"); + + if (!in_tec) { + reg = TrackerInfo::Reg_Barrel; + printf(" --> region = %d, barrel\n", reg); + } else { + reg = TrackerInfo::Reg_Transition_Neg; + printf(" --> region = %d, transition neg\n", reg); + } + } + } + + part.m_region[i] = reg; + if (part.m_phi_eta_foo) + part.m_phi_eta_foo(eoh[hot.layer].refHit(hot.index).phi(), eta); + } + } + + CMS_SA_ALLOW struct register_seed_partitioners { + register_seed_partitioners() { + IterationConfig::register_seed_partitioner("phase2:1", partitionSeeds1); + IterationConfig::register_seed_partitioner("phase2:1:debug", partitionSeeds1debug); + } + } rsp_instance; +} // namespace diff --git a/RecoTracker/MkFitCMS/src/MkStdSeqs.cc b/RecoTracker/MkFitCMS/src/MkStdSeqs.cc index 2c1ff9d008f0a..e1c2f5de653e9 100644 --- a/RecoTracker/MkFitCMS/src/MkStdSeqs.cc +++ b/RecoTracker/MkFitCMS/src/MkStdSeqs.cc @@ -1,7 +1,10 @@ +#include "RecoTracker/MkFitCore/interface/cms_common_macros.h" #include "RecoTracker/MkFitCMS/interface/MkStdSeqs.h" #include "RecoTracker/MkFitCore/interface/HitStructures.h" #include "RecoTracker/MkFitCore/interface/IterationConfig.h" +#include "RecoTracker/MkFitCore/interface/MkJob.h" +#include "RecoTracker/MkFitCore/interface/TrackStructures.h" #include "RecoTracker/MkFitCore/interface/binnor.h" @@ -81,22 +84,22 @@ namespace mkfit { //========================================================================= // Seed cleaning (multi-iter) //========================================================================= - int clean_cms_seedtracks_iter(TrackVec *seed_ptr, const IterationConfig &itrcfg, const BeamSpot &bspot) { + int clean_cms_seedtracks_iter(TrackVec &seeds, const IterationConfig &itrcfg, const BeamSpot &bspot) { using Algo = TrackBase::TrackAlgorithm; const float etamax_brl = Config::c_etamax_brl; const float dpt_common = Config::c_dpt_common; - const float dzmax_bh = itrcfg.m_params.c_dzmax_bh; - const float drmax_bh = itrcfg.m_params.c_drmax_bh; - const float dzmax_eh = itrcfg.m_params.c_dzmax_eh; - const float drmax_eh = itrcfg.m_params.c_drmax_eh; - const float dzmax_bl = itrcfg.m_params.c_dzmax_bl; - const float drmax_bl = itrcfg.m_params.c_drmax_bl; - const float dzmax_el = itrcfg.m_params.c_dzmax_el; - const float drmax_el = itrcfg.m_params.c_drmax_el; + const float dzmax_bh = itrcfg.sc_dzmax_bh; + const float drmax_bh = itrcfg.sc_drmax_bh; + const float dzmax_eh = itrcfg.sc_dzmax_eh; + const float drmax_eh = itrcfg.sc_drmax_eh; + const float dzmax_bl = itrcfg.sc_dzmax_bl; + const float drmax_bl = itrcfg.sc_drmax_bl; + const float dzmax_el = itrcfg.sc_dzmax_el; + const float drmax_el = itrcfg.sc_drmax_el; - const float ptmin_hpt = itrcfg.m_params.c_ptthr_hpt; + const float ptmin_hpt = itrcfg.sc_ptthr_hpt; const float dzmax2_inv_bh = 1.f / (dzmax_bh * dzmax_bh); const float drmax2_inv_bh = 1.f / (drmax_bh * drmax_bh); @@ -114,9 +117,8 @@ namespace mkfit { const float ptmax_merge_lowPtQuad = 0.2; const float etamin_merge_lowPtQuad = 1.5; - if (seed_ptr == nullptr) + if (seeds.empty()) return 0; - TrackVec &seeds = *seed_ptr; const int ns = seeds.size(); #ifdef DEBUG @@ -335,11 +337,24 @@ namespace mkfit { return seeds.size(); } + namespace { + CMS_SA_ALLOW struct register_seed_cleaners { + register_seed_cleaners() { + IterationConfig::register_seed_cleaner("phase1:default", clean_cms_seedtracks_iter); + } + } rsc_instance; + } // namespace + //========================================================================= // Duplicate cleaning //========================================================================= - void find_duplicates(TrackVec &tracks) { + void remove_duplicates(TrackVec &tracks) { + tracks.erase(std::remove_if(tracks.begin(), tracks.end(), [](auto track) { return track.getDuplicateValue(); }), + tracks.end()); + } + + void clean_duplicates(TrackVec &tracks, const IterationConfig &) { const auto ntracks = tracks.size(); float eta1, phi1, pt1, deta, dphi, dr2; @@ -348,10 +363,6 @@ namespace mkfit { } for (auto itrack = 0U; itrack < ntracks - 1; itrack++) { auto &track = tracks[itrack]; - using Algo = TrackBase::TrackAlgorithm; - auto const algo = track.algorithm(); - if (algo == Algo::pixelLessStep || algo == Algo::tobTecStep) - continue; eta1 = track.momEta(); phi1 = track.momPhi(); pt1 = track.pT(); @@ -421,18 +432,16 @@ namespace mkfit { } //end of else } //end of loop over track2 } //end of loop over track1 - } - void remove_duplicates(TrackVec &tracks) { - tracks.erase(std::remove_if(tracks.begin(), tracks.end(), [](auto track) { return track.getDuplicateValue(); }), - tracks.end()); + remove_duplicates(tracks); } //========================================================================= // SHARED HITS DUPLICATE CLEANING //========================================================================= - void find_duplicates_sharedhits(TrackVec &tracks, const float fraction) { + void clean_duplicates_sharedhits(TrackVec &tracks, const IterationConfig &itconf) { + const float fraction = itconf.dc_fracSharedHits; const auto ntracks = tracks.size(); std::vector ctheta(ntracks); @@ -492,15 +501,14 @@ namespace mkfit { } // end sharing hits loop } // end trk loop - tracks.erase(std::remove_if(tracks.begin(), tracks.end(), [](auto track) { return track.getDuplicateValue(); }), - tracks.end()); + remove_duplicates(tracks); } - void find_duplicates_sharedhits_pixelseed(TrackVec &tracks, - const float fraction, - const float drth_central, - const float drth_obarrel, - const float drth_forward) { + void clean_duplicates_sharedhits_pixelseed(TrackVec &tracks, const IterationConfig &itconf) { + const float fraction = itconf.dc_fracSharedHits; + const float drth_central = itconf.dc_drth_central; + const float drth_obarrel = itconf.dc_drth_obarrel; + const float drth_forward = itconf.dc_drth_forward; const auto ntracks = tracks.size(); std::vector ctheta(ntracks); @@ -586,39 +594,144 @@ namespace mkfit { } } //end loop one over tracks - //removal here - tracks.erase(std::remove_if(tracks.begin(), tracks.end(), [](auto track) { return track.getDuplicateValue(); }), - tracks.end()); + remove_duplicates(tracks); } + namespace { + CMS_SA_ALLOW struct register_duplicate_cleaners { + register_duplicate_cleaners() { + IterationConfig::register_duplicate_cleaner("phase1:clean_duplicates", clean_duplicates); + IterationConfig::register_duplicate_cleaner("phase1:clean_duplicates_sharedhits", + clean_duplicates_sharedhits); + IterationConfig::register_duplicate_cleaner("phase1:clean_duplicates_sharedhits_pixelseed", + clean_duplicates_sharedhits_pixelseed); + } + } rdc_instance; + } // namespace + //========================================================================= - // + // Quality filters //========================================================================= - void find_and_remove_duplicates(TrackVec &tracks, const IterationConfig &itconf) { -#ifdef DEBUG - std::cout << " find_and_remove_duplicates: input track size " << tracks.size() << std::endl; -#endif - if (itconf.m_requires_quality_filter && !(itconf.m_requires_dupclean_tight)) { - find_duplicates_sharedhits(tracks, itconf.m_params.fracSharedHits); - } else if (itconf.m_requires_dupclean_tight) { - find_duplicates_sharedhits_pixelseed(tracks, - itconf.m_params.fracSharedHits, - itconf.m_params.drth_central, - itconf.m_params.drth_obarrel, - itconf.m_params.drth_forward); - } else { - find_duplicates(tracks); - remove_duplicates(tracks); - } + // quality filter for n hits with seed hit "penalty" for strip-based seeds + // this implicitly separates triplets and doublet seeds with glued layers + template + bool qfilter_n_hits(const TRACK &t, const MkJob &j) { + int seedHits = t.getNSeedHits(); + int seedReduction = (seedHits <= 5) ? 2 : 3; + return t.nFoundHits() - seedReduction >= j.params_cur().minHitsQF; + } -#ifdef DEBUG - std::cout << " find_and_remove_duplicates: output track size " << tracks.size() << std::endl; - for (auto const &tk : tracks) { - std::cout << tk.parameters() << std::endl; - } -#endif + // simple hit-count quality filter; used with pixel-based seeds + template + bool qfilter_n_hits_pixseed(const TRACK &t, const MkJob &j) { + return t.nFoundHits() >= j.params_cur().minHitsQF; } + // layer-dependent quality filter + // includes ad hoc tuning for phase-1 + template + bool qfilter_n_layers(const TRACK &t, const MkJob &j) { + const BeamSpot &bspot = j.m_beam_spot; + const TrackerInfo &trk_inf = j.m_trk_info; + int enhits = t.nHitsByTypeEncoded(trk_inf); + int npixhits = t.nPixelDecoded(enhits); + int enlyrs = t.nLayersByTypeEncoded(trk_inf); + int npixlyrs = t.nPixelDecoded(enlyrs); + int nmatlyrs = t.nTotMatchDecoded(enlyrs); + int llyr = t.getLastFoundHitLyr(); + int lplyr = t.getLastFoundPixelHitLyr(); + float invpt = t.invpT(); + + // based on fr and eff vs pt (convert to native invpt) + float invptmin = 1.43; // min 1/pT (=1/0.7) for full filter on (npixhits<=3 .or. npixlyrs<=3) + float d0BS = t.d0BeamSpot(bspot.x, bspot.y); + float d0_max = 0.1; // 1 mm, max for somewhat prompt + + // next-to-outermost pixel layers (almost): BPIX3 or FPIX1 + bool endsInsidePix = (llyr == 2 || llyr == 18 || llyr == 45); + // not last pixel layers: BPIX[123] or FPIX[12] + bool lastInsidePix = ((0 <= lplyr && lplyr < 3) || (18 <= lplyr && lplyr < 20) || (45 <= lplyr && lplyr < 47)); + // reject short tracks missing last pixel layer except for prompt-looking + return !(((npixhits <= 3 || npixlyrs <= 3) && endsInsidePix && + (invpt < invptmin || (invpt >= invptmin && std::abs(d0BS) > d0_max))) || + ((npixlyrs <= 3 && nmatlyrs <= 6) && lastInsidePix && llyr != lplyr && std::abs(d0BS) > d0_max)); + } + + /// quality filter tuned for pixelLess iteration during forward search + // includes ad hoc tuning for phase-1 + template + bool qfilter_pixelLessFwd(const TRACK &t, const MkJob &j) { + const BeamSpot &bspot = j.m_beam_spot; + const TrackerInfo &tk_info = j.m_trk_info; + float d0BS = t.d0BeamSpot(bspot.x, bspot.y); + float d0_max = 0.05; // 0.5 mm, max for somewhat prompt + + int encoded; + encoded = t.nLayersByTypeEncoded(tk_info); + int nLyrs = t.nTotMatchDecoded(encoded); + encoded = t.nHitsByTypeEncoded(tk_info); + int nHits = t.nTotMatchDecoded(encoded); + + // to subtract stereo seed layers to count just r-phi seed layers (better pt err) + int seedReduction = (t.getNSeedHits() <= 5) ? 2 : 3; + + // based on fr and eff vs pt and eta (convert to native invpt and theta) + float invpt = t.invpT(); + float invptmin = 1.11; // =1/0.9 + + float thetasym = std::abs(t.theta() - Const::PIOver2); + float thetasymmin = 1.11; // -> |eta|=1.45 + + // accept longer tracks, reject too short and displaced + return (((t.nFoundHits() - seedReduction >= 4 && invpt < invptmin) || + (t.nFoundHits() - seedReduction >= 3 && invpt > invptmin && thetasym <= thetasymmin) || + (t.nFoundHits() - seedReduction >= 4 && invpt > invptmin && thetasym > thetasymmin)) && + !((nLyrs <= 4 || nHits <= 4) && std::abs(d0BS) > d0_max && invpt < invptmin)); + } + + /// quality filter tuned for pixelLess iteration during backward search + // includes ad hoc tuning for phase-1 + template + bool qfilter_pixelLessBkwd(const TRACK &t, const MkJob &j) { + const BeamSpot &bspot = j.m_beam_spot; + const TrackerInfo &tk_info = j.m_trk_info; + float d0BS = t.d0BeamSpot(bspot.x, bspot.y); + float d0_max = 0.1; // 1 mm + + int encoded; + encoded = t.nLayersByTypeEncoded(tk_info); + int nLyrs = t.nTotMatchDecoded(encoded); + encoded = t.nHitsByTypeEncoded(tk_info); + int nHits = t.nTotMatchDecoded(encoded); + + // based on fr and eff vs pt and eta (convert to native invpt and theta) + float invpt = t.invpT(); + float invptmin = 1.11; // =1/0.9 + + float thetasym = std::abs(t.theta() - Const::PIOver2); + float thetasymmin_l = 0.80; // -> |eta|=0.9 + float thetasymmin_h = 1.11; // -> |eta|=1.45 + + // reject too short or too displaced tracks + return !( + ((nLyrs <= 3 || nHits <= 3)) || + ((nLyrs <= 4 || nHits <= 4) && (invpt < invptmin || (thetasym > thetasymmin_l && std::abs(d0BS) > d0_max))) || + ((nLyrs <= 5 || nHits <= 5) && (invpt > invptmin && thetasym > thetasymmin_h && std::abs(d0BS) > d0_max))); + } + + namespace { + CMS_SA_ALLOW struct register_quality_filters { + register_quality_filters() { + IterationConfig::register_candidate_filter("phase1:qfilter_n_hits", qfilter_n_hits); + IterationConfig::register_candidate_filter("phase1:qfilter_n_hits_pixseed", + qfilter_n_hits_pixseed); + IterationConfig::register_candidate_filter("phase1:qfilter_n_layers", qfilter_n_layers); + IterationConfig::register_candidate_filter("phase1:qfilter_pixelLessFwd", qfilter_pixelLessFwd); + IterationConfig::register_candidate_filter("phase1:qfilter_pixelLessBkwd", qfilter_pixelLessBkwd); + } + } rqf_instance; + } // namespace + } // namespace StdSeq } // namespace mkfit diff --git a/RecoTracker/MkFitCMS/src/runFunctions.cc b/RecoTracker/MkFitCMS/src/runFunctions.cc index 71ecd647d7dd1..3f5106f773065 100644 --- a/RecoTracker/MkFitCMS/src/runFunctions.cc +++ b/RecoTracker/MkFitCMS/src/runFunctions.cc @@ -38,15 +38,15 @@ namespace mkfit { bool do_remove_duplicates) { IterationMaskIfcCmssw it_mask_ifc(trackerInfo, hit_masks); - MkJob job({trackerInfo, itconf, eoh, &it_mask_ifc}); + MkJob job({trackerInfo, itconf, eoh, eoh.refBeamSpot(), &it_mask_ifc}); builder.begin_event(&job, nullptr, __func__); - // Seed cleaning not done on pixelLess / tobTec iters - do_seed_clean = do_seed_clean && itconf.m_requires_dupclean_tight; + // Seed cleaning not done on all iterations. + do_seed_clean = do_seed_clean && itconf.m_seed_cleaner; if (do_seed_clean) - StdSeq::clean_cms_seedtracks_iter(&seeds, itconf, eoh.refBeamSpot()); + itconf.m_seed_cleaner(seeds, itconf, eoh.refBeamSpot()); // Check nans in seeds -- this should not be needed when Slava fixes // the track parameter coordinate transformation. @@ -61,19 +61,12 @@ namespace mkfit { builder.findTracksCloneEngine(); - using Algo = TrackBase::TrackAlgorithm; - if (itconf.m_requires_quality_filter && Algo(itconf.m_track_algorithm) != Algo::detachedTripletStep) { - if (Algo(itconf.m_track_algorithm) == Algo::pixelPairStep) { - builder.filter_comb_cands([&](const TrackCand &t) { return StdSeq::qfilter_n_hits_pixseed(t, 3); }); - } else if (Algo(itconf.m_track_algorithm) == Algo::pixelLessStep) { - builder.filter_comb_cands( - [&](const TrackCand &t) { return StdSeq::qfilter_pixelLessFwd(t, eoh.refBeamSpot(), trackerInfo); }); - } else { - builder.filter_comb_cands( - [&](const TrackCand &t) { return StdSeq::qfilter_n_hits(t, itconf.m_params.minHitsQF); }); - } + if (itconf.m_pre_bkfit_filter) { + builder.filter_comb_cands(itconf.m_pre_bkfit_filter); } + job.switch_to_backward(); + if (do_backward_fit) { if (itconf.m_backward_search) { builder.compactifyHitStorageForBestCand(itconf.m_backward_drop_seed_hits, itconf.m_backward_fit_min_hits); @@ -87,24 +80,17 @@ namespace mkfit { builder.endBkwSearch(); } - if (itconf.m_requires_quality_filter && (Algo(itconf.m_track_algorithm) == Algo::detachedTripletStep || - Algo(itconf.m_track_algorithm) == Algo::pixelLessStep)) { - if (Algo(itconf.m_track_algorithm) == Algo::detachedTripletStep) { - builder.filter_comb_cands( - [&](const TrackCand &t) { return StdSeq::qfilter_n_layers(t, eoh.refBeamSpot(), trackerInfo); }); - } else if (Algo(itconf.m_track_algorithm) == Algo::pixelLessStep) { - builder.filter_comb_cands( - [&](const TrackCand &t) { return StdSeq::qfilter_pixelLessBkwd(t, eoh.refBeamSpot(), trackerInfo); }); - } + if (itconf.m_post_bkfit_filter) { + builder.filter_comb_cands(itconf.m_post_bkfit_filter); } } - builder.filter_comb_cands([&](const TrackCand &t) { return StdSeq::qfilter_nan_n_silly(t); }); + builder.filter_comb_cands(StdSeq::qfilter_nan_n_silly); builder.export_best_comb_cands(out_tracks, true); - if (do_remove_duplicates) { - StdSeq::find_and_remove_duplicates(out_tracks, itconf); + if (do_remove_duplicates && itconf.m_duplicate_cleaner) { + itconf.m_duplicate_cleaner(out_tracks, itconf); } builder.end_event(); diff --git a/RecoTracker/MkFitCMS/standalone/Geoms/CMS-2017.cc b/RecoTracker/MkFitCMS/standalone/Geoms/CMS-phase1.cc similarity index 58% rename from RecoTracker/MkFitCMS/standalone/Geoms/CMS-2017.cc rename to RecoTracker/MkFitCMS/standalone/Geoms/CMS-phase1.cc index 4dc7529961699..b96178a23006f 100644 --- a/RecoTracker/MkFitCMS/standalone/Geoms/CMS-2017.cc +++ b/RecoTracker/MkFitCMS/standalone/Geoms/CMS-phase1.cc @@ -1,7 +1,7 @@ //------------------- -// CMS 2017 geometry +// CMS phase1 geometry //------------------- // Redirect to external geometry creator function. -#include "mkfit-geom-cms-2017/CMS-2017.cc" +#include "mkfit-geom-cms-phase1/CMS-phase1.cc" diff --git a/RecoTracker/MkFitCMS/standalone/Geoms/CMS-phase2.cc b/RecoTracker/MkFitCMS/standalone/Geoms/CMS-phase2.cc new file mode 100644 index 0000000000000..81be16df899f9 --- /dev/null +++ b/RecoTracker/MkFitCMS/standalone/Geoms/CMS-phase2.cc @@ -0,0 +1,7 @@ +//------------------- +// CMS phase2 geometry +//------------------- + +// Redirect to external geometry creator function. + +#include "mkfit-geom-cms-phase2/CMS-phase2.cc" diff --git a/RecoTracker/MkFitCMS/standalone/Geoms/Makefile b/RecoTracker/MkFitCMS/standalone/Geoms/Makefile index 436a7fe7419e1..28d92c8c5ac3e 100644 --- a/RecoTracker/MkFitCMS/standalone/Geoms/Makefile +++ b/RecoTracker/MkFitCMS/standalone/Geoms/Makefile @@ -4,7 +4,7 @@ CPPFLAGS := -I${SRCDIR} -I../mkFit-external ${CPPFLAGS} .PHONY: all clean distclean echo -SRCS := ${SACMS}/Geoms/CMS-2017.cc +SRCS := ${SACMS}/Geoms/CMS-phase1.cc ${SACMS}/Geoms/CMS-phase2.cc SRCB := $(notdir ${SRCS}) DEPS := $(SRCB:.cc=.d) OBJS := $(SRCB:.cc=.o) @@ -12,11 +12,12 @@ OBJS := $(SRCB:.cc=.o) TGTS := $(basename ${OBJS}) TGTS := $(addprefix ../, $(addsuffix .so, ${TGTS})) -GEO_2017_BIN = ../CMS-2017.bin +GEO_phase1_BIN = ../CMS-phase1.bin +GEO_phase2_BIN = ../CMS-phase2.bin vpath %.cc ${SACMS}/Geoms -all: ${TGTS} ${GEO_2017_BIN} +all: ${TGTS} ${GEO_phase1_BIN} ${GEO_phase2_BIN} %.o: %.cc %.d ${CXX} ${CPPFLAGS} ${CXXFLAGS} ${VEC_HOST} -c -o $@ $< @@ -27,8 +28,11 @@ all: ${TGTS} ${GEO_2017_BIN} ../%.so: %.o ${CXX} -shared -L.. -lMicCore -o $@ $< -${GEO_2017_BIN}: - curl http://xrd-cache-1.t2.ucsd.edu/matevz/PKF/CMS-2017.bin -o $@ +${GEO_phase1_BIN}: + curl http://xrd-cache-1.t2.ucsd.edu/matevz/PKF/CMS-phase1.bin -o $@ + +${GEO_phase2_BIN}: + curl http://xrd-cache-1.t2.ucsd.edu/matevz/PKF/CMS-phase2.bin -o $@ ifeq ($(filter clean distclean, ${MAKECMDGOALS}),) include ${DEPS} diff --git a/RecoTracker/MkFitCMS/standalone/Makefile b/RecoTracker/MkFitCMS/standalone/Makefile index 01b0945bdf82e..401e6118afd19 100644 --- a/RecoTracker/MkFitCMS/standalone/Makefile +++ b/RecoTracker/MkFitCMS/standalone/Makefile @@ -7,10 +7,11 @@ CMS_DIR := ${SRCDIR}/RecoTracker/MkFitCMS LIB_CMS := ../libMicCMS.so MAIN := ../mkFit WRMEMF := ../writeMemoryFile +DICTPCM := ../DictsDict_rdict.pcm TGTS := ${LIB_CMS} ${MAIN} ifdef WITH_ROOT -TGTS += ${WRMEMF} +TGTS += ${WRMEMF} ${DICTPCM} endif .PHONY: all clean distclean @@ -20,6 +21,9 @@ all: ${TGTS} SRCS := $(wildcard ${CMS_DIR}/src/*.cc) $(wildcard ${SACMS}/*.cc) ifdef WITH_ROOT SRCS += ${SACMS}/tkNtuple/WriteMemoryFile.cc +DictsDict.cc ${DICTPCM}: ${SACMS}/tkNtuple/DictsLinkDef.h + rootcint -v3 -f $@ $< + mv DictsDict_rdict.pcm ${DICTPCM} endif SRCB := $(notdir ${SRCS}) DEPS := $(SRCB:.cc=.d) @@ -34,7 +38,7 @@ include ${DEPS} endif clean-local: - -rm -f ${TGTS} *.d *.o *.om *.so + -rm -f ${TGTS} *.d *.o *.om *.so *.pcm -rm -rf main.dSYM -rm -rf plotting/*.so plotting/*.d plotting/*.pcm @@ -51,7 +55,7 @@ ${LIB_CMS}: ${CMS_OBJS} ${MAIN}: ${LIB_CMS} mkFit.o ${CXX} ${CXXFLAGS} ${VEC_HOST} ${LDFLAGS} mkFit.o -o $@ ${LDFLAGS_HOST} -ltbb -L.. -lMicCore -lMicCMS -Wl,-rpath=. -${WRMEMF}: WriteMemoryFile.o +${WRMEMF}: WriteMemoryFile.o DictsDict.o ${CXX} ${CXXFLAGS} ${LDFLAGS} $^ -o $@ ${LDFLAGS_HOST} -ltbb -L.. -lMicCore -Wl,-rpath=. ${OBJS}: %.o: %.cc %.d diff --git a/RecoTracker/MkFitCMS/standalone/MkStandaloneSeqs.cc b/RecoTracker/MkFitCMS/standalone/MkStandaloneSeqs.cc index cee972d963b9b..65a00def35bd6 100644 --- a/RecoTracker/MkFitCMS/standalone/MkStandaloneSeqs.cc +++ b/RecoTracker/MkFitCMS/standalone/MkStandaloneSeqs.cc @@ -29,20 +29,22 @@ namespace mkfit { eoh.setBeamSpot(ev.beamSpot_); } - void handle_duplicates(Event *event) { + void handle_duplicates(Event *) { + /* // Mark tracks as duplicates; if within CMSSW, remove duplicate tracks from fit or candidate track collection if (Config::removeDuplicates) { if (Config::quality_val || Config::sim_val || Config::cmssw_val) { - find_duplicates(event->candidateTracks_); + clean_duplicates(event->candidateTracks_); if (Config::backwardFit) - find_duplicates(event->fitTracks_); + clean_duplicates(event->fitTracks_); } // For the MEIF benchmarks and the stress tests, no validation flags are set so we will enter this block else { // Only care about the candidate tracks here; no need to run the duplicate removal on both candidate and fit tracks - find_duplicates(event->candidateTracks_); + clean_duplicates(event->candidateTracks_); } } + */ } //========================================================================= diff --git a/RecoTracker/MkFitCMS/standalone/buildtestMPlex.cc b/RecoTracker/MkFitCMS/standalone/buildtestMPlex.cc index 3bbd789a140de..6ce5433eb858d 100644 --- a/RecoTracker/MkFitCMS/standalone/buildtestMPlex.cc +++ b/RecoTracker/MkFitCMS/standalone/buildtestMPlex.cc @@ -90,7 +90,7 @@ namespace mkfit { void runBuildingTestPlexDumbCMSSW(Event &ev, const EventOfHits &eoh, MkBuilder &builder) { const IterationConfig &itconf = Config::ItrInfo[0]; - MkJob job({Config::TrkInfo, itconf, eoh}); + MkJob job({Config::TrkInfo, itconf, eoh, eoh.refBeamSpot()}); builder.begin_event(&job, &ev, __func__); @@ -131,7 +131,7 @@ namespace mkfit { ev.fill_hitmask_bool_vectors(itconf.m_track_algorithm, mask_ifc.m_mask_vector); - MkJob job({Config::TrkInfo, itconf, eoh, &mask_ifc}); + MkJob job({Config::TrkInfo, itconf, eoh, eoh.refBeamSpot(), &mask_ifc}); builder.begin_event(&job, &ev, __func__); @@ -160,14 +160,16 @@ namespace mkfit { // Hack, get the tracks out. ev.candidateTracks_ = builder.ref_tracks(); - // For best hit, the candidateTracks_ vector is the direct input to the backward fit so only need to do find_duplicates once + // For best hit, the candidateTracks_ vector is the direct input to the backward fit so only need to do clean_duplicates once if (Config::quality_val || Config::sim_val || Config::cmssw_val) { //Mark tracks as duplicates; if within CMSSW, remove duplicate tracks before backward fit - if (Config::removeDuplicates) { - StdSeq::find_duplicates(ev.candidateTracks_); - } + // CCCC if (Config::removeDuplicates) { + // CCCC StdSeq::clean_duplicates(ev.candidateTracks_); + // CCCC } } + job.switch_to_backward(); + // now do backwards fit... do we want to time this section? if (Config::backwardFit) { builder.backwardFitBH(); @@ -218,7 +220,7 @@ namespace mkfit { ev.fill_hitmask_bool_vectors(itconf.m_track_algorithm, mask_ifc.m_mask_vector); - MkJob job({Config::TrkInfo, itconf, eoh, &mask_ifc}); + MkJob job({Config::TrkInfo, itconf, eoh, eoh.refBeamSpot(), &mask_ifc}); builder.begin_event(&job, &ev, __func__); @@ -248,6 +250,8 @@ namespace mkfit { // first store candidate tracks builder.export_best_comb_cands(ev.candidateTracks_); + job.switch_to_backward(); + // now do backwards fit... do we want to time this section? if (Config::backwardFit) { // Using the TrackVec version until we home in on THE backward fit etc. @@ -259,7 +263,7 @@ namespace mkfit { check_nan_n_silly_bkfit(ev); } - StdSeq::handle_duplicates(&ev); + // CCCC StdSeq::handle_duplicates(&ev); if (Config::quality_val) { StdSeq::Quality qval; @@ -305,7 +309,7 @@ namespace mkfit { ev.fill_hitmask_bool_vectors(itconf.m_track_algorithm, mask_ifc.m_mask_vector); - MkJob job({Config::TrkInfo, itconf, eoh, &mask_ifc}); + MkJob job({Config::TrkInfo, itconf, eoh, eoh.refBeamSpot(), &mask_ifc}); builder.begin_event(&job, &ev, __func__); @@ -335,6 +339,8 @@ namespace mkfit { // first store candidate tracks - needed for BH backward fit and root_validation builder.export_best_comb_cands(ev.candidateTracks_); + job.switch_to_backward(); + // now do backwards fit... do we want to time this section? if (Config::backwardFit) { // a) TrackVec version: @@ -349,7 +355,7 @@ namespace mkfit { check_nan_n_silly_bkfit(ev); } - StdSeq::handle_duplicates(&ev); + // CCCC StdSeq::handle_duplicates(&ev); // validation section if (Config::quality_val) { @@ -421,7 +427,7 @@ namespace mkfit { ev.fill_hitmask_bool_vectors(itconf.m_track_algorithm, mask_ifc.m_mask_vector); - MkJob job({Config::TrkInfo, itconf, eoh, &mask_ifc}); + MkJob job({Config::TrkInfo, itconf, eoh, eoh.refBeamSpot(), &mask_ifc}); builder.begin_event(&job, &ev, __func__); @@ -440,10 +446,10 @@ namespace mkfit { } } - bool do_seed_clean = itconf.m_requires_dupclean_tight; + bool do_seed_clean = bool(itconf.m_seed_cleaner); if (do_seed_clean) - StdSeq::clean_cms_seedtracks_iter(&seeds, itconf, eoh.refBeamSpot()); + itconf.m_seed_cleaner(seeds, itconf, eoh.refBeamSpot()); builder.seed_post_cleaning(seeds); @@ -466,30 +472,24 @@ namespace mkfit { if (validation_on) seeds_used.insert(seeds_used.end(), seeds.begin(), seeds.end()); //cleaned seeds need to be stored somehow - using Algo = TrackBase::TrackAlgorithm; - if (itconf.m_requires_quality_filter && Algo(itconf.m_track_algorithm) != Algo::detachedTripletStep) { - if (Algo(itconf.m_track_algorithm) == Algo::pixelPairStep) { - builder.filter_comb_cands([&](const TrackCand &t) { return StdSeq::qfilter_n_hits_pixseed(t, 3); }); - } else if (Algo(itconf.m_track_algorithm) == Algo::pixelLessStep) { - builder.filter_comb_cands( - [&](const TrackCand &t) { return StdSeq::qfilter_pixelLessFwd(t, eoh.refBeamSpot(), Config::TrkInfo); }); - } else { - builder.filter_comb_cands( - [&](const TrackCand &t) { return StdSeq::qfilter_n_hits(t, itconf.m_params.minHitsQF); }); - } + if (itconf.m_pre_bkfit_filter) { + builder.filter_comb_cands(itconf.m_pre_bkfit_filter); } builder.select_best_comb_cands(); { builder.export_tracks(tmp_tvec); - StdSeq::find_and_remove_duplicates(tmp_tvec, itconf); + if (itconf.m_duplicate_cleaner) + itconf.m_duplicate_cleaner(builder.ref_tracks_nc(), itconf); ev.candidateTracks_.reserve(ev.candidateTracks_.size() + tmp_tvec.size()); for (auto &&t : tmp_tvec) ev.candidateTracks_.emplace_back(std::move(t)); tmp_tvec.clear(); } + job.switch_to_backward(); + // now do backwards fit... do we want to time this section? if (Config::backwardFit) { // a) TrackVec version: @@ -514,23 +514,17 @@ namespace mkfit { builder.endBkwSearch(); } - if (itconf.m_requires_quality_filter && (Algo(itconf.m_track_algorithm) == Algo::detachedTripletStep || - Algo(itconf.m_track_algorithm) == Algo::pixelLessStep)) { - if (Algo(itconf.m_track_algorithm) == Algo::detachedTripletStep) { - builder.filter_comb_cands( - [&](const TrackCand &t) { return StdSeq::qfilter_n_layers(t, eoh.refBeamSpot(), Config::TrkInfo); }); - } else if (Algo(itconf.m_track_algorithm) == Algo::pixelLessStep) { - builder.filter_comb_cands([&](const TrackCand &t) { - return StdSeq::qfilter_pixelLessBkwd(t, eoh.refBeamSpot(), Config::TrkInfo); - }); - } + if (itconf.m_post_bkfit_filter) { + builder.filter_comb_cands(itconf.m_post_bkfit_filter); } - builder.filter_comb_cands([&](const TrackCand &t) { return StdSeq::qfilter_nan_n_silly(t); }); + builder.filter_comb_cands(StdSeq::qfilter_nan_n_silly); builder.select_best_comb_cands(true); // true -> clear m_tracks as they were already filled once above - StdSeq::find_and_remove_duplicates(builder.ref_tracks_nc(), itconf); + if (itconf.m_duplicate_cleaner) + itconf.m_duplicate_cleaner(builder.ref_tracks_nc(), itconf); + builder.export_tracks(ev.fitTracks_); } @@ -538,7 +532,7 @@ namespace mkfit { } // MIMI - Fake back event pointer for final processing (that should be done elsewhere) - MkJob job({Config::TrkInfo, Config::ItrInfo[0], eoh}); + MkJob job({Config::TrkInfo, Config::ItrInfo[0], eoh, eoh.refBeamSpot()}); builder.begin_event(&job, &ev, __func__); if (validation_on) { diff --git a/RecoTracker/MkFitCMS/standalone/mkFit.cc b/RecoTracker/MkFitCMS/standalone/mkFit.cc index 442db10b8aab8..fc1d4cc206209 100644 --- a/RecoTracker/MkFitCMS/standalone/mkFit.cc +++ b/RecoTracker/MkFitCMS/standalone/mkFit.cc @@ -104,6 +104,8 @@ void initGeom() { Config::ItrInfo, Config::json_save_iters_fname_fmt, Config::json_save_iters_include_iter_info_preamble); } + Config::ItrInfo.setupStandardFunctionsFromNames(); + // Test functions for ConfigJsonPatcher // cj.test_Direct (Config::ItrInfo[0]); // cj.test_Patcher(Config::ItrInfo[0]); @@ -152,6 +154,14 @@ namespace { } const char* b2a(bool b) { return b ? "true" : "false"; } + + std::vector deadvectors; + + void init_deadvectors() { + deadvectors.resize(Config::TrkInfo.n_layers()); +#include "RecoTracker/MkFitCMS/standalone/deadmodules.h" + } + } // namespace //============================================================================== @@ -207,6 +217,13 @@ void test_standard() { printf("- reading seeds from file\n"); initGeom(); + if (Config::useDeadModules) { + init_deadvectors(); + } + + if (Config::nEvents <= 0) { + return; + } DataFile data_file; if (g_operation == "read") { @@ -310,8 +327,6 @@ void test_standard() { StdSeq::loadHitsAndBeamSpot(ev, eoh); - std::vector deadvectors(ev.layerHits_.size()); -#include "RecoTracker/MkFitCMS/standalone/deadmodules.h" if (Config::useDeadModules) { StdSeq::loadDeads(eoh, deadvectors); } @@ -584,13 +599,13 @@ int main(int argc, const char* argv[]) { "as reference [eff, FR, DR] (def: %s)\n" " --cmssw-val enable ROOT based validation for building and fitting with CMSSW tracks as " "reference [eff, FR, DR] (def: %s)\n" - " must enable: --geom CMS-2017 --read-cmssw-tracks\n" + " must enable: --geom CMS-phase1 --read-cmssw-tracks\n" " --cmssw-match-fw which cmssw track matching routine to use if validating against CMSSW tracks, " "forward built tracks only (def: %s)\n" - " must enable: --geom CMS-2017 --cmssw-val --read-cmssw-tracks\n" + " must enable: --geom CMS-phase1 --cmssw-val --read-cmssw-tracks\n" " --cmssw-match-bk which cmssw track matching routine to use if validating against CMSSW tracks, " "backward fit tracks only (def: %s)\n" - " must enable: --geom CMS-2017 --cmssw-val --read-cmssw-tracks --backward-fit " + " must enable: --geom CMS-phase1 --cmssw-val --read-cmssw-tracks --backward-fit " "--backward-fit-pca\n" " --inc-shorts include short reco tracks into FR (def: %s)\n" " --keep-hit-info keep vectors of hit idxs and branches in trees (def: %s)\n" @@ -611,16 +626,16 @@ int main(int argc, const char* argv[]) { " == --backward-fit --include-pca\n" " **Seed combo options\n" " --cmssw-simseeds use CMS geom with simtracks for seeds\n" - " == --geom CMS-2017 --seed-input %s\n" + " == --geom CMS-phase1 --seed-input %s\n" " --cmssw-stdseeds use CMS geom with CMSSW seeds uncleaned\n" - " == --geom CMS-2017 --seed-input %s --seed-cleaning %s\n" + " == --geom CMS-phase1 --seed-input %s --seed-cleaning %s\n" " --cmssw-n2seeds use CMS geom with CMSSW seeds cleaned with N^2 routine\n" - " == --geom CMS-2017 --seed-input %s --seed-cleaning %s\n" + " == --geom CMS-phase1 --seed-input %s --seed-cleaning %s\n" " --cmssw-pureseeds use CMS geom with pure CMSSW seeds (seeds which produced CMSSW reco tracks), " "enable read of CMSSW tracks\n" - " == --geom CMS-2017 --seed-input %s --seed-cleaning %s --read-cmssw-tracks\n" + " == --geom CMS-phase1 --seed-input %s --seed-cleaning %s --read-cmssw-tracks\n" " --cmssw-goodlabelseeds use CMS geom with CMSSW seeds with label() >= 0\n" - " == --geom CMS-2017 --seed-input %s --seed-cleaning %s\n" + " == --geom CMS-phase1 --seed-input %s --seed-cleaning %s\n" "\n" " **CMSSW validation combo options\n" " --cmssw-val-fhit-bhit use CMSSW validation with hit based matching (50 percent after seed) for forward " @@ -756,7 +771,8 @@ int main(int argc, const char* argv[]) { printf("List of options for string based inputs \n"); printf( "--geom \n" - " CMS-2017 \n" + " CMS-phase1 \n" + " CMS-phase2 \n" " CylCowWLids \n" "\n"); @@ -935,23 +951,23 @@ int main(int argc, const char* argv[]) { Config::backwardFit = true; Config::includePCA = true; } else if (*i == "--cmssw-simseeds") { - Config::geomPlugin = "CMS-2017"; + Config::geomPlugin = "CMS-phase1"; Config::seedInput = simSeeds; } else if (*i == "--cmssw-stdseeds") { - Config::geomPlugin = "CMS-2017"; + Config::geomPlugin = "CMS-phase1"; Config::seedInput = cmsswSeeds; Config::seedCleaning = noCleaning; } else if (*i == "--cmssw-n2seeds") { - Config::geomPlugin = "CMS-2017"; + Config::geomPlugin = "CMS-phase1"; Config::seedInput = cmsswSeeds; Config::seedCleaning = cleanSeedsN2; } else if (*i == "--cmssw-pureseeds") { - Config::geomPlugin = "CMS-2017"; + Config::geomPlugin = "CMS-phase1"; Config::seedInput = cmsswSeeds; Config::seedCleaning = cleanSeedsPure; Config::readCmsswTracks = true; } else if (*i == "--cmssw-goodlabelseeds") { - Config::geomPlugin = "CMS-2017"; + Config::geomPlugin = "CMS-phase1"; Config::seedInput = cmsswSeeds; Config::seedCleaning = cleanSeedsBadLabel; } else if (*i == "--cmssw-val-fhit-bhit") { diff --git a/RecoTracker/MkFitCMS/standalone/tkNtuple/Makefile b/RecoTracker/MkFitCMS/standalone/tkNtuple/Makefile deleted file mode 100644 index c391026a73c98..0000000000000 --- a/RecoTracker/MkFitCMS/standalone/tkNtuple/Makefile +++ /dev/null @@ -1,57 +0,0 @@ -ifndef ROOTSYS -$(error ROOTSYS is not set. Please set ROOT environment properly) -endif - -WITH_ROOT:=1 - -include ../Makefile.config - -.PHONY: all clean distclean echo - -all: default - -CPPEXTRA := -I.. ${USER_CPPFLAGS} ${DEFS} -LDEXTRA := ${USER_LDFLAGS} - -CPPFLAGS := -I${ROOTSYS}/include ${CPPEXTRA} ${CPPFLAGS} -CXXFLAGS += ${USER_CXXFLAGS} -LDFLAGS += -L${ROOTSYS}/lib -L../lib -lCore -lRIO -lTree -lMathCore ${LDEXTRA} - - -EXES := writeMemoryFile - -default: ${EXES} - -clean: - rm -f ${EXES} *.d *.o *.om *.so *Dict.* *.pcm - rm -rf writeMemoryFile.dSYM - -distclean: clean - rm -f *.optrpt - rm -f ${EXES} *.ah - -echo: - @echo "CXX = ${CXX}" - @echo "CPPFLAGS = ${CPPFLAGS}" - @echo "CXXFLAGS = ${CXXFLAGS}" - @echo "LDFLAGS = ${LDFLAGS}" - @echo "EXES = ${EXES}" - - -################################################################ - -SRCS := WriteMemoryFile.cc - -OBJS := $(SRCS:.cc=.o) - -#libDicts.so: DictsDict.o -# ${CXX} ${CXXFLAGS} ${LDFLAGS} $< -shared -o $@ ${LDFLAGS_HOST} - -DictsDict.cc: DictsLinkDef.h - rootcint -v3 -f $@ -c -p $< - -writeMemoryFile: ${OBJS} DictsDict.o - ${CXX} ${CXXFLAGS} ${LDFLAGS} $^ -L.. -lMicCore -Wl,-rpath=..,-rpath=.,-rpath=../lib -o $@ - -%.o: %.cc - ${CXX} ${CPPFLAGS} ${CXXFLAGS} -c -o $@ $< diff --git a/RecoTracker/MkFitCMS/standalone/tkNtuple/WriteMemoryFile.cc b/RecoTracker/MkFitCMS/standalone/tkNtuple/WriteMemoryFile.cc index 7242066a9e3a3..4c31f9446df1f 100644 --- a/RecoTracker/MkFitCMS/standalone/tkNtuple/WriteMemoryFile.cc +++ b/RecoTracker/MkFitCMS/standalone/tkNtuple/WriteMemoryFile.cc @@ -47,7 +47,7 @@ void printHelp(const char* av0) { "Options:\n" " --input input file\n" " --output output file\n" - " --geo binary TrackerInfo geometry (def: CMS-2017.bin)\n" + " --geo binary TrackerInfo geometry (def: CMS-phase1.bin)\n" " --verbosity print details (0 quiet, 1 print counts, 2 print all; def: 0)\n" " --maxevt maxevt events to write (-1 for everything in the file def: -1)\n" " --clean-sim-tracks apply sim track cleaning (def: no cleaning)\n" @@ -63,7 +63,7 @@ int main(int argc, char* argv[]) { std::string inputFileName; bool haveOutput = false; std::string outputFileName; - std::string geoFileName("CMS-2017.bin"); + std::string geoFileName("CMS-phase1.bin"); bool cleanSimTracks = false; bool writeAllEvents = false; @@ -136,15 +136,6 @@ int main(int argc, char* argv[]) { using namespace std; - TrackerInfo tkinfo; - tkinfo.read_bin_file(geoFileName); - LayerNumberConverter lnc(TkLayout::phase1); - const unsigned int nTotalLayers = lnc.nLayers(); - assert(nTotalLayers == (unsigned int)tkinfo.n_layers()); - - int nstot = 0; - std::vector nhitstot(nTotalLayers, 0); - TFile* f = TFile::Open(inputFileName.c_str()); if (f == 0) { fprintf(stderr, "Failed opening input root file '%s'\n", inputFileName.c_str()); @@ -153,6 +144,16 @@ int main(int argc, char* argv[]) { TTree* t = (TTree*)f->Get("trackingNtuple/tree"); + bool hasPh2hits = t->GetBranch("ph2_isLower") != nullptr; + TrackerInfo tkinfo; + tkinfo.read_bin_file(geoFileName); + LayerNumberConverter lnc(hasPh2hits ? TkLayout::phase2 : TkLayout::phase1); + const unsigned int nTotalLayers = lnc.nLayers(); + assert(nTotalLayers == (unsigned int)tkinfo.n_layers()); + + int nstot = 0; + std::vector nhitstot(nTotalLayers, 0); + unsigned long long event; t->SetBranchAddress("event", &event); @@ -399,26 +400,28 @@ int main(int argc, char* argv[]) { vector* glu_yz = 0; vector* glu_zz = 0; vector* glu_zx = 0; - t->SetBranchAddress("glu_isBarrel", &glu_isBarrel); - if (has910_det_lay) { - t->SetBranchAddress("glu_subdet", &glu_det); - t->SetBranchAddress("glu_layer", &glu_lay); - } else { - t->SetBranchAddress("glu_det", &glu_det); - t->SetBranchAddress("glu_lay", &glu_lay); + if (!hasPh2hits) { + t->SetBranchAddress("glu_isBarrel", &glu_isBarrel); + if (has910_det_lay) { + t->SetBranchAddress("glu_subdet", &glu_det); + t->SetBranchAddress("glu_layer", &glu_lay); + } else { + t->SetBranchAddress("glu_det", &glu_det); + t->SetBranchAddress("glu_lay", &glu_lay); + } + t->SetBranchAddress("glu_detId", &glu_detId); + t->SetBranchAddress("glu_monoIdx", &glu_monoIdx); + t->SetBranchAddress("glu_stereoIdx", &glu_stereoIdx); + t->SetBranchAddress("glu_x", &glu_x); + t->SetBranchAddress("glu_y", &glu_y); + t->SetBranchAddress("glu_z", &glu_z); + t->SetBranchAddress("glu_xx", &glu_xx); + t->SetBranchAddress("glu_xy", &glu_xy); + t->SetBranchAddress("glu_yy", &glu_yy); + t->SetBranchAddress("glu_yz", &glu_yz); + t->SetBranchAddress("glu_zz", &glu_zz); + t->SetBranchAddress("glu_zx", &glu_zx); } - t->SetBranchAddress("glu_detId", &glu_detId); - t->SetBranchAddress("glu_monoIdx", &glu_monoIdx); - t->SetBranchAddress("glu_stereoIdx", &glu_stereoIdx); - t->SetBranchAddress("glu_x", &glu_x); - t->SetBranchAddress("glu_y", &glu_y); - t->SetBranchAddress("glu_z", &glu_z); - t->SetBranchAddress("glu_xx", &glu_xx); - t->SetBranchAddress("glu_xy", &glu_xy); - t->SetBranchAddress("glu_yy", &glu_yy); - t->SetBranchAddress("glu_yz", &glu_yz); - t->SetBranchAddress("glu_zz", &glu_zz); - t->SetBranchAddress("glu_zx", &glu_zx); vector* str_isBarrel = 0; vector* str_isStereo = 0; @@ -438,36 +441,78 @@ int main(int argc, char* argv[]) { vector* str_chargePerCM = 0; vector* str_csize = 0; vector* str_usedMask = 0; - t->SetBranchAddress("str_isBarrel", &str_isBarrel); - t->SetBranchAddress("str_isStereo", &str_isStereo); - if (has910_det_lay) { - t->SetBranchAddress("str_subdet", &str_det); - t->SetBranchAddress("str_layer", &str_lay); - } else { - t->SetBranchAddress("str_det", &str_det); - t->SetBranchAddress("str_lay", &str_lay); - } - t->SetBranchAddress("str_detId", &str_detId); - t->SetBranchAddress("str_simType", &str_simType); - t->SetBranchAddress("str_x", &str_x); - t->SetBranchAddress("str_y", &str_y); - t->SetBranchAddress("str_z", &str_z); - t->SetBranchAddress("str_xx", &str_xx); - t->SetBranchAddress("str_xy", &str_xy); - t->SetBranchAddress("str_yy", &str_yy); - t->SetBranchAddress("str_yz", &str_yz); - t->SetBranchAddress("str_zz", &str_zz); - t->SetBranchAddress("str_zx", &str_zx); - t->SetBranchAddress("str_chargePerCM", &str_chargePerCM); - t->SetBranchAddress("str_clustSize", &str_csize); - if (writeHitIterMasks) { - t->SetBranchAddress("str_usedMask", &str_usedMask); - } - vector>* str_simHitIdx = 0; - t->SetBranchAddress("str_simHitIdx", &str_simHitIdx); vector>* str_chargeFraction = 0; - t->SetBranchAddress("str_chargeFraction", &str_chargeFraction); + if (!hasPh2hits) { + t->SetBranchAddress("str_isBarrel", &str_isBarrel); + t->SetBranchAddress("str_isStereo", &str_isStereo); + if (has910_det_lay) { + t->SetBranchAddress("str_subdet", &str_det); + t->SetBranchAddress("str_layer", &str_lay); + } else { + t->SetBranchAddress("str_det", &str_det); + t->SetBranchAddress("str_lay", &str_lay); + } + t->SetBranchAddress("str_detId", &str_detId); + t->SetBranchAddress("str_simType", &str_simType); + t->SetBranchAddress("str_x", &str_x); + t->SetBranchAddress("str_y", &str_y); + t->SetBranchAddress("str_z", &str_z); + t->SetBranchAddress("str_xx", &str_xx); + t->SetBranchAddress("str_xy", &str_xy); + t->SetBranchAddress("str_yy", &str_yy); + t->SetBranchAddress("str_yz", &str_yz); + t->SetBranchAddress("str_zz", &str_zz); + t->SetBranchAddress("str_zx", &str_zx); + t->SetBranchAddress("str_chargePerCM", &str_chargePerCM); + t->SetBranchAddress("str_clustSize", &str_csize); + if (writeHitIterMasks) { + t->SetBranchAddress("str_usedMask", &str_usedMask); + } + + t->SetBranchAddress("str_simHitIdx", &str_simHitIdx); + t->SetBranchAddress("str_chargeFraction", &str_chargeFraction); + } + + vector* ph2_isLower = 0; + vector* ph2_subdet = 0; + vector* ph2_layer = 0; + vector* ph2_detId = 0; + vector* ph2_simType = 0; + vector* ph2_x = 0; + vector* ph2_y = 0; + vector* ph2_z = 0; + vector* ph2_xx = 0; + vector* ph2_xy = 0; + vector* ph2_yy = 0; + vector* ph2_yz = 0; + vector* ph2_zz = 0; + vector* ph2_zx = 0; + vector* ph2_usedMask = 0; + vector>* ph2_simHitIdx = 0; + if (hasPh2hits && applyCCC) + std::cout << "WARNING: applyCCC is set for Phase2 inputs: applyCCC will be ignored" << std::endl; + if (hasPh2hits) { + t->SetBranchAddress("ph2_isLower", &ph2_isLower); + t->SetBranchAddress("ph2_subdet", &ph2_subdet); + t->SetBranchAddress("ph2_layer", &ph2_layer); + t->SetBranchAddress("ph2_detId", &ph2_detId); + t->SetBranchAddress("ph2_simType", &ph2_simType); + t->SetBranchAddress("ph2_x", &ph2_x); + t->SetBranchAddress("ph2_y", &ph2_y); + t->SetBranchAddress("ph2_z", &ph2_z); + t->SetBranchAddress("ph2_xx", &ph2_xx); + t->SetBranchAddress("ph2_xy", &ph2_xy); + t->SetBranchAddress("ph2_yy", &ph2_yy); + t->SetBranchAddress("ph2_yz", &ph2_yz); + t->SetBranchAddress("ph2_zz", &ph2_zz); + t->SetBranchAddress("ph2_zx", &ph2_zx); + if (writeHitIterMasks) { + t->SetBranchAddress("ph2_usedMask", &ph2_usedMask); + } + t->SetBranchAddress("ph2_simHitIdx", &ph2_simHitIdx); + } + vector ph2_chargeFraction_dummy(16, 0.f); // beam spot float bsp_x; @@ -523,10 +568,12 @@ int main(int argc, char* argv[]) { bs.beamWidthY = bsp_sigmay; //dxdz and dydz are not in the trackingNtuple at the moment - for (unsigned int istr = 0; istr < str_lay->size(); ++istr) { - if (str_chargePerCM->at(istr) < cutValueCCC) - numFailCCC++; - numTotalStr++; + if (!hasPh2hits) { + for (unsigned int istr = 0; istr < str_lay->size(); ++istr) { + if (str_chargePerCM->at(istr) < cutValueCCC) + numFailCCC++; + numTotalStr++; + } } auto nSims = sim_q->size(); @@ -675,6 +722,16 @@ int main(int argc, char* argv[]) { } break; } + case HitType::Phase2OT: { + int istr = ihIdx; + if (istr < 0) + continue; + int cmsswlay = lnc.convertLayerNumber( + ph2_subdet->at(istr), ph2_layer->at(istr), useMatched, ph2_isLower->at(istr), ph2_z->at(istr) > 0); + if (cmsswlay >= 0 && cmsswlay < static_cast(nTotalLayers)) + hitlay[cmsswlay]++; + break; + } case HitType::Invalid: break; //FIXME. Skip, really? default: @@ -729,28 +786,43 @@ int main(int argc, char* argv[]) { if (simTkIdxNt >= 0) nRecToSimHit++; } - if (useMatched) { - for (unsigned int iglu = 0; iglu < glu_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++iglu) { - if (glu_isBarrel->at(iglu) == 0) + if (hasPh2hits) { + for (unsigned int istr = 0; istr < ph2_layer->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++istr) { + int ilay = -1; + ilay = lnc.convertLayerNumber( + ph2_subdet->at(istr), ph2_layer->at(istr), useMatched, ph2_isLower->at(istr), ph2_z->at(istr) > 0); + if (useMatched && !ph2_isLower->at(istr)) + continue; + if (ilay == -1) continue; - int igluMono = glu_monoIdx->at(iglu); - int simTkIdxNt = - bestTkIdx(str_simHitIdx->at(igluMono), str_chargeFraction->at(igluMono), igluMono, HitType::Strip); + int simTkIdxNt = bestTkIdx(ph2_simHitIdx->at(istr), ph2_chargeFraction_dummy, istr, HitType::Phase2OT); + if (simTkIdxNt >= 0) + nRecToSimHit++; + } + } else { + if (useMatched) { + for (unsigned int iglu = 0; iglu < glu_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++iglu) { + if (glu_isBarrel->at(iglu) == 0) + continue; + int igluMono = glu_monoIdx->at(iglu); + int simTkIdxNt = + bestTkIdx(str_simHitIdx->at(igluMono), str_chargeFraction->at(igluMono), igluMono, HitType::Strip); + if (simTkIdxNt >= 0) + nRecToSimHit++; + } + } + for (unsigned int istr = 0; istr < str_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++istr) { + int ilay = -1; + ilay = lnc.convertLayerNumber( + str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0); + if (useMatched && str_isBarrel->at(istr) == 1 && str_isStereo->at(istr)) + continue; + if (ilay == -1) + continue; + int simTkIdxNt = bestTkIdx(str_simHitIdx->at(istr), str_chargeFraction->at(istr), istr, HitType::Strip); if (simTkIdxNt >= 0) nRecToSimHit++; } - } - for (unsigned int istr = 0; istr < str_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++istr) { - int ilay = -1; - ilay = lnc.convertLayerNumber( - str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0); - if (useMatched && str_isBarrel->at(istr) == 1 && str_isStereo->at(istr)) - continue; - if (ilay == -1) - continue; - int simTkIdxNt = bestTkIdx(str_simHitIdx->at(istr), str_chargeFraction->at(istr), istr, HitType::Strip); - if (simTkIdxNt >= 0) - nRecToSimHit++; } if (nRecToSimHit < cleanSimTrack_minRecHits) continue; @@ -766,8 +838,9 @@ int main(int argc, char* argv[]) { vector& seedTracks_ = EE.seedTracks_; vector> pixHitSeedIdx(pix_lay->size()); - vector> strHitSeedIdx(str_lay->size()); - vector> gluHitSeedIdx(glu_lay->size()); + vector> strHitSeedIdx(hasPh2hits ? 0 : str_lay->size()); + vector> gluHitSeedIdx(hasPh2hits ? 0 : glu_lay->size()); + vector> ph2HitSeedIdx(hasPh2hits ? ph2_layer->size() : 0); for (unsigned int is = 0; is < see_q->size(); ++is) { auto isAlgo = TrackAlgorithm(see_algo->at(is)); if (not allSeeds) @@ -843,6 +916,10 @@ int main(int argc, char* argv[]) { } break; } + case HitType::Phase2OT: { + ph2HitSeedIdx[hidx].push_back(seedTracks_.size()); + break; + } case HitType::Invalid: break; //FIXME. Skip, really? default: @@ -857,8 +934,9 @@ int main(int argc, char* argv[]) { vector& cmsswTracks_ = EE.cmsswTracks_; vector> pixHitRecIdx(pix_lay->size()); - vector> strHitRecIdx(str_lay->size()); - vector> gluHitRecIdx(glu_lay->size()); + vector> strHitRecIdx(hasPh2hits ? 0 : str_lay->size()); + vector> gluHitRecIdx(hasPh2hits ? 0 : glu_lay->size()); + vector> ph2HitRecIdx(hasPh2hits ? ph2_layer->size() : 0); for (unsigned int ir = 0; ir < trk_q->size(); ++ir) { //check the origin; redundant for initialStep ntuples if (not allSeeds) @@ -922,17 +1000,22 @@ int main(int argc, char* argv[]) { break; } case HitType::Strip: { - //cout << "pix=" << hidx << " track=" << cmsswTracks_.size() << endl; + //cout << "str=" << hidx << " track=" << cmsswTracks_.size() << endl; strHitRecIdx[hidx].push_back(cmsswTracks_.size()); break; } case HitType::Glued: { if (not useMatched) throw std::logic_error("Tracks have glued hits, but matchedHit load is not configured"); - //cout << "pix=" << hidx << " track=" << cmsswTracks_.size() << endl; + //cout << "glu=" << hidx << " track=" << cmsswTracks_.size() << endl; gluHitRecIdx[hidx].push_back(cmsswTracks_.size()); break; } + case HitType::Phase2OT: { + //cout << "ph2=" << hidx << " track=" << cmsswTracks_.size() << endl; + ph2HitRecIdx[hidx].push_back(cmsswTracks_.size()); + break; + } case HitType::Invalid: break; //FIXME. Skip, really? default: @@ -958,7 +1041,6 @@ int main(int argc, char* argv[]) { int simTkIdxNt = bestTkIdx(pix_simHitIdx->at(ipix), pix_chargeFraction->at(ipix), ipix, HitType::Pixel); int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1; //switch to index in simTracks_ - //cout << Form("pix lay=%i det=%i x=(%6.3f, %6.3f, %6.3f)",ilay+1,pix_det->at(ipix),pix_x->at(ipix),pix_y->at(ipix),pix_z->at(ipix)) << endl; SVector3 pos(pix_x->at(ipix), pix_y->at(ipix), pix_z->at(ipix)); SMatrixSym33 err; @@ -989,105 +1071,155 @@ int main(int argc, char* argv[]) { totHits++; } - if (useMatched) { - for (unsigned int iglu = 0; iglu < glu_lay->size(); ++iglu) { - if (glu_isBarrel->at(iglu) == 0) + if (hasPh2hits) { + vector ph2Idx(ph2_layer->size()); + for (unsigned int iph2 = 0; iph2 < ph2_layer->size(); ++iph2) { + int ilay = -1; + ilay = lnc.convertLayerNumber( + ph2_subdet->at(iph2), ph2_layer->at(iph2), useMatched, ph2_isLower->at(iph2), ph2_z->at(iph2) > 0); + if (useMatched && !ph2_isLower->at(iph2)) continue; - int igluMono = glu_monoIdx->at(iglu); - int simTkIdxNt = - bestTkIdx(str_simHitIdx->at(igluMono), str_chargeFraction->at(igluMono), igluMono, HitType::Strip); + if (ilay == -1) + continue; + + unsigned int imoduleid = tkinfo[ilay].short_id(ph2_detId->at(iph2)); + + int simTkIdxNt = bestTkIdx(ph2_simHitIdx->at(iph2), ph2_chargeFraction_dummy, iph2, HitType::Phase2OT); int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1; //switch to index in simTracks_ - int ilay = lnc.convertLayerNumber(glu_det->at(iglu), glu_lay->at(iglu), useMatched, -1, glu_z->at(iglu) > 0); - // cout << Form("glu lay=%i det=%i bar=%i x=(%6.3f, %6.3f, %6.3f)",ilay+1,glu_det->at(iglu),glu_isBarrel->at(iglu),glu_x->at(iglu),glu_y->at(iglu),glu_z->at(iglu)) << endl; - SVector3 pos(glu_x->at(iglu), glu_y->at(iglu), glu_z->at(iglu)); + SVector3 pos(ph2_x->at(iph2), ph2_y->at(iph2), ph2_z->at(iph2)); SMatrixSym33 err; - err.At(0, 0) = glu_xx->at(iglu); - err.At(1, 1) = glu_yy->at(iglu); - err.At(2, 2) = glu_zz->at(iglu); - err.At(0, 1) = glu_xy->at(iglu); - err.At(0, 2) = glu_zx->at(iglu); - err.At(1, 2) = glu_yz->at(iglu); + err.At(0, 0) = ph2_xx->at(iph2); + err.At(1, 1) = ph2_yy->at(iph2); + err.At(2, 2) = ph2_zz->at(iph2); + err.At(0, 1) = ph2_xy->at(iph2); + err.At(0, 2) = ph2_zx->at(iph2); + err.At(1, 2) = ph2_yz->at(iph2); if (simTkIdx >= 0) { simTracks_[simTkIdx].addHitIdx(layerHits_[ilay].size(), ilay, 0); } - for (unsigned int ir = 0; ir < gluHitSeedIdx[iglu].size(); ir++) { - //cout << "xxx iglu=" << iglu << " seed=" << gluHitSeedIdx[iglu][ir] << endl; - seedTracks_[gluHitSeedIdx[iglu][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known + for (unsigned int ir = 0; ir < ph2HitSeedIdx[iph2].size(); ir++) { + //cout << "xxx iph2=" << iph2 << " seed=" << ph2HitSeedIdx[iph2][ir] << endl; + seedTracks_[ph2HitSeedIdx[iph2][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known } - for (unsigned int ir = 0; ir < gluHitRecIdx[iglu].size(); ir++) { - //cout << "xxx iglu=" << iglu << " recTrack=" << gluHitRecIdx[iglu][ir] << endl; - cmsswTracks_[gluHitRecIdx[iglu][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known + for (unsigned int ir = 0; ir < ph2HitRecIdx[iph2].size(); ir++) { + //cout << "xxx iph2=" << iph2 << " recTrack=" << ph2HitRecIdx[iph2][ir] << endl; + cmsswTracks_[ph2HitRecIdx[iph2][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known } - - // QQQQ module-id-in-layer, adc and phi/theta spans are not done for matched hits. - // Will we ever use / need this? - assert(false && "Implement module-ids, cluster adc and spans for matched hits!"); - Hit hit(pos, err, totHits); + hit.setupAsStrip(imoduleid, 0, 1); layerHits_[ilay].push_back(hit); + if (writeHitIterMasks) + layerHitMasks_[ilay].push_back(ph2_usedMask->at(iph2)); MCHitInfo hitInfo(simTkIdx, ilay, layerHits_[ilay].size() - 1, totHits); simHitsInfo_.push_back(hitInfo); totHits++; } - } + } else { + if (useMatched) { + for (unsigned int iglu = 0; iglu < glu_lay->size(); ++iglu) { + if (glu_isBarrel->at(iglu) == 0) + continue; + int igluMono = glu_monoIdx->at(iglu); + int simTkIdxNt = + bestTkIdx(str_simHitIdx->at(igluMono), str_chargeFraction->at(igluMono), igluMono, HitType::Strip); + int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1; //switch to index in simTracks_ + + int ilay = lnc.convertLayerNumber(glu_det->at(iglu), glu_lay->at(iglu), useMatched, -1, glu_z->at(iglu) > 0); + // cout << Form("glu lay=%i det=%i bar=%i x=(%6.3f, %6.3f, %6.3f)",ilay+1,glu_det->at(iglu),glu_isBarrel->at(iglu),glu_x->at(iglu),glu_y->at(iglu),glu_z->at(iglu)) << endl; + SVector3 pos(glu_x->at(iglu), glu_y->at(iglu), glu_z->at(iglu)); + SMatrixSym33 err; + err.At(0, 0) = glu_xx->at(iglu); + err.At(1, 1) = glu_yy->at(iglu); + err.At(2, 2) = glu_zz->at(iglu); + err.At(0, 1) = glu_xy->at(iglu); + err.At(0, 2) = glu_zx->at(iglu); + err.At(1, 2) = glu_yz->at(iglu); + if (simTkIdx >= 0) { + simTracks_[simTkIdx].addHitIdx(layerHits_[ilay].size(), ilay, 0); + } + for (unsigned int ir = 0; ir < gluHitSeedIdx[iglu].size(); ir++) { + //cout << "xxx iglu=" << iglu << " seed=" << gluHitSeedIdx[iglu][ir] << endl; + seedTracks_[gluHitSeedIdx[iglu][ir]].addHitIdx( + layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known + } + for (unsigned int ir = 0; ir < gluHitRecIdx[iglu].size(); ir++) { + //cout << "xxx iglu=" << iglu << " recTrack=" << gluHitRecIdx[iglu][ir] << endl; + cmsswTracks_[gluHitRecIdx[iglu][ir]].addHitIdx( + layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known + } - vector strIdx; - strIdx.resize(str_lay->size()); - for (unsigned int istr = 0; istr < str_lay->size(); ++istr) { - int ilay = -1; - ilay = lnc.convertLayerNumber( - str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0); - if (useMatched && str_isBarrel->at(istr) == 1 && str_isStereo->at(istr)) - continue; - if (ilay == -1) - continue; + // QQQQ module-id-in-layer, adc and phi/theta spans are not done for matched hits. + // Will we ever use / need this? + assert(false && "Implement module-ids, cluster adc and spans for matched hits!"); - unsigned int imoduleid = tkinfo[ilay].short_id(str_detId->at(istr)); + Hit hit(pos, err, totHits); + layerHits_[ilay].push_back(hit); + MCHitInfo hitInfo(simTkIdx, ilay, layerHits_[ilay].size() - 1, totHits); + simHitsInfo_.push_back(hitInfo); + totHits++; + } + } - int simTkIdxNt = bestTkIdx(str_simHitIdx->at(istr), str_chargeFraction->at(istr), istr, HitType::Strip); - int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1; //switch to index in simTracks_ + vector strIdx; + strIdx.resize(str_lay->size()); + for (unsigned int istr = 0; istr < str_lay->size(); ++istr) { + int ilay = -1; + ilay = lnc.convertLayerNumber( + str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0); + if (useMatched && str_isBarrel->at(istr) == 1 && str_isStereo->at(istr)) + continue; + if (ilay == -1) + continue; - bool passCCC = applyCCC ? (str_chargePerCM->at(istr) > cutValueCCC) : true; + unsigned int imoduleid = tkinfo[ilay].short_id(str_detId->at(istr)); - //if (str_onTrack->at(istr)==0) continue;//do not consider hits that are not on track! - SVector3 pos(str_x->at(istr), str_y->at(istr), str_z->at(istr)); - SMatrixSym33 err; - err.At(0, 0) = str_xx->at(istr); - err.At(1, 1) = str_yy->at(istr); - err.At(2, 2) = str_zz->at(istr); - err.At(0, 1) = str_xy->at(istr); - err.At(0, 2) = str_zx->at(istr); - err.At(1, 2) = str_yz->at(istr); - if (simTkIdx >= 0) { - if (passCCC) - simTracks_[simTkIdx].addHitIdx(layerHits_[ilay].size(), ilay, 0); - else - simTracks_[simTkIdx].addHitIdx(-9, ilay, 0); - } - for (unsigned int ir = 0; ir < strHitSeedIdx[istr].size(); ir++) { - //cout << "xxx istr=" << istr << " seed=" << strHitSeedIdx[istr][ir] << endl; - if (passCCC) - seedTracks_[strHitSeedIdx[istr][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known - else - seedTracks_[strHitSeedIdx[istr][ir]].addHitIdx(-9, ilay, 0); - } - for (unsigned int ir = 0; ir < strHitRecIdx[istr].size(); ir++) { - //cout << "xxx istr=" << istr << " recTrack=" << strHitRecIdx[istr][ir] << endl; - if (passCCC) - cmsswTracks_[strHitRecIdx[istr][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known - else - cmsswTracks_[strHitRecIdx[istr][ir]].addHitIdx(-9, ilay, 0); - } - if (passCCC) { - Hit hit(pos, err, totHits); - hit.setupAsStrip(imoduleid, str_chargePerCM->at(istr), str_csize->at(istr)); - layerHits_[ilay].push_back(hit); - if (writeHitIterMasks) - layerHitMasks_[ilay].push_back(str_usedMask->at(istr)); - MCHitInfo hitInfo(simTkIdx, ilay, layerHits_[ilay].size() - 1, totHits); - simHitsInfo_.push_back(hitInfo); - totHits++; + int simTkIdxNt = bestTkIdx(str_simHitIdx->at(istr), str_chargeFraction->at(istr), istr, HitType::Strip); + int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1; //switch to index in simTracks_ + + bool passCCC = applyCCC ? (str_chargePerCM->at(istr) > cutValueCCC) : true; + + //if (str_onTrack->at(istr)==0) continue;//do not consider hits that are not on track! + SVector3 pos(str_x->at(istr), str_y->at(istr), str_z->at(istr)); + SMatrixSym33 err; + err.At(0, 0) = str_xx->at(istr); + err.At(1, 1) = str_yy->at(istr); + err.At(2, 2) = str_zz->at(istr); + err.At(0, 1) = str_xy->at(istr); + err.At(0, 2) = str_zx->at(istr); + err.At(1, 2) = str_yz->at(istr); + if (simTkIdx >= 0) { + if (passCCC) + simTracks_[simTkIdx].addHitIdx(layerHits_[ilay].size(), ilay, 0); + else + simTracks_[simTkIdx].addHitIdx(-9, ilay, 0); + } + for (unsigned int ir = 0; ir < strHitSeedIdx[istr].size(); ir++) { + //cout << "xxx istr=" << istr << " seed=" << strHitSeedIdx[istr][ir] << endl; + if (passCCC) + seedTracks_[strHitSeedIdx[istr][ir]].addHitIdx( + layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known + else + seedTracks_[strHitSeedIdx[istr][ir]].addHitIdx(-9, ilay, 0); + } + for (unsigned int ir = 0; ir < strHitRecIdx[istr].size(); ir++) { + //cout << "xxx istr=" << istr << " recTrack=" << strHitRecIdx[istr][ir] << endl; + if (passCCC) + cmsswTracks_[strHitRecIdx[istr][ir]].addHitIdx( + layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known + else + cmsswTracks_[strHitRecIdx[istr][ir]].addHitIdx(-9, ilay, 0); + } + if (passCCC) { + Hit hit(pos, err, totHits); + hit.setupAsStrip(imoduleid, str_chargePerCM->at(istr), str_csize->at(istr)); + layerHits_[ilay].push_back(hit); + if (writeHitIterMasks) + layerHitMasks_[ilay].push_back(str_usedMask->at(istr)); + MCHitInfo hitInfo(simTkIdx, ilay, layerHits_[ilay].size() - 1, totHits); + simHitsInfo_.push_back(hitInfo); + totHits++; + } } } @@ -1241,7 +1373,8 @@ int main(int argc, char* argv[]) { il, float(nhitstot[il]) / float(savedEvents)); //Includes those that failed the cluster charge cut - printf("Out of %i hits, %i failed the cut", numTotalStr, numFailCCC); + if (!hasPh2hits) + printf("Out of %i hits, %i failed the cut", numTotalStr, numFailCCC); //======================================================================== diff --git a/RecoTracker/MkFitCore/interface/Config.h b/RecoTracker/MkFitCore/interface/Config.h index f754055add89b..e885792fe8975 100644 --- a/RecoTracker/MkFitCore/interface/Config.h +++ b/RecoTracker/MkFitCore/interface/Config.h @@ -84,7 +84,7 @@ namespace mkfit { constexpr int Niter = 5; constexpr bool useTrigApprox = true; - // Config for Bfield. Note: for now the same for CMS-2017 and CylCowWLids. + // Config for Bfield. Note: for now the same for CMS-phase1 and CylCowWLids. constexpr float Bfield = 3.8112; constexpr float mag_c1 = 3.8114; constexpr float mag_b0 = -3.94991e-06; diff --git a/RecoTracker/MkFitCore/interface/Hit.h b/RecoTracker/MkFitCore/interface/Hit.h index f0534c5da773b..a19f516348558 100644 --- a/RecoTracker/MkFitCore/interface/Hit.h +++ b/RecoTracker/MkFitCore/interface/Hit.h @@ -199,7 +199,7 @@ namespace mkfit { static constexpr int kMinChargePerCM = 1620; static constexpr int kChargePerCMBits = 8; - static constexpr int kDetIdInLayerBits = 12; + static constexpr int kDetIdInLayerBits = 14; static constexpr int kClusterSizeBits = 5; static constexpr int kMaxClusterSize = (1 << kClusterSizeBits) - 1; diff --git a/RecoTracker/MkFitCore/interface/IterationConfig.h b/RecoTracker/MkFitCore/interface/IterationConfig.h index daac1b00fabe7..d045a1c770d98 100644 --- a/RecoTracker/MkFitCore/interface/IterationConfig.h +++ b/RecoTracker/MkFitCore/interface/IterationConfig.h @@ -9,9 +9,12 @@ namespace mkfit { + struct BeamSpot; class EventOfHits; class TrackerInfo; class Track; + class TrackCand; + class MkJob; typedef std::vector TrackVec; @@ -39,6 +42,7 @@ namespace mkfit { class IterationLayerConfig { public: + int m_layer = -1; // Selection limits. float m_select_min_dphi; float m_select_max_dphi; @@ -95,23 +99,8 @@ namespace mkfit { float chi2CutOverlap = 3.5; float pTCutOverlap = 0.0; - //seed cleaning params - float c_ptthr_hpt = 2.0; - //initial - float c_drmax_bh = 0.010; - float c_dzmax_bh = 0.005; - float c_drmax_eh = 0.020; - float c_dzmax_eh = 0.020; - float c_drmax_bl = 0.010; - float c_dzmax_bl = 0.005; - float c_drmax_el = 0.030; - float c_dzmax_el = 0.030; - + //quality filter params int minHitsQF = 4; - float fracSharedHits = 0.19; - float drth_central = 0.001; - float drth_obarrel = 0.001; - float drth_forward = 0.001; //min pT cut float minPtCut = 0.0; @@ -140,23 +129,49 @@ namespace mkfit { class IterationConfig { public: - using partition_seeds_foo = void(const TrackerInfo &, - const TrackVec &, - const EventOfHits &, - IterationSeedPartition &); + // Called directly. + using clean_seeds_cf = int(TrackVec &, const IterationConfig &, const BeamSpot &); + using clean_seeds_func = std::function; + // Called from MkBuilder::find_tracks_load_seeds(). + using partition_seeds_cf = void(const TrackerInfo &, + const TrackVec &, + const EventOfHits &, + IterationSeedPartition &); + using partition_seeds_func = std::function; + // Passed to MkBuilder::filter_comb_cands(). + using filter_candidates_cf = bool(const TrackCand &, const MkJob &); + using filter_candidates_func = std::function; + // Called directly. + using clean_duplicates_cf = void(TrackVec &, const IterationConfig &); + using clean_duplicates_func = std::function; int m_iteration_index = -1; int m_track_algorithm = -1; bool m_requires_seed_hit_sorting = false; - bool m_requires_quality_filter = false; - bool m_requires_dupclean_tight = false; bool m_backward_search = false; bool m_backward_drop_seed_hits = false; int m_backward_fit_min_hits = -1; // Min number of hits to keep when m_backward_drop_seed_hits is true + // seed cleaning params with good defaults (all configurable) + float sc_ptthr_hpt = 2.0; + float sc_drmax_bh = 0.010; + float sc_dzmax_bh = 0.005; + float sc_drmax_eh = 0.020; + float sc_dzmax_eh = 0.020; + float sc_drmax_bl = 0.010; + float sc_dzmax_bl = 0.005; + float sc_drmax_el = 0.030; + float sc_dzmax_el = 0.030; + + // duplicate cleaning params with good defaults (all configurable) + float dc_fracSharedHits = 0.19; + float dc_drth_central = 0.001; + float dc_drth_obarrel = 0.001; + float dc_drth_forward = 0.001; + // Iteration parameters (could be a ptr) IterationParams m_params; IterationParams m_backward_params; @@ -166,7 +181,17 @@ namespace mkfit { std::vector m_steering_params; std::vector m_layer_configs; - std::function m_partition_seeds; + // Standard functions + clean_seeds_func m_seed_cleaner; + partition_seeds_func m_seed_partitioner; + filter_candidates_func m_pre_bkfit_filter, m_post_bkfit_filter; + clean_duplicates_func m_duplicate_cleaner; + + // Names for Standard functions that get saved to / loaded from JSON. + std::string m_seed_cleaner_name; + std::string m_seed_partitioner_name; + std::string m_pre_bkfit_filter_name, m_post_bkfit_filter_name; + std::string m_duplicate_cleaner_name; //---------------------------------------------------------------------------- @@ -179,19 +204,27 @@ namespace mkfit { bool merge_seed_hits_during_cleaning() const { return m_backward_search && m_backward_drop_seed_hits; } - // -------- Setup function + // -------- Setup functions + + void setupStandardFunctionsFromNames(); void cloneLayerSteerCore(const IterationConfig &o) { // Clone common settings for an iteration. // m_iteration_index, m_track_algorithm, cleaning and bkw-search flags, // and IterationParams are not copied. + // Standard functions are also not copied, only their names so one should + // call setupStandardFunctionsFromNames() later on. m_n_regions = o.m_n_regions; m_region_order = o.m_region_order; m_steering_params = o.m_steering_params; m_layer_configs = o.m_layer_configs; - m_partition_seeds = o.m_partition_seeds; + m_seed_cleaner_name = o.m_seed_cleaner_name; + m_seed_partitioner_name = o.m_seed_partitioner_name; + m_pre_bkfit_filter_name = o.m_pre_bkfit_filter_name; + m_post_bkfit_filter_name = o.m_post_bkfit_filter_name; + m_duplicate_cleaner_name = o.m_duplicate_cleaner_name; } void set_iteration_index_and_track_algorithm(int idx, int trk_alg) { @@ -199,23 +232,11 @@ namespace mkfit { m_track_algorithm = trk_alg; } - void set_qf_flags() { - m_requires_seed_hit_sorting = true; - m_requires_quality_filter = true; - } - - void set_qf_params(int minHits, float sharedFrac) { - m_params.minHitsQF = minHits; - m_params.fracSharedHits = sharedFrac; - } - - void set_dupclean_flag() { m_requires_dupclean_tight = true; } - void set_dupl_params(float sharedFrac, float drthCentral, float drthObarrel, float drthForward) { - m_params.fracSharedHits = sharedFrac; - m_params.drth_central = drthCentral; - m_params.drth_obarrel = drthObarrel; - m_params.drth_forward = drthForward; + dc_fracSharedHits = sharedFrac; + dc_drth_central = drthCentral; + dc_drth_obarrel = drthObarrel; + dc_drth_forward = drthForward; } void set_seed_cleaning_params(float pt_thr, @@ -227,15 +248,15 @@ namespace mkfit { float drmax_eh, float dzmax_el, float drmax_el) { - m_params.c_ptthr_hpt = pt_thr; - m_params.c_drmax_bh = drmax_bh; - m_params.c_dzmax_bh = dzmax_bh; - m_params.c_drmax_eh = drmax_eh; - m_params.c_dzmax_eh = dzmax_eh; - m_params.c_drmax_bl = drmax_bl; - m_params.c_dzmax_bl = dzmax_bl; - m_params.c_drmax_el = drmax_el; - m_params.c_dzmax_el = dzmax_el; + sc_ptthr_hpt = pt_thr; + sc_drmax_bh = drmax_bh; + sc_dzmax_bh = dzmax_bh; + sc_drmax_eh = drmax_eh; + sc_dzmax_eh = dzmax_eh; + sc_drmax_bl = drmax_bl; + sc_dzmax_bl = dzmax_bl; + sc_drmax_el = drmax_el; + sc_dzmax_el = dzmax_el; } void set_num_regions_layers(int nreg, int nlay) { @@ -245,7 +266,20 @@ namespace mkfit { for (int i = 0; i < nreg; ++i) m_steering_params[i].m_region = i; m_layer_configs.resize(nlay); + for (int i = 0; i < nlay; ++i) + m_layer_configs[i].m_layer = i; } + + // Catalog of Standard functions + static void register_seed_cleaner(const std::string &name, clean_seeds_func func); + static void register_seed_partitioner(const std::string &name, partition_seeds_func func); + static void register_candidate_filter(const std::string &name, filter_candidates_func func); + static void register_duplicate_cleaner(const std::string &name, clean_duplicates_func func); + + static clean_seeds_func get_seed_cleaner(const std::string &name); + static partition_seeds_func get_seed_partitioner(const std::string &name); + static filter_candidates_func get_candidate_filter(const std::string &name); + static clean_duplicates_func get_duplicate_cleaner(const std::string &name); }; //============================================================================== @@ -264,11 +298,16 @@ namespace mkfit { IterationConfig &operator[](int i) { return m_iterations[i]; } const IterationConfig &operator[](int i) const { return m_iterations[i]; } + + void setupStandardFunctionsFromNames() { + for (auto &i : m_iterations) + i.setupStandardFunctionsFromNames(); + } }; //============================================================================== - // IterationConfig instances are created in Geoms/CMS-2017.cc, Create_CMS_2017(), + // IterationConfig instances are created in Geoms/CMS-phase1.cc, Create_CMS_phase1(), // filling the IterationsInfo object passed in by reference. //============================================================================== @@ -361,7 +400,7 @@ namespace mkfit { // Load a single iteration from JSON file. // This leaves IterationConfig data-members that are not registered // in JSON schema at their default values. - // The only such member is std::function m_partition_seeds. + // There are several std::function members like this. // Assumes JSON file has been saved WITHOUT iteration-info preamble. // Returns a unique_ptr to the cloned IterationConfig. std::unique_ptr load_File(const std::string &fname); diff --git a/RecoTracker/MkFitCore/interface/MkBuilder.h b/RecoTracker/MkFitCore/interface/MkBuilder.h index b55c15602c40a..caf0d3b8423c2 100644 --- a/RecoTracker/MkFitCore/interface/MkBuilder.h +++ b/RecoTracker/MkFitCore/interface/MkBuilder.h @@ -4,6 +4,7 @@ #include "RecoTracker/MkFitCore/interface/IterationConfig.h" #include "RecoTracker/MkFitCore/interface/HitStructures.h" #include "RecoTracker/MkFitCore/interface/TrackStructures.h" +#include "RecoTracker/MkFitCore/interface/MkJob.h" #include #include @@ -24,31 +25,6 @@ namespace mkfit { // MkJob //============================================================================== - class MkJob { - public: - const TrackerInfo &m_trk_info; - // Config &config; // If we want to get rid of namespace / global config - const IterationConfig &m_iter_config; - const EventOfHits &m_event_of_hits; - - const IterationMaskIfcBase *m_iter_mask_ifc = nullptr; - - int num_regions() const { return m_iter_config.m_n_regions; } - const auto regions_begin() const { return m_iter_config.m_region_order.begin(); } - const auto regions_end() const { return m_iter_config.m_region_order.end(); } - - const auto &steering_params(int i) { return m_iter_config.m_steering_params[i]; } - - const auto ¶ms() const { return m_iter_config.m_params; } - const auto ¶ms_bks() const { return m_iter_config.m_backward_params; } - - int max_max_cands() const { return std::max(params().maxCandsPerSeed, params_bks().maxCandsPerSeed); } - - const std::vector *get_mask_for_layer(int layer) { - return m_iter_mask_ifc ? m_iter_mask_ifc->get_mask_for_layer(layer) : nullptr; - } - }; - //============================================================================== // MkBuilder //============================================================================== @@ -56,7 +32,6 @@ namespace mkfit { class MkBuilder { public: using insert_seed_foo = void(const Track &, int, int); - using filter_track_cand_foo = bool(const TrackCand &); typedef std::vector> CandIdx_t; @@ -80,7 +55,7 @@ namespace mkfit { void import_seeds(const TrackVec &in_seeds, const bool seeds_sorted, std::function insert_seed); // filter for rearranging cands that will / will not do backward search. - int filter_comb_cands(std::function filter); + int filter_comb_cands(IterationConfig::filter_candidates_func filter); void find_min_max_hots_size(); diff --git a/RecoTracker/MkFitCore/interface/MkJob.h b/RecoTracker/MkFitCore/interface/MkJob.h new file mode 100644 index 0000000000000..3f9b0dbd56cc7 --- /dev/null +++ b/RecoTracker/MkFitCore/interface/MkJob.h @@ -0,0 +1,40 @@ +#ifndef RecoTracker_MkFitCore_interface_MkJob_h +#define RecoTracker_MkFitCore_interface_MkJob_h + +#include "RecoTracker/MkFitCore/interface/IterationConfig.h" + +namespace mkfit { + + class MkJob { + public: + const TrackerInfo &m_trk_info; + // Config &config; // If we want to get rid of namespace / global config + const IterationConfig &m_iter_config; + const EventOfHits &m_event_of_hits; + const BeamSpot &m_beam_spot; + + const IterationMaskIfcBase *m_iter_mask_ifc = nullptr; + + bool m_in_fwd = true; + void switch_to_backward() { m_in_fwd = false; } + + int num_regions() const { return m_iter_config.m_n_regions; } + const auto regions_begin() const { return m_iter_config.m_region_order.begin(); } + const auto regions_end() const { return m_iter_config.m_region_order.end(); } + + const auto &steering_params(int i) { return m_iter_config.m_steering_params[i]; } + + const auto ¶ms() const { return m_iter_config.m_params; } + const auto ¶ms_bks() const { return m_iter_config.m_backward_params; } + const auto ¶ms_cur() const { return m_in_fwd ? params() : params_bks(); } + + int max_max_cands() const { return std::max(params().maxCandsPerSeed, params_bks().maxCandsPerSeed); } + + const std::vector *get_mask_for_layer(int layer) { + return m_iter_mask_ifc ? m_iter_mask_ifc->get_mask_for_layer(layer) : nullptr; + } + }; + +} // namespace mkfit + +#endif diff --git a/RecoTracker/MkFitCore/src/IterationConfig.cc b/RecoTracker/MkFitCore/src/IterationConfig.cc index f0ddbd06390ea..5025d99cd3fcf 100644 --- a/RecoTracker/MkFitCore/src/IterationConfig.cc +++ b/RecoTracker/MkFitCore/src/IterationConfig.cc @@ -1,10 +1,15 @@ +#include "RecoTracker/MkFitCore/interface/cms_common_macros.h" #include "RecoTracker/MkFitCore/interface/IterationConfig.h" #include "RecoTracker/MkFitCore/interface/Config.h" #include "RecoTracker/MkFitCore/interface/Track.h" +//#define DEBUG +#include "Debug.h" + #include "nlohmann/json.hpp" #include +#include #include #include #include @@ -41,6 +46,7 @@ namespace mkfit { ) ITCONF_DEFINE_TYPE_NON_INTRUSIVE(mkfit::IterationLayerConfig, + /* int */ m_layer, /* float */ m_select_min_dphi, /* float */ m_select_max_dphi, /* float */ m_select_min_dq, @@ -66,47 +72,147 @@ namespace mkfit { /* float */ chi2Cut_min, /* float */ chi2CutOverlap, /* float */ pTCutOverlap, - /* float */ c_ptthr_hpt, - /* float */ c_drmax_bh, - /* float */ c_dzmax_bh, - /* float */ c_drmax_eh, - /* float */ c_dzmax_eh, - /* float */ c_drmax_bl, - /* float */ c_dzmax_bl, - /* float */ c_drmax_el, - /* float */ c_dzmax_el, /* int */ minHitsQF, - /* float */ fracSharedHits, - /* float */ drth_central, - /* float */ drth_obarrel, - /* float */ drth_forward - - ) - - ITCONF_DEFINE_TYPE_NON_INTRUSIVE( - mkfit::IterationConfig, - /* int */ m_iteration_index, - /* int */ m_track_algorithm, - /* bool */ m_requires_seed_hit_sorting, - /* bool */ m_requires_quality_filter, - /* bool */ m_requires_dupclean_tight, - /* bool */ m_backward_search, - /* bool */ m_backward_drop_seed_hits, - /* int */ m_backward_fit_min_hits, - /* mkfit::IterationParams */ m_params, - /* mkfit::IterationParams */ m_backward_params, - /* int */ m_n_regions, - /* vector */ m_region_order, - /* vector */ m_steering_params, - /* vector */ m_layer_configs - // /* function */ m_partition_seeds - ) + /* float */ minPtCut, + /* unsigned int */ maxClusterSize) + + ITCONF_DEFINE_TYPE_NON_INTRUSIVE(mkfit::IterationConfig, + /* int */ m_iteration_index, + /* int */ m_track_algorithm, + /* std::string */ m_seed_cleaner_name, + /* std::string */ m_seed_partitioner_name, + /* std::string */ m_pre_bkfit_filter_name, + /* std::string */ m_post_bkfit_filter_name, + /* std::string */ m_duplicate_cleaner_name, + /* bool */ m_requires_seed_hit_sorting, + /* bool */ m_backward_search, + /* bool */ m_backward_drop_seed_hits, + /* int */ m_backward_fit_min_hits, + /* float */ sc_ptthr_hpt, + /* float */ sc_drmax_bh, + /* float */ sc_dzmax_bh, + /* float */ sc_drmax_eh, + /* float */ sc_dzmax_eh, + /* float */ sc_drmax_bl, + /* float */ sc_dzmax_bl, + /* float */ sc_drmax_el, + /* float */ sc_dzmax_el, + /* float */ dc_fracSharedHits, + /* float */ dc_drth_central, + /* float */ dc_drth_obarrel, + /* float */ dc_drth_forward, + /* mkfit::IterationParams */ m_params, + /* mkfit::IterationParams */ m_backward_params, + /* int */ m_n_regions, + /* vector */ m_region_order, + /* vector */ m_steering_params, + /* vector */ m_layer_configs) ITCONF_DEFINE_TYPE_NON_INTRUSIVE(mkfit::IterationsInfo, /* vector */ m_iterations) // End AUTO code. + // Begin IterationConfig function catalogs + + namespace { + struct FuncCatalog { + std::map seed_cleaners; + std::map seed_partitioners; + std::map candidate_filters; + std::map duplicate_cleaners; + + std::mutex catalog_mutex; + }; + + FuncCatalog &get_catalog() { + CMS_SA_ALLOW static FuncCatalog func_catalog; + return func_catalog; + } + } // namespace + +#define GET_FC \ + auto &fc = get_catalog(); \ + const std::lock_guard lock(fc.catalog_mutex) + + void IterationConfig::register_seed_cleaner(const std::string &name, clean_seeds_func func) { + GET_FC; + fc.seed_cleaners.insert({name, func}); + } + void IterationConfig::register_seed_partitioner(const std::string &name, partition_seeds_func func) { + GET_FC; + fc.seed_partitioners.insert({name, func}); + } + void IterationConfig::register_candidate_filter(const std::string &name, filter_candidates_func func) { + GET_FC; + fc.candidate_filters.insert({name, func}); + } + void IterationConfig::register_duplicate_cleaner(const std::string &name, clean_duplicates_func func) { + GET_FC; + fc.duplicate_cleaners.insert({name, func}); + } + + namespace { + template + typename T::mapped_type resolve_func_name(const T &cont, const std::string &name, const char *func) { + if (name.empty()) { + return nullptr; + } + auto ii = cont.find(name); + if (ii == cont.end()) { + std::string es(func); + es += " '" + name + "' not found in function registry."; + throw std::runtime_error(es); + } + return ii->second; + } + } // namespace + + IterationConfig::clean_seeds_func IterationConfig::get_seed_cleaner(const std::string &name) { + GET_FC; + return resolve_func_name(fc.seed_cleaners, name, __func__); + } + IterationConfig::partition_seeds_func IterationConfig::get_seed_partitioner(const std::string &name) { + GET_FC; + return resolve_func_name(fc.seed_partitioners, name, __func__); + } + IterationConfig::filter_candidates_func IterationConfig::get_candidate_filter(const std::string &name) { + GET_FC; + return resolve_func_name(fc.candidate_filters, name, __func__); + } + IterationConfig::clean_duplicates_func IterationConfig::get_duplicate_cleaner(const std::string &name) { + GET_FC; + return resolve_func_name(fc.duplicate_cleaners, name, __func__); + } + +#undef GET_FC + + // End IterationConfig function catalogs + + void IterationConfig::setupStandardFunctionsFromNames() { + m_seed_cleaner = get_seed_cleaner(m_seed_cleaner_name); + dprintf(" Set seed_cleaner for '%s' %s\n", m_seed_cleaner_name.c_str(), m_seed_cleaner ? "SET" : "NOT SET"); + + m_seed_partitioner = get_seed_partitioner(m_seed_partitioner_name); + dprintf( + " Set seed_partitioner for '%s' %s\n", m_seed_partitioner_name.c_str(), m_seed_partitioner ? "SET" : "NOT SET"); + + m_pre_bkfit_filter = get_candidate_filter(m_pre_bkfit_filter_name); + dprintf( + " Set pre_bkfit_filter for '%s' %s\n", m_pre_bkfit_filter_name.c_str(), m_pre_bkfit_filter ? "SET" : "NOT SET"); + + m_post_bkfit_filter = get_candidate_filter(m_post_bkfit_filter_name); + dprintf(" Set post_bkfit_filter for '%s' %s\n", + m_post_bkfit_filter_name.c_str(), + m_post_bkfit_filter ? "SET" : "NOT SET"); + + m_duplicate_cleaner = get_duplicate_cleaner(m_duplicate_cleaner_name); + dprintf(" Set duplicate_cleaner for '%s' %s\n", + m_duplicate_cleaner_name.c_str(), + m_duplicate_cleaner ? "SET" : "NOT SET"); + ; + } + // ============================================================================ // ConfigJsonPatcher // ============================================================================ diff --git a/RecoTracker/MkFitCore/src/MkBuilder.cc b/RecoTracker/MkFitCore/src/MkBuilder.cc index 8a0f6af8056f7..de249104b50fa 100644 --- a/RecoTracker/MkFitCore/src/MkBuilder.cc +++ b/RecoTracker/MkFitCore/src/MkBuilder.cc @@ -20,6 +20,7 @@ //#define DEBUG #include "Debug.h" +//#define DEBUG_FINAL_FIT #include "oneapi/tbb/parallel_for.h" #include "oneapi/tbb/parallel_for_each.h" @@ -244,11 +245,11 @@ namespace mkfit { part.m_phi_eta_foo = [&](float phi, float eta) { phi_eta_binnor.register_entry_safe(phi, eta); }; phi_eta_binnor.begin_registration(size); - m_job->m_iter_config.m_partition_seeds(m_job->m_trk_info, in_seeds, m_job->m_event_of_hits, part); + m_job->m_iter_config.m_seed_partitioner(m_job->m_trk_info, in_seeds, m_job->m_event_of_hits, part); phi_eta_binnor.finalize_registration(); ranks.swap(phi_eta_binnor.m_ranks); } else { - m_job->m_iter_config.m_partition_seeds(m_job->m_trk_info, in_seeds, m_job->m_event_of_hits, part); + m_job->m_iter_config.m_seed_partitioner(m_job->m_trk_info, in_seeds, m_job->m_event_of_hits, part); } for (int i = 0; i < size; ++i) { @@ -294,13 +295,13 @@ namespace mkfit { m_seedEtaSeparators[2] - m_seedEtaSeparators[1], m_seedMinLastLayer[2], m_seedMaxLastLayer[2], m_seedEtaSeparators[3] - m_seedEtaSeparators[2], m_seedMinLastLayer[3], m_seedMaxLastLayer[3], m_seedEtaSeparators[4] - m_seedEtaSeparators[3], m_seedMinLastLayer[4], m_seedMaxLastLayer[4]); - dcall(print_seeds(m_event_of_comb_cands)); + // dcall(print_seeds(m_event_of_comb_cands)); // clang-format on } //------------------------------------------------------------------------------ - int MkBuilder::filter_comb_cands(std::function filter) { + int MkBuilder::filter_comb_cands(IterationConfig::filter_candidates_func filter) { EventOfCombCandidates &eoccs = m_event_of_comb_cands; int i = 0, place_pos = 0; @@ -308,7 +309,7 @@ namespace mkfit { std::vector removed_cnts(m_job->num_regions()); while (i < eoccs.size()) { - if (filter(eoccs[i].front())) { + if (filter(eoccs[i].front(), *m_job)) { if (place_pos != i) std::swap(eoccs[place_pos], eoccs[i]); ++place_pos; @@ -1178,14 +1179,14 @@ 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 = PropagationConfig::get_default(); -#ifdef DEBUG +#ifdef DEBUG_FINAL_FIT EventOfCombCandidates &eoccs = m_event_of_comb_cands; #endif for (int icand = start_cand; icand < end_cand; icand += NN) { const int end = std::min(icand + NN, end_cand); -#ifdef DEBUG +#ifdef DEBUG_FINAL_FIT printf("Pre Final fit for %d - %d\n", icand, end); for (int i = icand; i < end; ++i) { const TrackCand &t = eoccs[i][0]; @@ -1268,7 +1269,7 @@ namespace mkfit { // copy out full set of info at last propagated position mkfndr->bkFitOutputTracks(m_tracks, icand, end, prop_config.backward_fit_to_pca); -#ifdef DEBUG +#ifdef DEBUG_FINAL_FIT printf("Post Final fit for %d - %d\n", icand, end); for (int i = icand; i < end; ++i) { const TrackCand &t = eoccs[i][0]; @@ -1321,7 +1322,7 @@ namespace mkfit { for (int icand = start_cand; icand < end_cand; icand += step) { int end = std::min(icand + NN, end_cand); -#ifdef DEBUG +#ifdef DEBUG_FINAL_FIT printf("Pre Final fit for %d - %d\n", icand, end); for (int i = icand; i < end; ++i) { const TrackCand &t = eoccs[i][0]; @@ -1369,7 +1370,7 @@ namespace mkfit { mkfndr->bkFitOutputTracks(eoccs, icand, end, prop_config.backward_fit_to_pca); -#ifdef DEBUG +#ifdef DEBUG_FINAL_FIT printf("Post Final fit for %d - %d\n", icand, end); for (int i = icand; i < end; ++i) { const TrackCand &t = eoccs[i][0]; diff --git a/RecoTracker/MkFitCore/src/radix_sort.cc b/RecoTracker/MkFitCore/src/radix_sort.cc index 5a7f184632d27..132db04da6c4c 100644 --- a/RecoTracker/MkFitCore/src/radix_sort.cc +++ b/RecoTracker/MkFitCore/src/radix_sort.cc @@ -25,8 +25,8 @@ namespace mkfit { template void radix_sort::histo_loop(const std::vector& values, rank_t* histos) { // Create histograms (counters). Counters for all passes are created in one run. - ubyte_t* p = (ubyte_t*)values.data(); - ubyte_t* pe = p + (values.size() * c_NBytes); + const ubyte_t* p = (const ubyte_t*)values.data(); + const ubyte_t* pe = p + (values.size() * c_NBytes); std::array ha; for (rank_t j = 0; j < c_NBytes; ++j) ha[j] = &histos[j << 8]; @@ -51,7 +51,7 @@ namespace mkfit { rank_t* cur_count = &histos[j << 8]; // Get first byte - if that byte's counter equals nb, all values are the same. - ubyte_t first_entry_val = *(((ubyte_t*)values.data()) + j); + const ubyte_t first_entry_val = *(((const ubyte_t*)values.data()) + j); if (cur_count[first_entry_val] != nb) { // Create offsets link[0] = ranks2.data(); diff --git a/RecoTracker/MkFitCore/standalone/Event.cc b/RecoTracker/MkFitCore/standalone/Event.cc index 16722d02134ab..1a5d8259e16bb 100644 --- a/RecoTracker/MkFitCore/standalone/Event.cc +++ b/RecoTracker/MkFitCore/standalone/Event.cc @@ -849,8 +849,8 @@ namespace mkfit { //============================================================================== int DataFile::openRead(const std::string &fname, int expected_n_layers) { - constexpr int min_ver = 4; - constexpr int max_ver = 6; + constexpr int min_ver = 7; + constexpr int max_ver = 7; f_fp = fopen(fname.c_str(), "r"); assert(f_fp != 0 && "Opening of input file failed."); diff --git a/RecoTracker/MkFitCore/standalone/Event.h b/RecoTracker/MkFitCore/standalone/Event.h index aacc45daaa5a4..e5b59c258b59f 100644 --- a/RecoTracker/MkFitCore/standalone/Event.h +++ b/RecoTracker/MkFitCore/standalone/Event.h @@ -72,7 +72,7 @@ namespace mkfit { struct DataFileHeader { int f_magic = 0xBEEF; - int f_format_version = 6; + int f_format_version = 7; //last update with ph2 geom int f_sizeof_track = sizeof(Track); int f_sizeof_hit = sizeof(Hit); int f_sizeof_hot = sizeof(HitOnTrack); diff --git a/RecoTracker/MkFitCore/standalone/Makefile b/RecoTracker/MkFitCore/standalone/Makefile index 0a6b72ba67868..40e2ce6f91083 100644 --- a/RecoTracker/MkFitCore/standalone/Makefile +++ b/RecoTracker/MkFitCore/standalone/Makefile @@ -13,14 +13,13 @@ TGTS := ${LIB_CORE} all: ${TGTS} SRCS := $(wildcard ${CORE_DIR}/src/*.cc) \ - $(wildcard ${CORE_DIR}/src/Ice/*.cc) \ $(wildcard ${CORE_DIR}/src/Matriplex/*.cc) \ $(wildcard ${SADIR}/*.cc) SRCB := $(notdir ${SRCS}) DEPS := $(SRCB:.cc=.d) OBJS := $(SRCB:.cc=.o) -vpath %.cc ${CORE_DIR}/src ${CORE_DIR}/src/Ice ${CORE_DIR}/src/Matriplex ${SADIR} +vpath %.cc ${CORE_DIR}/src ${CORE_DIR}/src/Matriplex ${SADIR} AUTO_TGTS := diff --git a/RecoTracker/MkFitCore/standalone/plotting/Common.hh b/RecoTracker/MkFitCore/standalone/plotting/Common.hh index 445e275b942e2..171b9e86102c6 100644 --- a/RecoTracker/MkFitCore/standalone/plotting/Common.hh +++ b/RecoTracker/MkFitCore/standalone/plotting/Common.hh @@ -12,6 +12,7 @@ #include "TAxis.h" #include "TH1F.h" #include "TF1.h" +#include "TSystem.h" #include #include @@ -257,7 +258,8 @@ namespace { builds.emplace_back(buildsMap["BH"]); builds.emplace_back(buildsMap["CE"]); } else { - builds.emplace_back(buildsMap["STD"]); + if (!gSystem->Getenv("MKFIT_MIMI")) // MIMI does not support STD + builds.emplace_back(buildsMap["STD"]); builds.emplace_back(buildsMap["CE"]); } } else if (SUITE == forConf) { diff --git a/RecoTracker/MkFitCore/standalone/val_scripts/validation-cmssw-benchmarks-multiiter.sh b/RecoTracker/MkFitCore/standalone/val_scripts/validation-cmssw-benchmarks-multiiter.sh index b9eb25f6f8b13..ae9337193624e 100755 --- a/RecoTracker/MkFitCore/standalone/val_scripts/validation-cmssw-benchmarks-multiiter.sh +++ b/RecoTracker/MkFitCore/standalone/val_scripts/validation-cmssw-benchmarks-multiiter.sh @@ -17,6 +17,8 @@ source xeon_scripts/init-env.sh export MIMI="CE mimi" declare -a val_builds=(MIMI) nevents=250 +extraargs="--use-dead-modules" +numiters=10 ## Common file setup case ${inputBin} in @@ -40,6 +42,17 @@ case ${inputBin} in nevents=20000 sample=10mu ;; +"TTbar_phase2") + dir=/home/matevz/mic-dev + subdir= + file=ttbar-p2.bin + nevents=100 + extraargs="--geom CMS-phase2" + numiters=1 + # Pass MIMI flag to ROOT plotting functions -- some will insert STD + # that does not work with MIMI at this point. + export MKFIT_MIMI=1 + ;; *) echo "INPUT BIN IS UNKNOWN" exit 12 @@ -47,16 +60,26 @@ case ${inputBin} in esac ## Common executable setup -maxth=64 -maxvu=16 -maxev=32 +if [[ `lsb_release -si` == "Fedora" ]] +then + maxth=16 # 64 + maxvu=8 # 16 + maxev=16 # 32 + avxmakearg="AVX2:=1" +else + maxth=64 + maxvu=16 + maxev=32 + avxmakearg="AVX_512:=1" +fi + if [[ "${suite}" == "valMT1" ]] then maxth=1 maxev=1 fi seeds="--cmssw-n2seeds" -exe="./mkFit/mkFit --silent ${seeds} --num-thr ${maxth} --num-thr-ev ${maxev} --input-file ${dir}/${subdir}/${file} --num-events ${nevents} --remove-dup --use-dead-modules" +exe="./mkFit --silent ${seeds} --num-thr ${maxth} --num-thr-ev ${maxev} --input-file ${dir}/${subdir}/${file} --num-events ${nevents} --remove-dup ${extraargs}" ## Common output setup tmpdir="tmp" @@ -69,8 +92,8 @@ siminfo="--try-to-save-sim-info" bkfit="--backward-fit" ## validation options: SIMVAL == sim tracks as reference, CMSSWVAL == cmssw tracks as reference -SIMVAL="SIMVAL --sim-val ${siminfo} ${bkfit} ${style} --num-iters-cmssw 10" -SIMVAL_SEED="SIMVALSEED --sim-val ${siminfo} ${bkfit} --mtv-require-seeds --num-iters-cmssw 10" +SIMVAL="SIMVAL --sim-val ${siminfo} ${bkfit} ${style} --num-iters-cmssw ${numiters}" +SIMVAL_SEED="SIMVALSEED --sim-val ${siminfo} ${bkfit} --mtv-require-seeds --num-iters-cmssw ${numiters}" declare -a vals=(SIMVAL SIMVAL_SEED) @@ -98,11 +121,16 @@ SIMPLOTSEED10="SIMVALSEED iter10 0 10 0" SIMPLOT6="SIMVAL iter6 0 6 0" SIMPLOTSEED6="SIMVALSEED iter6 0 6 0" -declare -a plots=(SIMPLOT4 SIMPLOTSEED4 SIMPLOT22 SIMPLOTSEED22 SIMPLOT23 SIMPLOTSEED23 SIMPLOT5 SIMPLOTSEED5 SIMPLOT24 SIMPLOTSEED24 SIMPLOT7 SIMPLOTSEED7 SIMPLOT8 SIMPLOTSEED8 SIMPLOT9 SIMPLOTSEED9 SIMPLOT10 SIMPLOTSEED10 SIMPLOT6 SIMPLOTSEED6) +if [[ "${inputBin}" == "TTbar_phase2" ]] +then + declare -a plots=(SIMPLOT4 SIMPLOTSEED4) +else + declare -a plots=(SIMPLOT4 SIMPLOTSEED4 SIMPLOT22 SIMPLOTSEED22 SIMPLOT23 SIMPLOTSEED23 SIMPLOT5 SIMPLOTSEED5 SIMPLOT24 SIMPLOTSEED24 SIMPLOT7 SIMPLOTSEED7 SIMPLOT8 SIMPLOTSEED8 SIMPLOT9 SIMPLOTSEED9 SIMPLOT10 SIMPLOTSEED10 SIMPLOT6 SIMPLOTSEED6) +fi ## special cmssw dummy build -CMSSW="CMSSW cmssw SIMVAL --sim-val-for-cmssw ${siminfo} --read-cmssw-tracks ${style} --num-iters-cmssw 10" -CMSSW2="CMSSW cmssw SIMVALSEED --sim-val-for-cmssw ${siminfo} --read-cmssw-tracks --mtv-require-seeds --num-iters-cmssw 10" +CMSSW="CMSSW cmssw SIMVAL --sim-val-for-cmssw ${siminfo} --read-cmssw-tracks ${style} --num-iters-cmssw ${numiters}" +CMSSW2="CMSSW cmssw SIMVALSEED --sim-val-for-cmssw ${siminfo} --read-cmssw-tracks --mtv-require-seeds --num-iters-cmssw ${numiters}" ############### ## Functions ## @@ -154,8 +182,8 @@ function plotVal() ######################## ## Compile once -make clean -mVal="-j 32 WITH_ROOT:=1 AVX_512:=1" +make distclean +mVal="-j 32 WITH_ROOT:=1 ${avxmakearg}" make ${mVal} mkdir -p ${tmpdir} @@ -184,12 +212,10 @@ do echo ${!val} | while read -r vN vO done ## clean up -make clean ${mVal} mv tmp/valtree_*.root . rm -rf ${tmpdir} - ## Compute observables and make images for plot in "${plots[@]}" do echo ${!plot} | while read -r pN suff pO iter cancel @@ -223,15 +249,12 @@ do echo ${!plot} | while read -r pN suff pO iter cancel echo "Overlaying histograms for: ${base} ${vN}" if [[ "${suff}" == "all" ]] then - root -b -q -l plotting/makeValidation.C\(\"${base}\",\"_${pN}\",${pO},\"${suite}\"\) + root -b -q -l plotting/makeValidation.C\(\"${base}\",\"_${pN}\",${pO},\"${suite}\"\) else root -b -q -l plotting/makeValidation.C\(\"${base}\",\"_${pN}_${suff}\",${pO},\"${suite}\"\) fi done done -## Final cleanup -make distclean ${mVal} - ## Final message echo "Finished physics validation!" diff --git a/RecoTracker/MkFitCore/standalone/validation-desc.txt b/RecoTracker/MkFitCore/standalone/validation-desc.txt index 3cc06beb75e48..9112e3d33b1ac 100644 --- a/RecoTracker/MkFitCore/standalone/validation-desc.txt +++ b/RecoTracker/MkFitCore/standalone/validation-desc.txt @@ -6,6 +6,26 @@ EDIT HISTORY PREFACE: This file is a compendium on how the validation runs within mkFit, which makes use of the TTreeValidation class and other supporting macros. +2022-07-13 MT notes +=================== + +I have run the MIMI forPR validation in standalone build for the first time +since migration to CMSSW. I have not tried anything else. + +Since standalone build is now done out-of-source, I did the following +in the top-level standalone build directory (where mkFit executable is): + +ln -s ../RecoTracker/MkFitCore/standalone/xeon_scripts . +ln -s ../RecoTracker/MkFitCore/standalone/val_scripts . +ln -s ../RecoTracker/MkFitCore/standalone/plotting . +ln -s ../RecoTracker/MkFitCore/standalone/web . + +rm -rf valtree_*.root log_*.txt +val_scripts/validation-cmssw-benchmarks-multiiter.sh forPR --mtv-like-val TTbar_Phase2 + +web/collectBenchmarks-multi.sh Phase2-0 forPR + + =================== Table of Contents =================== diff --git a/RecoTracker/MkFitCore/standalone/xeon_scripts/init-env.sh b/RecoTracker/MkFitCore/standalone/xeon_scripts/init-env.sh index ce25278a5ec0f..e6795b0983314 100644 --- a/RecoTracker/MkFitCore/standalone/xeon_scripts/init-env.sh +++ b/RecoTracker/MkFitCore/standalone/xeon_scripts/init-env.sh @@ -1,4 +1,9 @@ #!/bin/bash -source /cvmfs/cms.cern.ch/slc7_amd64_gcc10/lcg/root/6.24.07-f52350f4e0b802edeb9a2551a7d00b92/etc/profile.d/init.sh -export TBB_GCC=/cvmfs/cms.cern.ch/slc7_amd64_gcc10/external/tbb/v2021.4.0-d0152ca29055e3a1bbf629673f6e97c4 - +if [[ `lsb_release -si` == "Fedora" ]] +then + source /cvmfs/cms.cern.ch/el8_amd64_gcc10/lcg/root/6.24.07-da610b2b7ed663a0a05d3605f3d83ceb/etc/profile.d/init.sh + export TBB_GCC=/cvmfs/cms.cern.ch/el8_amd64_gcc10/external/tbb/v2021.5.0-e966a5acb1e4d5fd7605074bafbb079c/ +else + source /cvmfs/cms.cern.ch/slc7_amd64_gcc10/lcg/root/6.24.07-f52350f4e0b802edeb9a2551a7d00b92/etc/profile.d/init.sh + export TBB_GCC=/cvmfs/cms.cern.ch/slc7_amd64_gcc10/external/tbb/v2021.4.0-d0152ca29055e3a1bbf629673f6e97c4 +fi