From 5a3b3bdeffc04565198214d280cf98e588174757 Mon Sep 17 00:00:00 2001
From: Shigeki Hirose <shigeki.hirose@cern.ch>
Date: Mon, 25 May 2020 13:52:44 +0000
Subject: [PATCH] Fix distCut value and track selection

---
 .../SCT_Monitoring/src/SCTHitEffMonAlg.cxx    | 84 ++++++++++++++++++-
 .../SCT_Monitoring/src/SCTHitEffMonAlg.h      |  6 +-
 2 files changed, 85 insertions(+), 5 deletions(-)

diff --git a/InnerDetector/InDetMonitoring/SCT_Monitoring/src/SCTHitEffMonAlg.cxx b/InnerDetector/InDetMonitoring/SCT_Monitoring/src/SCTHitEffMonAlg.cxx
index bf887f5fd3a..cc8c861360b 100644
--- a/InnerDetector/InDetMonitoring/SCT_Monitoring/src/SCTHitEffMonAlg.cxx
+++ b/InnerDetector/InDetMonitoring/SCT_Monitoring/src/SCTHitEffMonAlg.cxx
@@ -148,6 +148,24 @@ int SCTHitEffMonAlg::becIdxLayer2Index(const int becIdx, const int layer) const
   }
 }
 
