Candidate ranking à la CMSSW

Config.h

@@ -344,6 +344,10 @@ namespace Config
   constexpr float treg_eta[2] = {0.45,1.5};
   constexpr float track_ptlow = 0.9;
 
+  // sorting config (bonus,penalty)
+  constexpr float validHitBonus_ = 2.5;
+  constexpr float missingHitPenalty_ = 20.0;
+
   // Threading
   extern int    numThreadsFinder;
   extern int    numThreadsSimulation;

TTreeValidation.cc

@@ -880,7 +880,8 @@ void TTreeValidation::mapRefTkToRecoTks(const TrackVec& evt_tracks, TrackExtraVe
       {
 	tmpMatches.emplace_back(evt_tracks[label]);
       }
-      std::sort(tmpMatches.begin(), tmpMatches.end(), sortByHitsChi2); // sort the tracks
+      //std::sort(tmpMatches.begin(), tmpMatches.end(), sortByHitsChi2); // sort the tracks
+      std::sort(tmpMatches.begin(), tmpMatches.end(), sortByScoreCand); // sort the tracks
       for (auto itrack = 0; itrack < (int) tmpMatches.size(); itrack++) // loop over sorted tracks, now set the vector of sorted labels match
       {
 	refTkMatches.second[itrack] = tmpMatches[itrack].label();

Track.h

@@ -25,6 +25,17 @@ inline int calculateCharge(const float hit0_x, const float hit0_y,
   return ((hit2_y-hit0_y)*(hit2_x-hit1_x)>(hit2_y-hit1_y)*(hit2_x-hit0_x)?1:-1);
 }
 
+struct IdxChi2List
+{
+public:
+  int   trkIdx; // candidate index
+  int   hitIdx; // hit index
+  int   nhits;  // number of hits (used for sorting)
+  int   nholes;  // number of holes (used for sorting)
+  float pt;   // pt (used for sorting)
+  float chi2;   // total chi2 (used for sorting)
+};
+
 struct ReducedTrack // used for cmssw reco track validation
 {
 public:
@@ -467,6 +478,49 @@ inline bool sortByHitsChi2(const Track & cand1, const Track & cand2)
   return cand1.nFoundHits()>cand2.nFoundHits();
 }
 
+inline bool sortByScoreLoop(const int nfoundhits[2], 
+			const int nmisshits[2], 
+			const float chi2[2],
+			const float pt[2])
+{
+  float score[2] = {0.f,0.f};
+  for(int c=0; c<2; ++c){
+    score[c] = Config::validHitBonus_*nfoundhits[c] - Config::missingHitPenalty_*nmisshits[c] - chi2[c];
+    if(pt[c]<0.9f) score[c] -= 0.5f*Config::validHitBonus_*nfoundhits[c];
+    else if(nfoundhits[c]>8) score[c] += Config::validHitBonus_*nfoundhits[c];
+  }
+  return score[0]>score[1];
+}
+
+inline bool sortByScoreCand(const Track& cand1, const Track& cand2)
+{
+  int nfoundhits[2] = {cand1.nFoundHits(),cand2.nFoundHits()};
+  int nmisshits[2] = {cand1.nTotalHits()-cand1.nFoundHits(),cand2.nTotalHits()-cand2.nFoundHits()};
+  float chi2[2] = {cand1.chi2(),cand2.chi2()};
+  float pt[2] = {cand1.pT(),cand2.pT()};
+  return sortByScoreLoop(nfoundhits,nmisshits,chi2,pt);
+}
+
+inline bool sortByScoreStruct(const IdxChi2List& cand1, const IdxChi2List& cand2)
+{
+
+  int nfoundhits[2] = {cand1.nhits,cand2.nhits};
+  int nmisshits[2] = {cand1.nholes,cand2.nholes};
+  float chi2[2] = {cand1.chi2,cand2.chi2};
+  float pt[2] = {cand1.pt,cand2.pt};
+  return sortByScoreLoop(nfoundhits,nmisshits,chi2,pt);
+}
+
+inline bool sortByScoreCandPair(const std::pair<Track, TrackState>& cand1, const std::pair<Track, TrackState>& cand2)
+{
+  int nfoundhits[2] = {cand1.first.nFoundHits(),cand2.first.nFoundHits()};
+  int nmisshits[2] = {cand1.first.nTotalHits()-cand1.first.nFoundHits(),cand2.first.nTotalHits()-cand2.first.nFoundHits()};
+  float chi2[2] = {cand1.first.chi2(),cand2.first.chi2()};
+  float pt[2] = {cand1.first.pT(),cand2.first.pT()};
+  return sortByScoreLoop(nfoundhits,nmisshits,chi2,pt);
+}
+
+
 template <typename Vector>
 inline void squashPhiGeneral(Vector& v)
 {

buildtest.cc

@@ -50,12 +50,14 @@ void processCandidates(const BinInfoMap & segmentMap, Event& ev, candvec& candid
     }
     if (tmp_candidates.size()>static_cast<size_t>(Config::maxCandsPerSeed)) {
       dprint("huge size=" << tmp_candidates.size() << " keeping best "<< Config::maxCandsPerSeed << " only");
-      std::partial_sort(tmp_candidates.begin(),tmp_candidates.begin()+(Config::maxCandsPerSeed-1),tmp_candidates.end(),sortByHitsChi2);
+      //std::partial_sort(tmp_candidates.begin(),tmp_candidates.begin()+(Config::maxCandsPerSeed-1),tmp_candidates.end(),sortByHitsChi2);
+      std::partial_sort(tmp_candidates.begin(),tmp_candidates.begin()+(Config::maxCandsPerSeed-1),tmp_candidates.end(),sortByScoreCand);
       tmp_candidates.resize(Config::maxCandsPerSeed); // thread local, so ok not thread safe
     } else if (tmp_candidates.size()==0) {
       // save the best candidate from the previous iteration and then swap in
       // the empty new candidate list; seed will be skipped on future iterations
-      auto best = std::min_element(candidates.begin(),candidates.end(),sortByHitsChi2);
+      //auto best = std::min_element(candidates.begin(),candidates.end(),sortByHitsChi2);
+      auto best = std::min_element(candidates.begin(),candidates.end(),sortByScoreCand);
 #ifdef TBB
       std::lock_guard<std::mutex> evtguard(evtlock); // should be rare
 #endif
@@ -113,7 +115,8 @@ void buildTracksBySeeds(const BinInfoMap & segmentMap, Event& ev)
     if (cand.size()>0) {
       // only save one track candidate per seed, one with lowest chi2
       //std::partial_sort(cand.begin(),cand.begin()+1,cand.end(),sortByHitsChi2);
-      auto best = std::min_element(cand.begin(),cand.end(),sortByHitsChi2);
+      //auto best = std::min_element(cand.begin(),cand.end(),sortByHitsChi2);
+      auto best = std::min_element(cand.begin(),cand.end(),sortByScoreCand);
       best->setLabel(evt_track_candidates.size());
       dprint("Saving track from seed " << iseed << " label " << best->label() << " hits " 
                                        << best->nFoundHits() << " parameters " << best->parameters() << std::endl);
@@ -166,7 +169,8 @@ void buildTracksByLayers(const BinInfoMap & segmentMap, Event& ev)
     if (cand.size()>0) {
       // only save one track candidate per seed, one with lowest chi2
       //std::partial_sort(cand.begin(),cand.begin()+1,cand.end(),sortByHitsChi2);
-      auto best = std::min_element(cand.begin(),cand.end(),sortByHitsChi2);
+      //auto best = std::min_element(cand.begin(),cand.end(),sortByHitsChi2);
+      auto best = std::min_element(cand.begin(),cand.end(),sortByScoreCand);
       best->setLabel(evt_track_candidates.size());
       evt_track_candidates.push_back(*best);
       ev.candidateTracksExtra_.emplace_back(evt_seeds_extra[iseed].seedID());

mkFit/CandCloner.cc

@@ -7,12 +7,18 @@
 
 namespace
 {
-inline bool sortCandListByHitsChi2(const mkfit::MkFinder::IdxChi2List& cand1,
-                                   const mkfit::MkFinder::IdxChi2List& cand2)
+inline bool sortCandListByHitsChi2(const mkfit::IdxChi2List& cand1,
+                                   const mkfit::IdxChi2List& cand2)
 {
   if (cand1.nhits == cand2.nhits) return cand1.chi2 < cand2.chi2;
   return cand1.nhits > cand2.nhits;
 }
+
+inline bool sortCandListByScore(const mkfit::IdxChi2List& cand1,
+				const mkfit::IdxChi2List& cand2)
+{
+  return mkfit::sortByScoreStruct(cand1, cand2);
+}
 }
 
 namespace mkfit {
@@ -30,7 +36,7 @@ void CandCloner::ProcessSeedRange(int is_beg, int is_end)
   //1) sort the candidates
   for (int is = is_beg; is < is_end; ++is)
   {
-    std::vector<MkFinder::IdxChi2List>& hitsForSeed = m_hits_to_add[is];
+    std::vector<IdxChi2List>& hitsForSeed = m_hits_to_add[is];
 
     std::vector<CombCandidate> &cands = mp_event_of_comb_candidates->m_candidates;
 
@@ -52,7 +58,8 @@ void CandCloner::ProcessSeedRange(int is_beg, int is_end)
     if ( ! hitsForSeed.empty())
     {
       //sort the new hits
-      std::sort(hitsForSeed.begin(), hitsForSeed.end(), sortCandListByHitsChi2);
+      //std::sort(hitsForSeed.begin(), hitsForSeed.end(), sortCandListByHitsChi2);
+      std::sort(hitsForSeed.begin(), hitsForSeed.end(), sortCandListByScore);
 
       int num_hits = std::min((int) hitsForSeed.size(), Config::maxCandsPerSeed);
 
@@ -61,7 +68,7 @@ void CandCloner::ProcessSeedRange(int is_beg, int is_end)
 
       for (int ih = 0; ih < num_hits; ih++)
       {
-        const MkFinder::IdxChi2List &h2a = hitsForSeed[ih];
+        const IdxChi2List &h2a = hitsForSeed[ih];
 
         cv.push_back( cands[ m_start_seed + is ][ h2a.trkIdx ] );
 

mkFit/CandCloner.h

@@ -81,7 +81,7 @@ class CandCloner
     // Do nothing, "secondary" state vars updated when work completed/assigned.
   }
 
-  void add_cand(int idx, const MkFinder::IdxChi2List& cand_info)
+  void add_cand(int idx, const IdxChi2List& cand_info)
   {
     m_hits_to_add[idx].push_back(cand_info);
 
@@ -166,7 +166,7 @@ class CandCloner
   // eventually, protected or private
 
   int  m_idx_max, m_idx_max_prev;
-  std::vector<std::vector<MkFinder::IdxChi2List>> m_hits_to_add;
+  std::vector<std::vector<IdxChi2List>> m_hits_to_add;
 
   EventOfCombCandidates           *mp_event_of_comb_candidates;
   std::vector<std::pair<int,int>> *mp_kalman_update_list;

mkFit/MkBuilder.cc

@@ -165,6 +165,11 @@ namespace
 
     return cand1.nFoundHits() > cand2.nFoundHits();
   }
+
+  bool sortCandByScore(const Track & cand1, const Track & cand2)
+  {
+    return mkfit::sortByScoreCand(cand1,cand2);    
+  }
 }
 
 //------------------------------------------------------------------------------
@@ -1747,7 +1752,8 @@ void MkBuilder::FindTracksStandard()
         for (int is = 0; is < n_seeds; ++is)
         {
           dprint("dump seed n " << is << " with input candidates=" << tmp_cands[is].size());
-          std::sort(tmp_cands[is].begin(), tmp_cands[is].end(), sortCandByHitsChi2);
+          //std::sort(tmp_cands[is].begin(), tmp_cands[is].end(), sortCandByHitsChi2);
+          std::sort(tmp_cands[is].begin(), tmp_cands[is].end(), sortCandByScore);
 
           if (tmp_cands[is].size() > static_cast<size_t>(Config::maxCandsPerSeed))
           {
@@ -1799,7 +1805,8 @@ void MkBuilder::FindTracksStandard()
       {
         std::vector<Track>& finalcands = eoccs[iseed];
         if (finalcands.size() == 0) continue;
-        std::sort(finalcands.begin(), finalcands.end(), sortCandByHitsChi2);
+        //std::sort(finalcands.begin(), finalcands.end(), sortCandByHitsChi2);
+        std::sort(finalcands.begin(), finalcands.end(), sortCandByScore);
       }
     }); // end parallel-for over chunk of seeds within region
   }); // end of parallel-for-each over eta regions
@@ -1976,7 +1983,8 @@ void MkBuilder::find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr,
         dprintf("Extra layer %d start_seed %d, is %d, N=%d  (orig size = %d) -- ",
                 curr_layer, start_seed, is, (int) extra_cands[is].size(), (int) eoccs.m_candidates[start_seed + is].size());
 
-        std::sort(extra_cands[is].begin(), extra_cands[is].end(), sortCandByHitsChi2);
+        //std::sort(extra_cands[is].begin(), extra_cands[is].end(), sortCandByHitsChi2);
+        std::sort(extra_cands[is].begin(), extra_cands[is].end(), sortCandByScore);
 
         auto src_i  = extra_cands[is].begin();
         auto src_e  = extra_cands[is].end();
@@ -1995,7 +2003,8 @@ void MkBuilder::find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr,
         int cnt = 0;
         while (dest_i != dest_e && src_i != src_e)
         {
-          while (dest_i != dest_e && sortCandByHitsChi2(*dest_i, *src_i)) ++dest_i;
+          //while (dest_i != dest_e && sortCandByHitsChi2(*dest_i, *src_i)) ++dest_i;
+          while (dest_i != dest_e && sortCandByScore(*dest_i, *src_i)) ++dest_i;
 
           if (dest_i != dest_e)
           {
@@ -2021,7 +2030,8 @@ void MkBuilder::find_tracks_in_layers(CandCloner &cloner, MkFinder *mkfndr,
   {
     std::vector<Track>& finalcands = eoccs.m_candidates[iseed];
     if (finalcands.size() == 0) continue;
-    std::sort(finalcands.begin(), finalcands.end(), sortCandByHitsChi2);
+    //std::sort(finalcands.begin(), finalcands.end(), sortCandByHitsChi2);
+    std::sort(finalcands.begin(), finalcands.end(), sortCandByScore);
   }
 }
 

mkFit/MkFinder.cc

@@ -71,7 +71,7 @@ void MkFinder::InputTracksAndHitIdx(const std::vector<CombCandidate>     & track
 }
 
 void MkFinder::InputTracksAndHitIdx(const std::vector<CombCandidate>                       & tracks,
-                                    const std::vector<std::pair<int,MkFinder::IdxChi2List>>& idxs,
+                                    const std::vector<std::pair<int,IdxChi2List>>& idxs,
                                     int beg, int end, bool inputProp)
 {
   // Assign track parameters to initial state and copy hit values in.
@@ -901,6 +901,8 @@ void MkFinder::FindCandidatesCloneEngine(const LayerOfHits &layer_of_hits, CandC
           tmpList.trkIdx = CandIdx(itrack, 0, 0);
           tmpList.hitIdx = XHitArr.At(itrack, hit_cnt, 0);
           tmpList.nhits  = NFoundHits(itrack,0,0) + 1;
+          tmpList.nholes  = num_invalid_hits(itrack);
+          tmpList.pt = std::abs(1.0f/Par[iP].At(itrack,3,0));
           tmpList.chi2   = Chi2(itrack, 0, 0) + chi2;
           cloner.add_cand(SeedIdx(itrack, 0, 0) - offset, tmpList);
           // hitsToAdd[SeedIdx(itrack, 0, 0)-offset].push_back(tmpList);
@@ -934,6 +936,8 @@ void MkFinder::FindCandidatesCloneEngine(const LayerOfHits &layer_of_hits, CandC
     tmpList.trkIdx = CandIdx(itrack, 0, 0);
     tmpList.hitIdx = fake_hit_idx;
     tmpList.nhits  = NFoundHits(itrack,0,0);
+    tmpList.nholes  = num_invalid_hits(itrack)+1;
+    tmpList.pt = std::abs(1.0f/Par[iP].At(itrack,3,0));
     tmpList.chi2   = Chi2(itrack, 0, 0);
     cloner.add_cand(SeedIdx(itrack, 0, 0) - offset, tmpList);
     dprint("adding invalid hit " << fake_hit_idx);

mkFit/MkFinder.h

@@ -54,14 +54,6 @@ class MkFinder : public MkBase
   using MPlexHitIdx = Matriplex::Matriplex<int, MPlexHitIdxMax, 1, NN>;
   using MPlexQHoT   = Matriplex::Matriplex<HitOnTrack, 1, 1, NN>;
 
-  struct IdxChi2List
-  {
-    int   trkIdx; // candidate index
-    int   hitIdx; // hit index
-    int   nhits;  // number of hits (used for sorting)
-    float chi2;   // total chi2 (used for sorting)
-  };
-
   //----------------------------------------------------------------------------
 
   MPlexQF    Chi2;

mkFit/buildtestMPlex.cc

@@ -30,6 +30,12 @@ inline bool sortByHitsChi2(const std::pair<Track, TrackState>& cand1,
   return cand1.first.nFoundHits() > cand2.first.nFoundHits();
 }
 
+inline bool sortByScore(const std::pair<Track, TrackState>& cand1,
+			const std::pair<Track, TrackState>& cand2)
+{
+  return sortByScoreCandPair(cand1,cand2);
+}
+
 inline bool sortByPhi(const Hit& hit1, const Hit& hit2)
 {
   return std::atan2(hit1.y(),hit1.x()) < std::atan2(hit2.y(),hit2.x());