Adapt pixel cpe algo to better handle broken clusters

CondFormats/SiPixelTransient/interface/SiPixelTemplate.h

@@ -285,6 +285,9 @@ class SiPixelTemplate {
   // initialize the rest;
   static void postInit(std::vector<SiPixelTemplateStore>& thePixelTemp_);
 
+  // Interpolate with y Gaussian Parameter interpolation to be used with goodEdge reconstruction algorithm
+  bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx, bool goodEdgeAlgo);
+
   // Interpolate input alpha and beta angles to produce a working template for each individual hit.
   bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx);
 

CondFormats/SiPixelTransient/interface/SiPixelUtils.h

@@ -2,21 +2,23 @@
 #define CondFormats_SiPixelTransient_SiPixelUtils_h
 
 namespace siPixelUtils {
-  float generic_position_formula(int size,                    //!< Size of this projection.
-                                 int q_f,                     //!< Charge in the first pixel.
-                                 int q_l,                     //!< Charge in the last pixel.
-                                 float upper_edge_first_pix,  //!< As the name says.
-                                 float lower_edge_last_pix,   //!< As the name says.
-                                 float lorentz_shift,         //!< L-width
-                                 float theThickness,          //detector thickness
-                                 float cot_angle,             //!< cot of alpha_ or beta_
-                                 float pitch,                 //!< thePitchX or thePitchY
-                                 float pitchfraction_first,   //!< true if the first is big
-                                 float pitchfraction_last,    //!< true if the last is big
-                                 float eff_charge_cut_low,    //!< Use edge if > W_eff (in pix) &&&
-                                 float eff_charge_cut_high,   //!< Use edge if < W_eff (in pix) &&&
-                                 float size_cut               //!< Use edge when size == cuts
-  );
+  float generic_position_formula(
+      int size,                     //!< Size of this projection.
+      int q_f,                      //!< Charge in the first pixel.
+      int q_l,                      //!< Charge in the last pixel.
+      float upper_edge_first_pix,   //!< As the name says.
+      float lower_edge_last_pix,    //!< As the name says.
+      float lorentz_shift,          //!< L-width
+      float theThickness,           //detector thickness
+      float cot_angle,              //!< cot of alpha_ or beta_
+      float pitch,                  //!< thePitchX or thePitchY
+      float pitchfraction_first,    //!< true if the first is big
+      float pitchfraction_last,     //!< true if the last is big
+      float eff_charge_cut_low,     //!< Use edge if > W_eff (in pix) &&&
+      float eff_charge_cut_high,    //!< Use edge if < W_eff (in pix) &&&
+      float size_cut,               //!< Use edge when size == cuts
+      float delta_length_cut = 2.,  //!< if charge len - cls size > this (in pix), use one-sided reco
+      bool goodEdgeAlgo = false);
 }  // namespace siPixelUtils
 
 #endif

CondFormats/SiPixelTransient/src/SiPixelTemplate.cc

@@ -82,6 +82,7 @@
 //  V10.21 - Address runtime issues in pushfile() for gcc 7.X due to using tempfile as char string + misc. cleanup [Petar]
 //  V10.22 - Move templateStore to the heap, fix variable name in pushfile() [Petar]
 //  V10.24 - Add sideload() + associated gymnastics [Petar and Oz]
+//  V10.25 - Restore y-residual Gaussian parameters [Morris]
 
 //  Created by Morris Swartz on 10/27/06.
 //
@@ -1324,15 +1325,18 @@ void SiPixelTemplate::postInit(std::vector<SiPixelTemplateStore>& thePixelTemp_)
 //! \param locBx - (input) the sign of this quantity is used to determine whether to flip cot(alpha/beta)<0 quantities from cot(alpha/beta)>0 (FPix only)
 //!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0
 //!                    for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0
+//! \param goodEdgeAlgo - (input) Flag to turn on the y Gaussian Parameter interpolation to be used with goodEdge reconstruction algorithm
 // ************************************************************************************************************
-bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx) {
+bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx, bool goodEdgeAlgo) {
   // Interpolate for a new set of track angles
 
+#ifndef SI_PIXEL_TEMPLATE_STANDALONE
   //check for nan's
   if (!edm::isFinite(cotalpha) || !edm::isFinite(cotbeta)) {
     success_ = false;
     return success_;
   }
+#endif
 
   // Local variables
   int i, j;
@@ -1553,16 +1557,19 @@ bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float l
     }
 
     for (i = 0; i < 4; ++i) {
-      yavg_[i] = (1.f - yratio_) * thePixelTemp_[index_id_].enty[ilow].yavg[i] +
-                 yratio_ * thePixelTemp_[index_id_].enty[ihigh].yavg[i];
-      if (flip_y_) {
-        yavg_[i] = -yavg_[i];
-      }
       yavg_[i] = (1.f - yratio_) * enty0_->yavg[i] + yratio_ * enty1_->yavg[i];
       if (flip_y_) {
         yavg_[i] = -yavg_[i];
       }
       yrms_[i] = (1.f - yratio_) * enty0_->yrms[i] + yratio_ * enty1_->yrms[i];
+
+      if (goodEdgeAlgo) {  // restore y Gaussian Parameter interpolation
+        ygx0_[i] = (1.f - yratio_) * enty0_->ygx0[i] + yratio_ * enty1_->ygx0[i];
+        if (flip_y_) {
+          ygx0_[i] = -ygx0_[i];
+        }
+        ygsig_[i] = (1.f - yratio_) * enty0_->ygsig[i] + yratio_ * enty1_->ygsig[i];
+      }  //if(goodEdgeAlgo)
       chi2yavg_[i] = (1.f - yratio_) * enty0_->chi2yavg[i] + yratio_ * enty1_->chi2yavg[i];
       chi2ymin_[i] = (1.f - yratio_) * enty0_->chi2ymin[i] + yratio_ * enty1_->chi2ymin[i];
       chi2xavg[i] = (1.f - yratio_) * enty0_->chi2xavg[i] + yratio_ * enty1_->chi2xavg[i];
@@ -1869,6 +1876,23 @@ bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float l
   return success_;
 }  // interpolate
 
+// ************************************************************************************************************
+//! Interpolate input alpha and beta angles to produce a working template for each individual hit.
+//! \param id - (input) index of the template to use
+//! \param cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
+//! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
+//! \param locBz - (input) the sign of this quantity is used to determine whether to flip cot(beta)<0 quantities from cot(beta)>0 (FPix only)
+//!                    for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0
+//!                    for Phase 1 FPix IP-related tracks, see next comment
+//! \param locBx - (input) the sign of this quantity is used to determine whether to flip cot(alpha/beta)<0 quantities from cot(alpha/beta)>0 (FPix only)
+//!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0
+//!                    for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0
+// ************************************************************************************************************
+bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx) {
+  // Interpolate for a new set of track angles, but without y Gaussian Parameter interpolation to be used with goodEdge reconstruction algorithm
+  return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx, false);
+}
+
 // ************************************************************************************************************
 //! Interpolate input alpha and beta angles to produce a working template for each individual hit.
 //! \param id - (input) index of the template to use
@@ -1888,7 +1912,7 @@ bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta) {
   if (cotalpha < 0.f) {
     locBz = -locBx;
   }
-  return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx);
+  return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx, false);
 }
 
 // ************************************************************************************************************
@@ -1903,7 +1927,7 @@ bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float l
 
   // Local variables
   float locBx = 1.f;
-  return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx);
+  return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx, false);
 }
 
 // *************************************************************************************************************************************

CondFormats/SiPixelTransient/src/SiPixelUtils.cc

