small improvement in fake track rejection

Detectors/MUON/MCH/Tracking/README.md

@@ -97,9 +97,11 @@ attached to the new track. Clusters that drive the track parameters outside of a
 are discarded.
 - Improve the tracks: run the smoother to recompute the local chi2 at each cluster, remove the worst cluster if it does
 not pass a stricter chi2 cut, refit the track and repeat the procedure until all clusters pass the cut or one of them
-cannot be removed (the track must contain at least 1 cluster per station), in which case the track is removed.
+cannot be removed without violating the tracking conditions (by default, the track must contain at least 1 cluster per
+station and both chambers fired on station 4 or 5), in which case the track is removed.
 - Remove connected tracks in station 3, 4 and 5. If two tracks share at least one cluster in these stations, remove the
-one with the smallest number of clusters or with the highest chi2 in case of equality, assuming it is a fake track.
+one with the smallest number of fired chambers or with the highest chi2/(ndf-1) in case of equality, assuming it is a 
+fake track.
 
 In all stations, the search for compatible clusters is done in a way to consider every possibilities, i.e. every
 combinations of 1 to 4 clusters, while skipping the already tested combinations. This includes subsets of previously

Detectors/MUON/MCH/Tracking/include/MCHTracking/Track.h

@@ -39,12 +39,12 @@ class Track
   Track(Track&&) = delete;
   Track& operator=(Track&&) = delete;
 
-  /// Return a reference to the track parameters at vertex
-  const TrackParam& getParamAtVertex() const { return mParamAtVertex; }
-
   /// Return the number of attached clusters
   int getNClusters() const { return mParamAtClusters.size(); }
 
+  /// Return the number of degrees of freedom of the track
+  int getNDF() const { return 2 * getNClusters() - 5; }
+
   /// Return a reference to the track parameters at first cluster
   const TrackParam& first() const { return mParamAtClusters.front(); }
   /// Return a reference to the track parameters at last cluster
@@ -72,7 +72,7 @@ class Track
 
   bool isBetter(const Track& track) const;
 
-  void tagRemovableClusters(uint8_t requestedStationMask);
+  void tagRemovableClusters(uint8_t requestedStationMask, bool request2ChInSameSt45);
 
   void setCurrentParam(const TrackParam& param, int chamber);
   TrackParam& getCurrentParam();
@@ -98,7 +98,6 @@ class Track
   void print() const;
 
  private:
-  TrackParam mParamAtVertex{};                 ///< track parameters at vertex
   std::list<TrackParam> mParamAtClusters{};    ///< list of track parameters at each cluster
   std::unique_ptr<TrackParam> mCurrentParam{}; ///< current track parameters used during tracking
   int mCurrentChamber = -1;                    ///< current chamber on which the current parameters are given

Detectors/MUON/MCH/Tracking/src/Track.cxx

@@ -28,8 +28,7 @@ using namespace std;
 
 //__________________________________________________________________________
 Track::Track(const Track& track)
-  : mParamAtVertex(track.mParamAtVertex),
-    mParamAtClusters(track.mParamAtClusters),
+  : mParamAtClusters(track.mParamAtClusters),
     mCurrentParam(nullptr),
     mCurrentChamber(-1),
     mConnected(track.mConnected),
@@ -133,63 +132,61 @@ bool Track::isBetter(const Track& track) const
 }
 
 //__________________________________________________________________________