+int SCTHitEffMonAlg::getWaferIndex(const int barrel_ec, const int layer_disk, const int side) const {
+  int waferIndex = -1;
+  if (barrel_ec == BARREL) {
+    // corresponds to the waferIndex of B3 side0
+    waferIndex = 0; 
+  } else if (barrel_ec == ENDCAP_A) {
+    // corresponds to the waferIndex of EA0 side0
+    waferIndex = N_BARRELS*N_SIDES;
+  } else if (barrel_ec == ENDCAP_C) {
+    // corresponds to the waferIndex of EC0 side0
+    waferIndex = N_BARRELS*N_SIDES + N_ENDCAPS*N_SIDES;
+  } else {
+    ATH_MSG_WARNING("The barrel_bc index" << barrel_ec << " is not defined.");
+    return waferIndex;
+  }
+  return waferIndex + layer_disk * N_SIDES + side;
+}
+
 double SCTHitEffMonAlg::getResidual(const Identifier& surfaceID,
                                     const Trk::TrackParameters* trkParam,
                                     const InDet::SCT_ClusterContainer* p_sctclcontainer) const {
@@ -346,6 +364,8 @@ StatusCode SCTHitEffMonAlg::fillHistograms(const EventContext& ctx) const {
 
   // Loop over original track collection
   for (const Trk::Track* pthisTrack: *tracks) {
+
+    // First, go through all cuts in this block
     ATH_MSG_VERBOSE("Starting new track");
     if (pthisTrack==nullptr) {
       continue;
@@ -395,7 +415,17 @@ StatusCode SCTHitEffMonAlg::fillHistograms(const EventContext& ctx) const {
       0, 0, 0
     };
     int pixelNHits{0};
+    int pixelNHoles{0};
     int trtNHits{0};
+
+    int sctNHitsPerRegion[N_LAYERS_TOTAL*N_SIDES] = {0};
+    int sctNHolesPerRegion[N_LAYERS_TOTAL*N_SIDES] = {0};
+    // Above two variables hold the number of hits for each SCT disk / layer.
+    // [N_LAYERS_TOTAL*N_SIDES(= 44)] indicates the waferIndex defined as below.
+    //  0- 7: B3 side0, B3 side1, B4 side0, ... B6 side1
+    //  8-25: EA0 side0, EA1 side1, ... EA8 side1
+    // 26-43: EC0 side0, EC1 side1, ... EC8 side1
+
     std::map < Identifier, double > mapOfTrackHitResiduals;
     double zmin = std::numeric_limits<float>::max();
     double zmax = -std::numeric_limits<float>::max();
@@ -405,13 +435,23 @@ StatusCode SCTHitEffMonAlg::fillHistograms(const EventContext& ctx) const {
     float max_layerSide{-1.};
     Identifier surfaceID;
 
-    // loop over all hits on track
+    // Loop over all TSOS (track state on surface) on track to check number of hits / holes on Pixel, SCT and TRT
     for (const Trk::TrackStateOnSurface* tsos: *(trackWithHoles->trackStateOnSurfaces())) {
       surfaceID = surfaceOnTrackIdentifier(tsos);
 
       if (not surfaceID.is_valid()) {
         continue;
       }
+      
+      // Check waferIndex; if the default value of -1 is kept, the corresponding TSOS is not associated with SCT.
+      int waferIndex = -1;
+      // Calculate waferIndex
+      if (m_sctId->is_sct(surfaceID)) {
+        waferIndex = getWaferIndex(m_sctId->barrel_ec(surfaceID),
+                                   m_sctId->layer_disk(surfaceID),
+                                   m_sctId->side(surfaceID));
+      }
+      
       if (tsos->type(Trk::TrackStateOnSurface::Measurement) or tsos->type(Trk::TrackStateOnSurface::Outlier)) {
         if (m_pixelId->is_pixel(surfaceID)) {
           pixelNHits++;
@@ -422,6 +462,13 @@ StatusCode SCTHitEffMonAlg::fillHistograms(const EventContext& ctx) const {
         if (m_sctId->is_sct(surfaceID)) {
           NHits[bec2Index(m_sctId->barrel_ec(surfaceID))]++;
           mapOfTrackHitResiduals[surfaceID] = getResidual(surfaceID, tsos->trackParameters(), &*p_sctclcontainer);
+          sctNHitsPerRegion[waferIndex]++;
+        }
+      } else if (tsos->type(Trk::TrackStateOnSurface::Hole)) {
+        if (m_pixelId->is_pixel(surfaceID)) {
+          pixelNHoles++;
+        } else if (m_sctId->is_sct(surfaceID)) {
+          sctNHolesPerRegion[waferIndex]++;
         }
       }
 
@@ -460,6 +507,7 @@ StatusCode SCTHitEffMonAlg::fillHistograms(const EventContext& ctx) const {
       layersCrossedByTrack[i].resize(n_layers[i] * 2, false);
     }
 
+    // Loop over all TSOS again; this time, to extract SCT-related hits and holes.
     for (const Trk::TrackStateOnSurface* tsos: *(trackWithHoles->trackStateOnSurfaces())) {
       ATH_MSG_VERBOSE("Starting new hit");
       surfaceID = surfaceOnTrackIdentifier(tsos);
@@ -468,11 +516,38 @@ StatusCode SCTHitEffMonAlg::fillHistograms(const EventContext& ctx) const {
         continue;
       }
 
-      unsigned int isub{bec2Index(m_sctId->barrel_ec(surfaceID))};
-      ATH_MSG_VERBOSE("New SCT candidate: " << m_sctId->print_to_string(surfaceID));
 
       int side{m_sctId->side(surfaceID)};
       int layer{m_sctId->layer_disk(surfaceID)};
+      int bec{m_sctId->barrel_ec(surfaceID)};
+      unsigned int isub{bec2Index(bec)};
+      ATH_MSG_VERBOSE("New SCT candidate: " << m_sctId->print_to_string(surfaceID));
+
+      int waferIndex = getWaferIndex(bec, layer, side);
+
+      Int_t sctNHitsExceptThisWafer{0};
+      Int_t sctNHolesExceptThisWafer{0};
+
+      for (Int_t i=0; i<N_LAYERS_TOTAL*N_SIDES; i++) {
+        if (i != waferIndex) {
+          sctNHitsExceptThisWafer  += sctNHitsPerRegion[i];
+          sctNHolesExceptThisWafer += sctNHolesPerRegion[i];
+        }
+      }
+
+      // The track is required to satisfy:
+      // - Number of Si hits to be >= 8
+      // - Number of Si holes to be <= 1
+      // without counting on this TSOS object. (avoid tracking bias.)
+      if ((unsigned int)(sctNHitsExceptThisWafer + pixelNHits) < m_minSiHits) {
+        ATH_MSG_VERBOSE("This track is rejected due to the number of hits: " << sctNHitsExceptThisWafer * pixelNHits);
+        continue;
+      }
+      if ((unsigned int)(sctNHolesExceptThisWafer + pixelNHoles) > m_maxSiHoles) {
+        ATH_MSG_VERBOSE("This track is rejected due to the number of holes: " << sctNHolesExceptThisWafer * pixelNHoles);
+        continue;
+      }
+
       std::string etaPhiSuffix = "_" + std::to_string(layer) + "_" + std::to_string(side);
       const int detIndex{becIdxLayer2Index(isub, layer)};
       if (detIndex == -1) {
@@ -486,7 +561,7 @@ StatusCode SCTHitEffMonAlg::fillHistograms(const EventContext& ctx) const {
       float dedicated_layerPlusHalfSide{static_cast<float>(layer) + static_cast<float>((side + 1) % 2) * 0.5f};
       const Trk::TrackParameters* trkParamOnSurface{tsos->trackParameters()};
       double trackHitResidual{getResidual(surfaceID, trkParamOnSurface, &*p_sctclcontainer)};
-
+      
       float distCut{m_effdistcut};
 
       if (tsos->type(Trk::TrackStateOnSurface::Measurement) or tsos->type(Trk::TrackStateOnSurface::Outlier)) {
@@ -679,3 +754,4 @@ StatusCode SCTHitEffMonAlg::fillHistograms(const EventContext& ctx) const {
   
   return StatusCode::SUCCESS;
 }
+
diff --git a/InnerDetector/InDetMonitoring/SCT_Monitoring/src/SCTHitEffMonAlg.h b/InnerDetector/InDetMonitoring/SCT_Monitoring/src/SCTHitEffMonAlg.h
index 8d0308676a0..2b58a1289a5 100644
--- a/InnerDetector/InDetMonitoring/SCT_Monitoring/src/SCTHitEffMonAlg.h
+++ b/InnerDetector/InDetMonitoring/SCT_Monitoring/src/SCTHitEffMonAlg.h
@@ -68,6 +68,7 @@ class SCTHitEffMonAlg : public AthMonitorAlgorithm {
   ///Convert a layer/disk number (0-21) to a layer number (0-8 for endcaps, 0-3 for barrel)
   int layerIndex2layer(const int index) const;
   int becIdxLayer2Index(const int becIdx, const int layer) const;
+  int getWaferIndex(const int barrel_bc, const int layer_disk, const int side) const;
 
   std::string m_path;
 
@@ -99,9 +100,11 @@ class SCTHitEffMonAlg : public AthMonitorAlgorithm {
   FloatProperty m_maxChi2{this, "MaxChi2", 3.};
   FloatProperty m_maxD0{this, "Maxd0", 10., "mm of D0"};
   FloatProperty m_minPt{this, "MinPt", 1000., "minimu pt in MeV/c"};
-  FloatProperty m_effdistcut{this, "effDistanceCut", 2., "mm"};
+  FloatProperty m_effdistcut{this, "effDistanceCut", 0.2, "mm"};
   FloatProperty m_maxZ0sinTheta{this, "MaxZ0sinTheta", 0.};
   UnsignedIntegerProperty m_maxTracks{this, "MaxTracks", 500};
+  UnsignedIntegerProperty m_minSiHits{this, "MinimumNumberOfSiHits", 8, "Threshold for number of Si hits. Count Si hits excluding hits in the wafer under investigation to reduce track selection bias"};
+  UnsignedIntegerProperty m_maxSiHoles{this, "MaximumNumberOfSiHoles", 1, "Threshold for number of Si holes. Count Si holes excluding holes in the wafer under investigation to reduce track selection bias"};
 
   BooleanProperty m_insideOutOnly{this, "InsideOutOnly", false};
   BooleanProperty m_isCosmic{this, "IsCosmic", false};
@@ -112,6 +115,7 @@ class SCTHitEffMonAlg : public AthMonitorAlgorithm {
   BooleanProperty m_requireGuardRing{this, "RequireGuardRing", false, "should be returned to true"};
   BooleanProperty m_vetoBadChips{this, "VetoBadChips", true};
   BooleanProperty m_useIDGlobal{this, "useIDGlobal", false};  
+  
 };
 
 #endif // SCTHITEFFMONALG_H
-- 
GitLab