@@ -13,21 +13,23 @@ namespace siPixelUtils {
   //!  X and Y, in the interest of the simplicity of the code, all parameters
   //!  are passed by the caller.
   //-----------------------------------------------------------------------------
-  float generic_position_formula(int size,                    //!< Size of this projection.
-                                 int q_f,                     //!< Charge in the first pixel.
-                                 int q_l,                     //!< Charge in the last pixel.
-                                 float upper_edge_first_pix,  //!< As the name says.
-                                 float lower_edge_last_pix,   //!< As the name says.
-                                 float lorentz_shift,         //!< L-shift at half thickness
-                                 float theThickness,          //detector thickness
-                                 float cot_angle,             //!< cot of alpha_ or beta_
-                                 float pitch,                 //!< thePitchX or thePitchY
-                                 float pitchfraction_first,
-                                 float pitchfraction_last,
-                                 float eff_charge_cut_low,   //!< Use edge if > w_eff  &&&
-                                 float eff_charge_cut_high,  //!< Use edge if < w_eff  &&&
-                                 float size_cut              //!< Use edge when size == cuts
-  ) {
+  float generic_position_formula(
+      int size,                    //!< Size of this projection.
+      int q_f,                     //!< Charge in the first pixel.
+      int q_l,                     //!< Charge in the last pixel.
+      float upper_edge_first_pix,  //!< As the name says.
+      float lower_edge_last_pix,   //!< As the name says.
+      float lorentz_shift,         //!< L-shift at half thickness
+      float theThickness,          //detector thickness
+      float cot_angle,             //!< cot of alpha_ or beta_
+      float pitch,                 //!< thePitchX or thePitchY
+      float pitchfraction_first,
+      float pitchfraction_last,
+      float eff_charge_cut_low,   //!< Use edge if > w_eff  &&&
+      float eff_charge_cut_high,  //!< Use edge if < w_eff  &&&
+      float size_cut,             //!< Use edge when size == cuts
+      float delta_length_cut,     //!< if charge len - cls size > this (in pix), use one-sided reco
+      bool goodEdgeAlgo) {
     float geom_center = 0.5f * (upper_edge_first_pix + lower_edge_last_pix);
 
     //--- The case of only one pixel in this projection is separate.  Note that
@@ -50,17 +52,32 @@ namespace siPixelUtils {
 
     //--- The `effective' charge width -- particle's path in first and last pixels only
     float w_eff = std::abs(w_pred) - w_inner;
+    float delta = w_eff - 0.5 * sum_of_edge * pitch;
 
     //--- If the observed charge width is inconsistent with the expectations
     //--- based on the track, do *not* use w_pred-w_innner.  Instead, replace
     //--- it with an *average* effective charge width, which is the average
     //--- length of the edge pixels.
     //
-    //  bool usedEdgeAlgo = false;
     if ((size >= size_cut) || ((w_eff / pitch < eff_charge_cut_low) | (w_eff / pitch > eff_charge_cut_high))) {
       w_eff = pitch * 0.5f * sum_of_edge;  // ave. length of edge pixels (first+last) (cm)
-                                           //  usedEdgeAlgo = true;
-    }
+
+      if (goodEdgeAlgo) {
+        if (delta / pitch > delta_length_cut) {
+          //  observed cluster is much shorter than expected, use one-sided reco
+          if (w_pred > 0.f) {
+            // x1,x2 are centers of the first last last pixel coordinates
+            float x1 = upper_edge_first_pix - 0.5 * pitchfraction_first * pitch;
+            float hit_pos = x1 + 0.5 * w_pred;
+            return hit_pos;
+          } else {
+            float x2 = lower_edge_last_pix + 0.5 * pitchfraction_last * pitch;
+            float hit_pos = x2 + 0.5 * w_pred;
+            return hit_pos;
+          }
+        }  //if (delta / pitch > delta_length_cut)
+      }  //if(goodEdgeAlgo)
+    }  //if(size >= size_cut) ||...
 
     //--- Finally, compute the position in this projection
     float q_diff = q_l - q_f;

Configuration/ProcessModifiers/python/siPixelGoodEdgeAlgo_cff.py

@@ -0,0 +1,4 @@
+import FWCore.ParameterSet.Config as cms
+
+# This modifier enables the good edge algorithm in pixel hit reconstruction that handles broken/truncated pixel cluster caused by radiation damage
+siPixelGoodEdgeAlgo = cms.Modifier()

Configuration/PyReleaseValidation/README.md

@@ -115,4 +115,4 @@ The offsets currently in use are:
 * 0.112: Activate OuterTracker inefficiency (PS-p: bias rails inefficiency; PS-s and SS: 1% bad strips)
 * 0.113: Activate OuterTracker inefficiency (PS-p: bias rails inefficiency; PS-s and SS: 5% bad strips)
 * 0.114: Activate OuterTracker inefficiency (PS-p: bias rails inefficiency; PS-s and SS: 10% bad strips)
-* 0.141: Activate emulation of the signal shape of the InnerTracker FE chip (CROC)
+* 0.141: Activate emulation of the signal shape of the InnerTracker FE chip (CROC)

Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py

@@ -640,6 +640,40 @@ def condition(self, fragment, stepList, key, hasHarvest):
     offset = 0.18,
 )
 
+# pixel GoodEdgeAlgo CPE workflows
+class UpgradeWorkflow_siPixelGoodEdgeAlgo(UpgradeWorkflow):
+    def setup_(self, step, stepName, stepDict, k, properties):
+        if 'HLTOnly'==step:
+            stepDict[stepName][k] = merge([{'--customise' : 'RecoLocalTracker/Configuration/customizeHLT.customiseFor47966'}, stepDict[step][k]])
+        if 'Reco' in step:
+            stepDict[stepName][k] = merge([{'--procModifiers': 'siPixelGoodEdgeAlgo'}, stepDict[step][k]])
+    def condition(self, fragment, stepList, key, hasHarvest):
+        result = (fragment=="QCD_Pt_1800_2400_14" or fragment=="TTbar_14TeV" ) and any(y in key for y in ['2025'])
+        return result
+upgradeWFs['siPixelGoodEdgeAlgo'] = UpgradeWorkflow_siPixelGoodEdgeAlgo(
+    steps = [
+        'HLTOnly',
+        'DigiTrigger',
+        'Reco',
+        'RecoFakeHLT',
+        'RecoGlobal',
+        'RecoGlobalFakeHLT',
+        'RecoNano',
+        'RecoNanoFakeHLT',
+    ],
+    PU = [
+        'HLTOnly',
+        'DigiTrigger',
+        'Reco',
+        'RecoFakeHLT',
+        'RecoGlobal',
+        'RecoGlobalFakeHLT',
+        'RecoNano',
+        'RecoNanoFakeHLT',
+    ],
+    suffix = '_siPixelGoodEdgeAlgo',
+    offset = 0.181,
+)
 
 #Workflow to enable displacedRegionalStep tracking iteration
 class UpgradeWorkflow_displacedRegional(UpgradeWorkflowTracking):
@@ -3705,4 +3739,4 @@ def __init__(self, howMuch, dataset):
     ('Hydjet_Quenched_MinBias_5362GeV_cfi', UpgradeFragment(U2000by1,'HydjetQMinBias_5362GeV')),
     ('Hydjet_Quenched_MinBias_5519GeV_cfi', UpgradeFragment(U2000by1,'HydjetQMinBias_5519GeV')),
     ('SingleMuPt15Eta0_0p4_cfi', UpgradeFragment(Kby(9,100),'SingleMuPt15Eta0p_0p4')),
-])
+])