-void Track::tagRemovableClusters(uint8_t requestedStationMask)
+void Track::tagRemovableClusters(uint8_t requestedStationMask, bool request2ChInSameSt45)
 {
-  /// Identify clusters that can be removed from the track,
-  /// with the only requirements to have at least 1 cluster per requested station
-  /// and at least 2 chambers over 4 in stations 4 & 5 that contain cluster(s)
-
-  int previousCh(-1), previousSt(-1), nChHitInSt45(0);
-  TrackParam* previousParam(nullptr);
+  /// Identify clusters that can be removed from the track, with the requirement
+  /// to have enough chambers fired to fulfill the tracking criteria
 
+  // count the number of clusters in each chamber and the number of chambers fired on stations 4 and 5
+  int nClusters[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
   for (auto& param : *this) {
+    ++nClusters[param.getClusterPtr()->getChamberId()];
+  }
+  int nChFiredInSt4 = (nClusters[6] > 0) ? 1 : 0;
+  if (nClusters[7] > 0) {
+    ++nChFiredInSt4;
+  }
+  int nChFiredInSt5 = (nClusters[8] > 0) ? 1 : 0;
+  if (nClusters[9] > 0) {
+    ++nChFiredInSt5;
+  }
+  int nChFiredInSt45 = nChFiredInSt4 + nChFiredInSt5;
 
-    int currentCh = param.getClusterPtr()->getChamberId();
-    int currentSt = currentCh / 2;
-
-    // set the cluster as removable if the station is not requested or if it is not alone in the station
-    if (((1 << currentSt) & requestedStationMask) == 0) {
-      param.setRemovable(true);
-    } else if (currentSt == previousSt) {
-      previousParam->setRemovable(true);
-      param.setRemovable(true);
-    } else {
-      param.setRemovable(false);
-      previousSt = currentSt;
-      previousParam = &param;
-    }
+  bool removable[10] = {false, false, false, false, false, false, false, false, false, false};
 
-    // count the number of chambers in station 4 & 5 that contain cluster(s)
-    if (currentCh > 5 && currentCh != previousCh) {
-      ++nChHitInSt45;
-      previousCh = currentCh;
+  // for station 1, 2 and 3, there must be at least one cluster per requested station
+  for (int iCh = 0; iCh < 6; iCh += 2) {
+    if (nClusters[iCh] + nClusters[iCh + 1] > 1 || (requestedStationMask & (1 << (iCh / 2))) == 0) {
+      removable[iCh] = removable[iCh + 1] = true;
     }
   }
 
-  // if there are less than 3 chambers containing cluster(s) in station 4 & 5
-  if (nChHitInSt45 < 3) {
-
-    previousCh = -1;
-    previousParam = nullptr;
-
-    for (auto itParam = this->rbegin(); itParam != this->rend(); ++itParam) {
-
-      int currentCh = itParam->getClusterPtr()->getChamberId();
-
-      if (currentCh < 6) {
-        break;
-      }
-
-      // set the cluster as not removable unless it is not alone in the chamber
-      if (currentCh == previousCh) {
-        previousParam->setRemovable(true);
-        itParam->setRemovable(true);
-      } else {
-        itParam->setRemovable(false);
-        previousCh = currentCh;
-        previousParam = &*itParam;
-      }
+  // for station 4 and 5, there must be at least one cluster per requested station and
+  // at least 2 chambers fired (on the same station or not depending on the requirement)
+  if (nChFiredInSt45 == 4) {
+    removable[6] = removable[7] = removable[8] = removable[9] = true;
+  } else if (nChFiredInSt45 == 3) {
+    if (nChFiredInSt4 == 2 && request2ChInSameSt45) {
+      removable[6] = (nClusters[6] > 1);
+      removable[7] = (nClusters[7] > 1);
+    } else if (nClusters[6] + nClusters[7] > 1 || (requestedStationMask & 0x8) == 0) {
+      removable[6] = removable[7] = true;
     }
+    if (nChFiredInSt5 == 2 && request2ChInSameSt45) {
+      removable[8] = (nClusters[8] > 1);
+      removable[9] = (nClusters[9] > 1);
+    } else if (nClusters[8] + nClusters[9] > 1 || (requestedStationMask & 0x10) == 0) {
+      removable[8] = removable[9] = true;
+    }
+  } else {
+    for (int iCh = 6; iCh < 10; ++iCh) {
+      removable[iCh] = (nClusters[iCh] > 1);
+    }
+  }
+
+  // tag the removable clusters
+  for (auto& param : *this) {
+    param.setRemovable(removable[param.getClusterPtr()->getChamberId()]);
   }
 }
 

Detectors/MUON/MCH/Tracking/src/TrackFinder.cxx

@@ -1122,7 +1122,7 @@ void TrackFinder::improveTracks()
       }
 
       // Identify removable clusters
-      itTrack->tagRemovableClusters(requestedStationMask());
+      itTrack->tagRemovableClusters(requestedStationMask(), !mMoreCandidates);
 
       // Look for the cluster with the worst local chi2
       double worstLocalChi2(-1.);
@@ -1176,7 +1176,7 @@ void TrackFinder::removeConnectedTracks(int stMin, int stMax)
 {
   /// Find and remove tracks sharing 1 cluster or more in station(s) [stMin, stMax]
   /// For each couple of connected tracks, one removes the one with the smallest
-  /// number of clusters or with the highest chi2 value in case of equality
+  /// number of fired chambers or with the highest chi2/(ndf-1) value in case of equality
 
   if (mTracks.size() < 2) {
     return;
@@ -1186,30 +1186,36 @@ void TrackFinder::removeConnectedTracks(int stMin, int stMax)
   int chMax = 2 * stMax + 1;
   int nPlane = 2 * (chMax - chMin + 1);
 
-  // first loop to fill the array of cluster Ids
-  std::vector<uint32_t> ClIds{};
-  ClIds.resize(nPlane * mTracks.size());
+  // first loop to fill the arrays of cluster Ids and number of fired chambers
+  std::vector<uint32_t> ClIds(nPlane * mTracks.size());
+  std::vector<uint8_t> nFiredCh(mTracks.size());
+  int previousCh(-1);
   int iTrack(0);
   for (auto itTrack = mTracks.begin(); itTrack != mTracks.end(); ++itTrack, ++iTrack) {
     for (auto itParam = itTrack->rbegin(); itParam != itTrack->rend(); ++itParam) {
       int ch = itParam->getClusterPtr()->getChamberId();
-      if (ch > chMax) {
-        continue;
-      } else if (ch < chMin) {
-        break;
+      if (ch != previousCh) {
+        ++nFiredCh[iTrack];
+        previousCh = ch;
+      }
+      if (ch >= chMin && ch <= chMax) {
+        ClIds[nPlane * iTrack + 2 * (ch - chMin) + itParam->getClusterPtr()->getDEId() % 2] = itParam->getClusterPtr()->getUniqueId();
       }
-      ClIds[nPlane * iTrack + 2 * (ch - chMin) + itParam->getClusterPtr()->getDEId() % 2] = itParam->getClusterPtr()->getUniqueId();
     }
   }
 
   // second loop to tag the tracks to remove
+  int iTrack1 = mTracks.size() - 1;
   int iindex = ClIds.size() - 1;
-  for (auto itTrack1 = mTracks.rbegin(); itTrack1 != mTracks.rend(); ++itTrack1, iindex -= nPlane) {
+  for (auto itTrack1 = mTracks.rbegin(); itTrack1 != mTracks.rend(); ++itTrack1, iindex -= nPlane, --iTrack1) {
+    int iTrack2 = iTrack1 - 1;
     int jindex = iindex - nPlane;
-    for (auto itTrack2 = std::next(itTrack1); itTrack2 != mTracks.rend(); ++itTrack2) {
+    for (auto itTrack2 = std::next(itTrack1); itTrack2 != mTracks.rend(); ++itTrack2, --iTrack2) {
       for (int iPlane = nPlane; iPlane > 0; --iPlane) {
         if (ClIds[iindex] > 0 && ClIds[iindex] == ClIds[jindex]) {
-          if (itTrack2->isBetter(*itTrack1)) {
+          if ((nFiredCh[iTrack2] > nFiredCh[iTrack1]) ||
+              ((nFiredCh[iTrack2] == nFiredCh[iTrack1]) &&
+               (itTrack2->first().getTrackChi2() / (itTrack2->getNDF() - 1) < itTrack1->first().getTrackChi2() / (itTrack1->getNDF() - 1)))) {
             itTrack1->connected();
           } else {
             itTrack2->connected();

Detectors/MUON/MCH/Tracking/src/TrackFinderOriginal.cxx

@@ -1221,7 +1221,7 @@ void TrackFinderOriginal::improveTracks()
       }
 
       // Identify removable clusters
-      itTrack->tagRemovableClusters(requestedStationMask());
+      itTrack->tagRemovableClusters(requestedStationMask(), false);
 
       // Look for the cluster with the worst local chi2
       double worstLocalChi2(-1.);