Adapt pixel cpe algo to better handle broken clusters

CondFormats/SiPixelTransient/interface/SiPixelTemplate.h

@@ -285,8 +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
   // 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);
+  bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx, bool goodEdgeAlgo = false);
 
   // 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);

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];

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

@@ -116,3 +116,4 @@ The offsets currently in use are:
 * 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.186: Run-3 goodEdgeAlgo CPE

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.186,
+)
 
 #Workflow to enable displacedRegionalStep tracking iteration
 class UpgradeWorkflow_displacedRegional(UpgradeWorkflowTracking):

RecoLocalTracker/Configuration/python/customizeHLT.py

@@ -0,0 +1,16 @@
+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,