RecoLocalTracker/Configuration/python/customizeHLT.py

@@ -0,0 +1,15 @@
+import FWCore.ParameterSet.Config as cms
+
+# helper functions
+from HLTrigger.Configuration.common import *
+
+
+# Enable good edge CPE algorithm in pixel local reconstruction
+def customiseFor47966(process):
+    for prod in esproducers_by_type(process, 'PixelCPEGenericESProducer', 'PixelCPEFastParamsESProducerAlpakaPhase1@alpaka'):
+        if not hasattr(prod, 'GoodEdgeAlgo'):
+            prod.GoodEdgeAlgo = cms.bool(True)
+        else:
+            prod.GoodEdgeAlgo = True
+
+    return process 

RecoLocalTracker/SiPixelRecHits/interface/PixelCPEClusterRepair.h

@@ -102,6 +102,7 @@ class PixelCPEClusterRepair : public PixelCPEBase {
   std::vector<SiPixelTemplateStore2D> thePixelTemp2D_;
 
   int speed_;
+  bool goodEdgeAlgo_;
 
   bool UseClusterSplitter_;
 

RecoLocalTracker/SiPixelRecHits/interface/PixelCPEFastParamsHost.h

@@ -26,7 +26,8 @@ class PixelCPEFastParamsHost : public PixelCPEGenericBase {
                          const TrackerTopology& ttopo,
                          const SiPixelLorentzAngle* lorentzAngle,
                          const SiPixelGenErrorDBObject* genErrorDBObject,
-                         const SiPixelLorentzAngle* lorentzAngleWidth);
+                         const SiPixelLorentzAngle* lorentzAngleWidth,
+                         const bool goodEdgeAlgo);
 
   // non-copyable
   PixelCPEFastParamsHost(PixelCPEFastParamsHost const&) = delete;
@@ -61,6 +62,7 @@ class PixelCPEFastParamsHost : public PixelCPEGenericBase {
   void fillParamsForDevice();
 
   Buffer buffer_;
+  bool goodEdgeAlgo_;
 };
 
 #endif  // RecoLocalTracker_SiPixelRecHits_interface_PixelCPEFastParamsHost_h

RecoLocalTracker/SiPixelRecHits/interface/PixelCPEGeneric.h

@@ -81,12 +81,14 @@ class PixelCPEGeneric : public PixelCPEGenericBase {
   float the_eff_charge_cut_highY;
   float the_size_cutX;
   float the_size_cutY;
+  float delta_length_cut;
 
   bool inflate_errors;
   bool inflate_all_errors_no_trk_angle;
 
   bool DoCosmics_;
   bool IrradiationBiasCorrection_;
+  bool goodEdgeAlgo_;
   bool isPhase2_;
   bool NoTemplateErrorsWhenNoTrkAngles_;
 

RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateReco.h

@@ -79,6 +79,7 @@ class PixelCPETemplateReco : public PixelCPEBase {
   const std::vector<SiPixelTemplateStore> *thePixelTemp_;
 
   int speed_;
+  bool goodEdgeAlgo_;
 
   bool UseClusterSplitter_;
 

RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateReco.h

@@ -99,7 +99,26 @@ namespace SiPixelTemplateReco {
                       std::vector<std::pair<int, int> >& zeropix,
                       float& probQ,
                       int& nypix,
-                      int& nxpix);
+                      int& nxpix,
+                      bool goodEdgeAlgo);
+
+  int PixelTempReco1D(int id,
+                      float cotalpha,
+                      float cotbeta,
+                      float locBz,
+                      float locBx,
+                      ClusMatrix& cluster,
+                      SiPixelTemplate& templ,
+                      float& yrec,
+                      float& sigmay,
+                      float& proby,
+                      float& xrec,
+                      float& sigmax,
+                      float& probx,
+                      int& qbin,
+                      int speed,
+                      float& probQ,
+                      bool goodEdgeAlgo);
 
   int PixelTempReco1D(int id,
                       float cotalpha,