SiPixelRecHits/src/PixelCPEClusterRepair.cc

0001 // Include our own header first
0002 #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEClusterRepair.h"
0003
0004 // Geometry services
0005 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0006 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0007
0008 // MessageLogger
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 // Magnetic field
0012 #include "MagneticField/Engine/interface/MagneticField.h"
0013
0014 // Commented for now (3/10/17) until we figure out how to resuscitate 2D template splitter
0015 /// #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateSplit.h"
0016
0017 #include <vector>
0018 #include "boost/multi_array.hpp"
0019 #include <boost/regex.hpp>
0020 #include <map>
0021
0022 #include <iostream>
0023
0024 using namespace SiPixelTemplateReco;
0025 //using namespace SiPixelTemplateSplit;
0026 using namespace std;
0027
0028 namespace {
0029   constexpr float micronsToCm = 1.0e-4;
0030   constexpr int cluster_matrix_size_x = 13;
0031   constexpr int cluster_matrix_size_y = 21;
0032 }  // namespace
0033
0034 //-----------------------------------------------------------------------------
0035 //  Constructor.
0036 //
0037 //-----------------------------------------------------------------------------
0038 PixelCPEClusterRepair::PixelCPEClusterRepair(edm::ParameterSet const& conf,
0039                                              const MagneticField* mag,
0040                                              const TrackerGeometry& geom,
0041                                              const TrackerTopology& ttopo,
0042                                              const SiPixelLorentzAngle* lorentzAngle,
0043                                              const SiPixelTemplateDBObject* templateDBobject,
0044                                              const SiPixel2DTemplateDBObject* templateDBobject2D)
0045     : PixelCPEBase(conf, mag, geom, ttopo, lorentzAngle, nullptr, templateDBobject, nullptr, 1) {
0046   LogDebug("PixelCPEClusterRepair::(constructor)") << endl;
0047
0048   //--- Parameter to decide between DB or text file template access
0049   if (LoadTemplatesFromDB_) {
0050     // Initialize template store to the selected ID [Morris, 6/25/08]
0051     if (!SiPixelTemplate::pushfile(*templateDBobject_, thePixelTemp_))
0052       throw cms::Exception("PixelCPEClusterRepair")
0053           << "\nERROR: Templates not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject version "
0054           << (*templateDBobject_).version() << "\n\n";
0055
0056     // Initialize template store to the selected ID [Morris, 6/25/08]
0057     if (!SiPixelTemplate2D::pushfile(*templateDBobject2D, thePixelTemp2D_))
0058       throw cms::Exception("PixelCPEClusterRepair")
0059           << "\nERROR: Templates not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject2D version "
0060           << (*templateDBobject2D).version() << "\n\n";
0061   } else {
0062     LogDebug("PixelCPEClusterRepair") << "Loading templates for barrel and forward from ASCII files." << endl;
0063     //--- (Archaic) Get configurable template IDs.  This code executes only if we loading pixels from ASCII
0064     //    files, and then they are mandatory.
0065     barrelTemplateID_ = conf.getParameter<int>("barrelTemplateID");
0066     forwardTemplateID_ = conf.getParameter<int>("forwardTemplateID");
0067     templateDir_ = conf.getParameter<int>("directoryWithTemplates");
0068
0069     if (!SiPixelTemplate::pushfile(barrelTemplateID_, thePixelTemp_, templateDir_))
0070       throw cms::Exception("PixelCPEClusterRepair")
0071           << "\nERROR: Template ID " << barrelTemplateID_
0072           << " not loaded correctly from text file. Reconstruction will fail.\n\n";
0073
0074     if (!SiPixelTemplate::pushfile(forwardTemplateID_, thePixelTemp_, templateDir_))
0075       throw cms::Exception("PixelCPEClusterRepair")
0076           << "\nERROR: Template ID " << forwardTemplateID_
0077           << " not loaded correctly from text file. Reconstruction will fail.\n\n";
0078   }
0079
0080   speed_ = conf.getParameter<int>("speed");
0081   LogDebug("PixelCPEClusterRepair::PixelCPEClusterRepair:") << "Template speed = " << speed_ << "\n";
0082
0083   // this returns the magnetic field value in kgauss (1T = 10 kgauss)
0084   int theMagField = mag->nominalValue();
0085
0086   if (theMagField >= 36 && theMagField < 39) {
0087     LogDebug("PixelCPEClusterRepair::PixelCPEClusterRepair:")
0088         << "Magnetic field value is: " << theMagField << " kgauss. Algorithm is being run \n";
0089
0090     templateDBobject2D_ = templateDBobject2D;
0091     fill2DTemplIDs();
0092   }
0093
0094   UseClusterSplitter_ = conf.getParameter<bool>("UseClusterSplitter");
0095
0096   maxSizeMismatchInY_ = conf.getParameter<double>("MaxSizeMismatchInY");
0097   minChargeRatio_ = conf.getParameter<double>("MinChargeRatio");
0098
0099   // read sub-detectors and apply rule to recommend 2D
0100   // can be:
0101   //     XYX (XYZ = PXB, PXE)
0102   //     XYZ n (XYZ as above, n = layer, wheel or disk = 1 .. 6 ;)
0103   std::vector<std::string> str_recommend2D = conf.getParameter<std::vector<std::string>>("Recommend2D");
0104   recommend2D_.reserve(str_recommend2D.size());
0105   for (auto& str : str_recommend2D) {
0106     recommend2D_.push_back(str);
0107   }
0108
0109   // do not recommend 2D if theMagField!=3.8T
0110   if (theMagField < 36 || theMagField > 39) {
0111     recommend2D_.clear();
0112   }
0113
0114   // run CR on damaged clusters (and not only on edge hits)
0115   runDamagedClusters_ = conf.getParameter<bool>("RunDamagedClusters");
0116 }
0117
0118 //-----------------------------------------------------------------------------
0119 //  Fill all 2D template IDs
0120 //-----------------------------------------------------------------------------
0121 void PixelCPEClusterRepair::fill2DTemplIDs() {
0122   auto const& dus = geom_.detUnits();
0123   unsigned m_detectors = dus.size();
0124   for (unsigned int i = 1; i < 7; ++i) {
0125     LogDebug("PixelCPEClusterRepair:: LookingForFirstStrip")
0126         << "Subdetector " << i << " GeomDetEnumerator " << GeomDetEnumerators::tkDetEnum[i] << " offset "
0127         << geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) << " is it strip? "
0128         << (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) != dus.size()
0129                 ? dus[geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i])]->type().isOuterTracker()
0130                 : false);
0131     if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) != dus.size() &&
0132         dus[geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i])]->type().isOuterTracker()) {
0133       if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) < m_detectors)
0134         m_detectors = geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]);
0135     }
0136   }
0137   LogDebug("LookingForFirstStrip") << " Chosen offset: " << m_detectors;
0138
0139   m_DetParams.resize(m_detectors);
0140   LogDebug("PixelCPEClusterRepair::fillDetParams():") << "caching " << m_detectors << " pixel detectors" << endl;
0141   //Loop through detectors and store 2D template ID
0142   bool printed_info = false;
0143   for (unsigned i = 0; i != m_detectors; ++i) {
0144     auto& p = m_DetParams[i];
0145
0146     p.detTemplateId2D = templateDBobject2D_->getTemplateID(p.theDet->geographicalId());
0147     if (p.detTemplateId != p.detTemplateId2D && !printed_info) {
0148       edm::LogWarning("PixelCPEClusterRepair")
0149           << "different template ID between 1D and 2D " << p.detTemplateId << " " << p.detTemplateId2D << endl;
0150       printed_info = true;
0151     }
0152   }
0153 }
0154
0155 //-----------------------------------------------------------------------------
0156 //  Clean up.
0157 //-----------------------------------------------------------------------------
0158 PixelCPEClusterRepair::~PixelCPEClusterRepair() {}
0159
0160 std::unique_ptr<PixelCPEBase::ClusterParam> PixelCPEClusterRepair::createClusterParam(const SiPixelCluster& cl) const {
0161   return std::make_unique<ClusterParamTemplate>(cl);
0162 }
0163
0164 //------------------------------------------------------------------
0165 //  Public methods mandated by the base class.
0166 //------------------------------------------------------------------
0167
0168 //------------------------------------------------------------------
0169 //  The main call to the template code.
0170 //------------------------------------------------------------------
0171 LocalPoint PixelCPEClusterRepair::localPosition(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0172   ClusterParamTemplate& theClusterParam = static_cast<ClusterParamTemplate&>(theClusterParamBase);
0173   bool filled_from_2d = false;
0174
0175   if (!GeomDetEnumerators::isTrackerPixel(theDetParam.thePart))
0176     throw cms::Exception("PixelCPEClusterRepair::localPosition :") << "A non-pixel detector type in here?";
0177
0178   int ID1 = -9999;
0179   int ID2 = -9999;
0180   if (LoadTemplatesFromDB_) {
0181     ID1 = theDetParam.detTemplateId;
0182     ID2 = theDetParam.detTemplateId2D;
0183   } else {  // from asci file
0184     if (!GeomDetEnumerators::isEndcap(theDetParam.thePart))
0185       ID1 = ID2 = barrelTemplateID_;  // barrel
0186     else
0187       ID1 = ID2 = forwardTemplateID_;  // forward
0188   }
0189
0190   // &&& PM, note for later: PixelCPEBase calculates minInX,Y, and maxInX,Y
0191   //     Why can't we simply use that and save time with row_offset, col_offset
0192   //     and mrow = maxInX-minInX, mcol = maxInY-minInY ... Except that we
0193   //     also need to take into account cluster_matrix_size_x,y.
0194
0195   //--- Preparing to retrieve ADC counts from the SiPixeltheClusterParam.theCluster
0196   //    The pixels from minPixelRow() will go into clust_array_2d[0][*],
0197   //    and the pixels from minPixelCol() will go into clust_array_2d[*][0].
0198   int row_offset = theClusterParam.theCluster->minPixelRow();
0199   int col_offset = theClusterParam.theCluster->minPixelCol();
0200
0201   //--- Store the coordinates of the center of the (0,0) pixel of the array that
0202   //    gets passed to PixelTempReco2D.  Will add these values to the output of TemplReco2D
0203   float tmp_x = float(row_offset) + 0.5f;
0204   float tmp_y = float(col_offset) + 0.5f;
0205
0206   //--- Store these offsets (to be added later) in a LocalPoint after tranforming
0207   //    them from measurement units (pixel units) to local coordinates (cm)
0208   LocalPoint lp;
0209   if (theClusterParam.with_track_angle)
0210     //--- Update call with trk angles needed for bow/kink corrections  // Gavril
0211     lp = theDetParam.theTopol->localPosition(MeasurementPoint(tmp_x, tmp_y), theClusterParam.loc_trk_pred);
0212   else {
0213     edm::LogError("PixelCPEClusterRepair") << "@SUB = PixelCPEClusterRepair::localPosition"
0214                                            << "Should never be here. PixelCPEClusterRepair should always be called "
0215                                               "with track angles. This is a bad error !!! ";
0216     lp = theDetParam.theTopol->localPosition(MeasurementPoint(tmp_x, tmp_y));
0217   }
0218
0219   //--- Compute the size of the matrix which will be passed to TemplateReco.
0220   //    We'll later make  clustMatrix[ mrow ][ mcol ]
0221   int mrow = 0, mcol = 0;
0222   for (int i = 0; i != theClusterParam.theCluster->size(); ++i) {
0223     auto pix = theClusterParam.theCluster->pixel(i);
0224     int irow = int(pix.x);
0225     int icol = int(pix.y);
0226     mrow = std::max(mrow, irow);
0227     mcol = std::max(mcol, icol);
0228   }
0229   mrow -= row_offset;
0230   mrow += 1;
0231   mrow = std::min(mrow, cluster_matrix_size_x);
0232   mcol -= col_offset;
0233   mcol += 1;
0234   mcol = std::min(mcol, cluster_matrix_size_y);
0235   assert(mrow > 0);
0236   assert(mcol > 0);
0237
0238   //--- Make and fill the bool arrays flagging double pixels
0239   bool xdouble[mrow], ydouble[mcol];
0240   // x directions (shorter), rows
0241   for (int irow = 0; irow < mrow; ++irow)
0242     xdouble[irow] = theDetParam.theRecTopol->isItBigPixelInX(irow + row_offset);
0243   //
0244   // y directions (longer), columns
0245   for (int icol = 0; icol < mcol; ++icol)
0246     ydouble[icol] = theDetParam.theRecTopol->isItBigPixelInY(icol + col_offset);
0247
0248   //--- C-style matrix.  We'll need it in either case.
0249   float clustMatrix[mrow][mcol];
0250   float clustMatrix2[mrow][mcol];
0251
0252   //--- Prepare struct that passes pointers to TemplateReco.  It doesn't own anything.
0253   SiPixelTemplateReco::ClusMatrix clusterPayload{&clustMatrix[0][0], xdouble, ydouble, mrow, mcol};
0254   SiPixelTemplateReco2D::ClusMatrix clusterPayload2d{&clustMatrix2[0][0], xdouble, ydouble, mrow, mcol};
0255
0256   //--- Copy clust's pixels (calibrated in electrons) into clustMatrix;
0257   memset(clustMatrix, 0, sizeof(float) * mrow * mcol);  // Wipe it clean.
0258   for (int i = 0; i != theClusterParam.theCluster->size(); ++i) {
0259     auto pix = theClusterParam.theCluster->pixel(i);
0260     int irow = int(pix.x) - row_offset;
0261     int icol = int(pix.y) - col_offset;
0262     // &&& Do we ever get pixels that are out of bounds ???  Need to check.
0263     if ((irow < mrow) & (icol < mcol))
0264       clustMatrix[irow][icol] = float(pix.adc);
0265   }
0266   // &&& Save for later: fillClustMatrix( float * clustMatrix );
0267
0268   //--- Save a copy of clustMatrix into clustMatrix2
0269   memcpy(clustMatrix2, clustMatrix, sizeof(float) * mrow * mcol);
0270
0271   //--- Set both return statuses, since we may call only one.
0272   theClusterParam.ierr = 0;
0273   theClusterParam.ierr2 = 0;
0274
0275   //--- Should we run the 2D reco?
0276   checkRecommend2D(theDetParam, theClusterParam, clusterPayload, ID1);
0277   if (theClusterParam.recommended2D_) {
0278     //--- Call the Template Reco 2d with cluster repair
0279     filled_from_2d = true;
0280     callTempReco2D(theDetParam, theClusterParam, clusterPayload2d, ID2, lp);
0281   } else {
0282     //--- Call the vanilla Template Reco
0283     callTempReco1D(theDetParam, theClusterParam, clusterPayload, ID1, lp);
0284     filled_from_2d = false;
0285   }
0286
0287   //--- Make sure cluster repair returns all info about the hit back up to caller
0288   //--- Necessary because it copied the base class so it does not modify it
0289   theClusterParamBase.isOnEdge_ = theClusterParam.isOnEdge_;
0290   theClusterParamBase.hasBadPixels_ = theClusterParam.hasBadPixels_;
0291   theClusterParamBase.spansTwoROCs_ = theClusterParam.spansTwoROCs_;
0292   theClusterParamBase.hasFilledProb_ = theClusterParam.hasFilledProb_;
0293   theClusterParamBase.qBin_ = theClusterParam.qBin_;
0294   theClusterParamBase.probabilityQ_ = theClusterParam.probabilityQ_;
0295   theClusterParamBase.filled_from_2d = filled_from_2d;
0296   if (filled_from_2d) {
0297     theClusterParamBase.probabilityX_ = theClusterParam.templProbXY_;
0298     theClusterParamBase.probabilityY_ = 0.;
0299   } else {
0300     theClusterParamBase.probabilityX_ = theClusterParam.probabilityX_;
0301     theClusterParamBase.probabilityY_ = theClusterParam.probabilityY_;
0302   }
0303
0304   return LocalPoint(theClusterParam.templXrec_, theClusterParam.templYrec_);
0305 }
0306
0307 //------------------------------------------------------------------
0308 //  Helper function to aggregate call & handling of Template Reco
0309 //------------------------------------------------------------------
0310 void PixelCPEClusterRepair::callTempReco1D(DetParam const& theDetParam,
0311                                            ClusterParamTemplate& theClusterParam,
0312                                            SiPixelTemplateReco::ClusMatrix& clusterPayload,
0313                                            int ID,
0314                                            LocalPoint& lp) const {
0315   SiPixelTemplate templ(thePixelTemp_);
0316
0317   // Output:
0318   float nonsense = -99999.9f;  // nonsense init value
0319   theClusterParam.templXrec_ = theClusterParam.templYrec_ = theClusterParam.templSigmaX_ =
0320       theClusterParam.templSigmaY_ = nonsense;
0321   // If the template recontruction fails, we want to return 1.0 for now
0322   theClusterParam.probabilityX_ = theClusterParam.probabilityY_ = theClusterParam.probabilityQ_ = 1.f;
0323   theClusterParam.qBin_ = 0;
0324   // We have a boolean denoting whether the reco failed or not
0325   theClusterParam.hasFilledProb_ = false;
0326
0327   // In case of template reco failure, these are the lorentz drift corrections
0328   // to be applied
0329   float lorentzshiftX = 0.5f * theDetParam.lorentzShiftInCmX;
0330   float lorentzshiftY = 0.5f * theDetParam.lorentzShiftInCmY;
0331
0332   // ******************************************************************
0333   //--- Call normal TemplateReco
0334   //
0335   float locBz = theDetParam.bz;
0336   float locBx = theDetParam.bx;
0337   //
0338   const bool deadpix = false;
0339   std::vector<std::pair<int, int>> zeropix;
0340   int nypix = 0, nxpix = 0;
0341   //
0342   theClusterParam.ierr = PixelTempReco1D(ID,
0343                                          theClusterParam.cotalpha,
0344                                          theClusterParam.cotbeta,
0345                                          locBz,
0346                                          locBx,
0347                                          clusterPayload,
0348                                          templ,
0349                                          theClusterParam.templYrec_,
0350                                          theClusterParam.templSigmaY_,
0351                                          theClusterParam.probabilityY_,
0352                                          theClusterParam.templXrec_,
0353                                          theClusterParam.templSigmaX_,
0354                                          theClusterParam.probabilityX_,
0355                                          theClusterParam.qBin_,
0356                                          speed_,
0357                                          deadpix,
0358                                          zeropix,
0359                                          theClusterParam.probabilityQ_,
0360                                          nypix,
0361                                          nxpix);
0362   // ******************************************************************
0363
0364   //--- Check exit status
0365   if UNLIKELY (theClusterParam.ierr != 0) {
0366     LogDebug("PixelCPEClusterRepair::localPosition")
0367         << "reconstruction failed with error " << theClusterParam.ierr << "\n";
0368
0369     theClusterParam.probabilityX_ = theClusterParam.probabilityY_ = theClusterParam.probabilityQ_ = 0.f;
0370     theClusterParam.qBin_ = 0;
0371     //
0372     // Template reco has failed, compute position estimates based on cluster center of gravity + Lorentz drift
0373     // Future improvement would be to call generic reco instead
0374
0375     if (theClusterParam.with_track_angle) {
0376       theClusterParam.templXrec_ =
0377           theDetParam.theTopol->localX(theClusterParam.theCluster->x(), theClusterParam.loc_trk_pred) + lorentzshiftX;
0378       theClusterParam.templYrec_ =
0379           theDetParam.theTopol->localY(theClusterParam.theCluster->y(), theClusterParam.loc_trk_pred) + lorentzshiftY;
0380     } else {
0381       edm::LogError("PixelCPEClusterRepair") << "@SUB = PixelCPEClusterRepair::localPosition"
0382                                              << "Should never be here. PixelCPEClusterRepair should always be called "
0383                                                 "with track angles. This is a bad error !!! ";
0384
0385       theClusterParam.templXrec_ = theDetParam.theTopol->localX(theClusterParam.theCluster->x()) + lorentzshiftX;
0386       theClusterParam.templYrec_ = theDetParam.theTopol->localY(theClusterParam.theCluster->y()) + lorentzshiftY;
0387     }
0388   } else {
0389     //--- Template Reco succeeded.  The probabilities are filled.
0390     theClusterParam.hasFilledProb_ = true;
0391
0392     //--- Go from microns to centimeters
0393     theClusterParam.templXrec_ *= micronsToCm;
0394     theClusterParam.templYrec_ *= micronsToCm;
0395
0396     //--- Go back to the module coordinate system
0397     theClusterParam.templXrec_ += lp.x();
0398     theClusterParam.templYrec_ += lp.y();
0399   }
0400   return;
0401 }
0402
0403 //------------------------------------------------------------------
0404 //  Helper function to aggregate call & handling of Template 2D fit
0405 //------------------------------------------------------------------
0406 void PixelCPEClusterRepair::callTempReco2D(DetParam const& theDetParam,
0407                                            ClusterParamTemplate& theClusterParam,
0408                                            SiPixelTemplateReco2D::ClusMatrix& clusterPayload,
0409                                            int ID,
0410                                            LocalPoint& lp) const {
0411   SiPixelTemplate2D templ2d(thePixelTemp2D_);
0412
0413   // Output:
0414   float nonsense = -99999.9f;  // nonsense init value
0415   theClusterParam.templXrec_ = theClusterParam.templYrec_ = theClusterParam.templSigmaX_ =
0416       theClusterParam.templSigmaY_ = nonsense;
0417   // If the template recontruction fails, we want to return 1.0 for now
0418   theClusterParam.probabilityX_ = theClusterParam.probabilityY_ = theClusterParam.probabilityQ_ = 1.f;
0419   theClusterParam.qBin_ = 0;
0420   // We have a boolean denoting whether the reco failed or not
0421   theClusterParam.hasFilledProb_ = false;
0422
0423   // In case of template reco failure, these are the lorentz drift corrections
0424   // to be applied
0425   float lorentzshiftX = 0.5f * theDetParam.lorentzShiftInCmX;
0426   float lorentzshiftY = 0.5f * theDetParam.lorentzShiftInCmY;
0427
0428   // ******************************************************************
0429   //--- Call 2D TemplateReco
0430   //
0431   float locBz = theDetParam.bz;
0432   float locBx = theDetParam.bx;
0433
0434   //--- Input:
0435   //   edgeflagy - (input) flag to indicate the present of edges in y:
0436   //           0=none (or interior gap),1=edge at small y, 2=edge at large y, 3=edge at either end
0437   //
0438   //   edgeflagx - (input) flag to indicate the present of edges in x:
0439   //           0=none, 1=edge at small x, 2=edge at large x
0440   //
0441   //   These two variables are calculated in setTheClu() and stored in edgeTypeX_ and edgeTypeY_
0442   //
0443   //--- Output:
0444   //   deltay - (output) template y-length - cluster length [when > 0, possibly missing end]
0445   //   npixels - ???     &&& Ask Morris
0446
0447   float deltay = 0;  // return param
0448   int npixels = 0;   // return param
0449
0450   //For now, turn off edgeX_ flags
0451   theClusterParam.edgeTypeX_ = 0;
0452
0453   if (clusterPayload.mrow > 4) {
0454     // The cluster is too big, the 2D reco will perform horribly.
0455     // Better to return immediately in error
0456     theClusterParam.ierr2 = 8;
0457
0458   } else {
0459     theClusterParam.ierr2 = PixelTempReco2D(ID,
0460                                             theClusterParam.cotalpha,
0461                                             theClusterParam.cotbeta,
0462                                             locBz,
0463                                             locBx,
0464                                             theClusterParam.edgeTypeY_,
0465                                             theClusterParam.edgeTypeX_,
0466                                             clusterPayload,
0467                                             templ2d,
0468                                             theClusterParam.templYrec_,
0469                                             theClusterParam.templSigmaY_,
0470                                             theClusterParam.templXrec_,
0471                                             theClusterParam.templSigmaX_,
0472                                             theClusterParam.templProbXY_,
0473                                             theClusterParam.probabilityQ_,
0474                                             theClusterParam.qBin_,
0475                                             deltay,
0476                                             npixels);
0477   }
0478   // ******************************************************************
0479
0480   //--- Check exit status
0481   if UNLIKELY (theClusterParam.ierr2 != 0) {
0482     LogDebug("PixelCPEClusterRepair::localPosition")
0483         << "2D reconstruction failed with error " << theClusterParam.ierr2 << "\n";
0484
0485     theClusterParam.probabilityX_ = theClusterParam.probabilityY_ = theClusterParam.probabilityQ_ = 0.f;
0486     theClusterParam.qBin_ = 0;
0487
0488     // 2D Template reco has failed, compute position estimates based on cluster center of gravity + Lorentz drift
0489     // Future improvement would be to call generic reco instead
0490
0491     if (theClusterParam.with_track_angle) {
0492       theClusterParam.templXrec_ =
0493           theDetParam.theTopol->localX(theClusterParam.theCluster->x(), theClusterParam.loc_trk_pred) + lorentzshiftX;
0494       theClusterParam.templYrec_ =
0495           theDetParam.theTopol->localY(theClusterParam.theCluster->y(), theClusterParam.loc_trk_pred) + lorentzshiftY;
0496     } else {
0497       edm::LogError("PixelCPEClusterRepair") << "@SUB = PixelCPEClusterRepair::localPosition"
0498                                              << "Should never be here. PixelCPEClusterRepair should always be called "
0499                                                 "with track angles. This is a bad error !!! ";
0500
0501       theClusterParam.templXrec_ = theDetParam.theTopol->localX(theClusterParam.theCluster->x()) + lorentzshiftX;
0502       theClusterParam.templYrec_ = theDetParam.theTopol->localY(theClusterParam.theCluster->y()) + lorentzshiftY;
0503     }
0504
0505   } else {
0506     //--- Template Reco succeeded.
0507     theClusterParam.hasFilledProb_ = true;
0508
0509     //--- Go from microns to centimeters
0510     theClusterParam.templXrec_ *= micronsToCm;
0511     theClusterParam.templYrec_ *= micronsToCm;
0512
0513     //--- Go back to the module coordinate system
0514     theClusterParam.templXrec_ += lp.x();
0515     theClusterParam.templYrec_ += lp.y();
0516   }
0517   return;
0518 }
0519 //---------------------------------------------------------------------------------
0520 // Helper function to see if we should recommend 2D reco to be run on this
0521 // cluster
0522 // ---------------------------------------------------------------------------------
0523 void PixelCPEClusterRepair::checkRecommend2D(DetParam const& theDetParam,
0524                                              ClusterParamTemplate& theClusterParam,
0525                                              SiPixelTemplateReco::ClusMatrix& clusterPayload,
0526                                              int ID) const
0527
0528 {
0529   // recommend2D is false by default
0530   theClusterParam.recommended2D_ = false;
0531
0532   DetId id = (theDetParam.theDet->geographicalId());
0533
0534   bool recommend = false;
0535   for (auto& rec : recommend2D_) {
0536     recommend = rec.recommend(id, ttopo_);
0537     if (recommend)
0538       break;
0539   }
0540
0541   // only run on those layers recommended by configuration
0542   if (!recommend) {
0543     theClusterParam.recommended2D_ = false;
0544     return;
0545   }
0546
0547   if (theClusterParam.edgeTypeY_) {
0548     // For clusters that have edges in Y, run 2d reco.
0549     // Flags already set by CPEBase
0550     theClusterParam.recommended2D_ = true;
0551     return;
0552   }
0553   // The 1d pixel template
0554   SiPixelTemplate templ(thePixelTemp_);
0555   if (!templ.interpolate(ID, theClusterParam.cotalpha, theClusterParam.cotbeta, theDetParam.bz, theDetParam.bx)) {
0556     //error setting up template, return false
0557     theClusterParam.recommended2D_ = false;
0558     return;
0559   }
0560
0561   //length of the cluster taking into account double sized pixels
0562   float nypix = clusterPayload.mcol;
0563   for (int i = 0; i < clusterPayload.mcol; i++) {
0564     if (clusterPayload.ydouble[i])
0565       nypix += 1.;
0566   }
0567
0568   // templ.clsleny() is the expected length of the cluster along y axis.
0569   // templ.qavg() is the expected total charge of the cluster
0570   // theClusterParam.theCluster->charge() is the total charge of this cluster
0571   float nydiff = templ.clsleny() - nypix;
0572   float qratio = theClusterParam.theCluster->charge() / templ.qavg();
0573
0574   if (nydiff > maxSizeMismatchInY_ && qratio < minChargeRatio_) {
0575     // If the cluster is shorter than expected and has less charge, likely
0576     // due to truncated cluster, try 2D reco
0577
0578     theClusterParam.recommended2D_ = true;
0579     theClusterParam.hasBadPixels_ = true;
0580
0581     // if not RunDamagedClusters flag, don't try to fix any clusters
0582     if (!runDamagedClusters_) {
0583       theClusterParam.recommended2D_ = false;
0584     }
0585
0586     // Figure out what edge flags to set for truncated cluster
0587     // Truncated clusters usually come from dead double columns
0588     //
0589     // If cluster is of even length,  either both of or neither of beginning and ending
0590     // edge are on a double column, so we cannot figure out the likely edge of
0591     // truncation, let the 2D algorithm try extending on both sides (option 3)
0592     if (theClusterParam.theCluster->sizeY() % 2 == 0)
0593       theClusterParam.edgeTypeY_ = 3;
0594     else {
0595       //If the cluster is of odd length, only one of the edges can end on
0596       //a double column, this is the likely edge of truncation
0597       //Double columns always start on even indexes
0598       int min_col = theClusterParam.theCluster->minPixelCol();
0599       if (min_col % 2 == 0) {
0600         //begining edge is at a double column (end edge cannot be,
0601         //because odd length) so likely truncated at small y (option 1)
0602         theClusterParam.edgeTypeY_ = 1;
0603       } else {
0604         //end edge is at a double column (beginning edge cannot be,
0605         //because odd length) so likely truncated at large y (option 2)
0606         theClusterParam.edgeTypeY_ = 2;
0607       }
0608     }
0609   }
0610 }
0611
0612 //------------------------------------------------------------------
0613 //  localError() relies on localPosition() being called FIRST!!!
0614 //------------------------------------------------------------------
0615 LocalError PixelCPEClusterRepair::localError(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0616   ClusterParamTemplate& theClusterParam = static_cast<ClusterParamTemplate&>(theClusterParamBase);
0617
0618   //--- Default is the maximum error used for edge clusters.
0619   //--- (never used, in fact: let comment it out, shut up the complains of the static analyzer, and save a few CPU cycles)
0620   float xerr = 0.0f, yerr = 0.0f;
0621
0622   // Check if the errors were already set at the clusters splitting level
0623   if (theClusterParam.theCluster->getSplitClusterErrorX() > 0.0f &&
0624       theClusterParam.theCluster->getSplitClusterErrorX() < clusterSplitMaxError_ &&
0625       theClusterParam.theCluster->getSplitClusterErrorY() > 0.0f &&
0626       theClusterParam.theCluster->getSplitClusterErrorY() < clusterSplitMaxError_) {
0627     xerr = theClusterParam.theCluster->getSplitClusterErrorX() * micronsToCm;
0628     yerr = theClusterParam.theCluster->getSplitClusterErrorY() * micronsToCm;
0629
0630     //cout << "Errors set at cluster splitting level : " << endl;
0631     //cout << "xerr = " << xerr << endl;
0632     //cout << "yerr = " << yerr << endl;
0633   } else {
0634     // If errors are not split at the cluster splitting level, set the errors here
0635
0636     //--- Check status of both template calls.
0637     if UNLIKELY ((theClusterParam.ierr != 0) || (theClusterParam.ierr2 != 0)) {
0638       // If reconstruction fails the hit position is calculated from cluster center of gravity
0639       // corrected in x by average Lorentz drift. Assign huge errors.
0640       //
0641       if UNLIKELY (!GeomDetEnumerators::isTrackerPixel(theDetParam.thePart))
0642         throw cms::Exception("PixelCPEClusterRepair::localPosition :") << "A non-pixel detector type in here?";
0643
0644       // Assign better errors based on the residuals for failed template cases
0645       if (GeomDetEnumerators::isBarrel(theDetParam.thePart)) {
0646         xerr = 55.0f * micronsToCm;  // &&& get errors from elsewhere?
0647         yerr = 36.0f * micronsToCm;
0648       } else {
0649         xerr = 42.0f * micronsToCm;
0650         yerr = 39.0f * micronsToCm;
0651       }
0652     }
0653     // Use special edge hit errors (derived from observed RMS's) for edge hits that we did not run the 2D reco on
0654     //
0655     else if (!theClusterParamBase.filled_from_2d && (theClusterParam.edgeTypeX_ || theClusterParam.edgeTypeY_)) {
0656       // for edge pixels assign errors according to observed residual RMS
0657       if (theClusterParam.edgeTypeX_ && !theClusterParam.edgeTypeY_) {
0658         xerr = xEdgeXError_ * micronsToCm;
0659         yerr = xEdgeYError_ * micronsToCm;
0660       } else if (!theClusterParam.edgeTypeX_ && theClusterParam.edgeTypeY_) {
0661         xerr = yEdgeXError_ * micronsToCm;
0662         yerr = yEdgeYError_ * micronsToCm;
0663       } else if (theClusterParam.edgeTypeX_ && theClusterParam.edgeTypeY_) {
0664         xerr = bothEdgeXError_ * micronsToCm;
0665         yerr = bothEdgeYError_ * micronsToCm;
0666       }
0667     } else {
0668       xerr = theClusterParam.templSigmaX_ * micronsToCm;
0669       yerr = theClusterParam.templSigmaY_ * micronsToCm;
0670       // &&& should also check ierr (saved as class variable) and return
0671       // &&& nonsense (another class static) if the template fit failed.
0672     }
0673   }
0674
0675   if (theVerboseLevel > 9) {
0676     LogDebug("PixelCPEClusterRepair") << " Sizex = " << theClusterParam.theCluster->sizeX()
0677                                       << " Sizey = " << theClusterParam.theCluster->sizeY()
0678                                       << " Edgex = " << theClusterParam.edgeTypeX_
0679                                       << " Edgey = " << theClusterParam.edgeTypeY_ << " ErrX  = " << xerr
0680                                       << " ErrY  = " << yerr;
0681   }
0682
0683   if (!(xerr > 0.0f))
0684     throw cms::Exception("PixelCPEClusterRepair::localError")
0685         << "\nERROR: Negative pixel error xerr = " << xerr << "\n\n";
0686
0687   if (!(yerr > 0.0f))
0688     throw cms::Exception("PixelCPEClusterRepair::localError")
0689         << "\nERROR: Negative pixel error yerr = " << yerr << "\n\n";
0690
0691   return LocalError(xerr * xerr, 0, yerr * yerr);
0692 }
0693
0694 PixelCPEClusterRepair::Rule::Rule(const std::string& str) {
0695   static const boost::regex rule("([A-Z]+)(\\s+(\\d+))?");
0696   boost::cmatch match;
0697   // match and check it works
0698   if (!regex_match(str.c_str(), match, rule))
0699     throw cms::Exception("Configuration") << "Rule '" << str << "' not understood.\n";
0700
0701   // subdet
0702   subdet_ = -1;
0703   if (strncmp(match[1].first, "PXB", 3) == 0)
0704     subdet_ = PixelSubdetector::PixelBarrel;
0705   else if (strncmp(match[1].first, "PXE", 3) == 0)
0706     subdet_ = PixelSubdetector::PixelEndcap;
0707   if (subdet_ == -1) {
0708     throw cms::Exception("PixelCPEClusterRepair::Configuration")
0709         << "Detector '" << match[1].first << "' not understood. Should be PXB, PXE.\n";
0710   }
0711   // layer (if present)
0712   if (match[3].first != match[3].second) {
0713     layer_ = atoi(match[3].first);
0714   } else {
0715     layer_ = 0;
0716   }
0717 }  //end Rule::Rule
0718
0719 void PixelCPEClusterRepair::fillPSetDescription(edm::ParameterSetDescription& desc) {
0720   desc.add<int>("barrelTemplateID", 0);
0721   desc.add<int>("forwardTemplateID", 0);
0722   desc.add<int>("directoryWithTemplates", 0);
0723   desc.add<int>("speed", -2);
0724   desc.add<bool>("UseClusterSplitter", false);
0725   desc.add<double>("MaxSizeMismatchInY", 0.3);
0726   desc.add<double>("MinChargeRatio", 0.8);
0727   desc.add<std::vector<std::string>>("Recommend2D", {"PXB 2", "PXB 3", "PXB 4"});
0728   desc.add<bool>("RunDamagedClusters", false);
0729 }