Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-25 05:16:13 

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 std::vector<SiPixelTemplateStore>* templateStore,
0044                                              const SiPixelTemplateDBObject* templateDBobject,
0045                                              const SiPixel2DTemplateDBObject* templateDBobject2D)
0046     : PixelCPEBase(conf, mag, geom, ttopo, lorentzAngle, nullptr, templateDBobject, nullptr, 1) {
0047   LogDebug("PixelCPEClusterRepair::(constructor)") << endl;
0048 
0049   //--- Parameter to decide between DB or text file template access
0050   if (LoadTemplatesFromDB_) {
0051     thePixelTemp_ = templateStore;
0052     // Initialize template store to the selected ID [Morris, 6/25/08]
0053     if (!SiPixelTemplate2D::pushfile(*templateDBobject2D, thePixelTemp2D_))
0054       throw cms::Exception("PixelCPEClusterRepair")
0055           << "\nERROR: Templates not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject2D version "
0056           << (*templateDBobject2D).version() << "\n\n";
0057   } else {
0058     LogDebug("PixelCPEClusterRepair") << "Loading templates for barrel and forward from ASCII files." << endl;
0059     //--- (Archaic) Get configurable template IDs.  This code executes only if we loading pixels from ASCII
0060     //    files, and then they are mandatory.
0061     barrelTemplateID_ = conf.getParameter<int>("barrelTemplateID");
0062     forwardTemplateID_ = conf.getParameter<int>("forwardTemplateID");
0063     templateDir_ = conf.getParameter<int>("directoryWithTemplates");
0064     thePixelTemp_ = &thePixelTempCache_;
0065 
0066     if (!SiPixelTemplate::pushfile(barrelTemplateID_, thePixelTempCache_, templateDir_))
0067       throw cms::Exception("PixelCPEClusterRepair")
0068           << "\nERROR: Template ID " << barrelTemplateID_
0069           << " not loaded correctly from text file. Reconstruction will fail.\n\n";
0070 
0071     if (!SiPixelTemplate::pushfile(forwardTemplateID_, thePixelTempCache_, templateDir_))
0072       throw cms::Exception("PixelCPEClusterRepair")
0073           << "\nERROR: Template ID " << forwardTemplateID_
0074           << " not loaded correctly from text file. Reconstruction will fail.\n\n";
0075   }
0076 
0077   speed_ = conf.getParameter<int>("speed");
0078   goodEdgeAlgo_ = conf.getParameter<bool>("GoodEdgeAlgo");
0079   LogDebug("PixelCPEClusterRepair::PixelCPEClusterRepair:") << "Template speed = " << speed_ << "\n";
0080   LogDebug("PixelCPEClusterRepair::PixelCPEClusterRepair:") << "GoodEdgeAlgo = " << goodEdgeAlgo_ << "\n";
0081 
0082   // this returns the magnetic field value in kgauss (1T = 10 kgauss)
0083   int theMagField = mag->nominalValue();
0084 
0085   if (theMagField >= 36 && theMagField < 39) {
0086     LogDebug("PixelCPEClusterRepair::PixelCPEClusterRepair:")
0087         << "Magnetic field value is: " << theMagField << " kgauss. Algorithm is being run \n";
0088 
0089     templateDBobject2D_ = templateDBobject2D;
0090     fill2DTemplIDs();
0091   }
0092 
0093   UseClusterSplitter_ = conf.getParameter<bool>("UseClusterSplitter");
0094 
0095   maxSizeMismatchInY_ = conf.getParameter<double>("MaxSizeMismatchInY");
0096   minChargeRatio_ = conf.getParameter<double>("MinChargeRatio");
0097 
0098   // read sub-detectors and apply rule to recommend 2D
0099   // can be:
0100   //     XYX (XYZ = PXB, PXE)
0101   //     XYZ n (XYZ as above, n = layer, wheel or disk = 1 .. 6 ;)
0102   std::vector<std::string> str_recommend2D = conf.getParameter<std::vector<std::string>>("Recommend2D");
0103   recommend2D_.reserve(str_recommend2D.size());
0104   for (auto& str : str_recommend2D) {
0105     recommend2D_.push_back(str);
0106   }
0107 
0108   // do not recommend 2D if theMagField!=3.8T
0109   if (theMagField < 36 || theMagField > 39) {
0110     recommend2D_.clear();
0111   }
0112 
0113   // run CR on damaged clusters (and not only on edge hits)
0114   runDamagedClusters_ = conf.getParameter<bool>("RunDamagedClusters");
0115 }
0116 
0117 //-----------------------------------------------------------------------------
0118 //  Fill all 2D template IDs
0119 //-----------------------------------------------------------------------------
0120 void PixelCPEClusterRepair::fill2DTemplIDs() {
0121   auto const& dus = geom_.detUnits();
0122   unsigned m_detectors = dus.size();
0123   for (unsigned int i = 1; i < 7; ++i) {
0124     LogDebug("PixelCPEClusterRepair:: LookingForFirstStrip")
0125         << "Subdetector " << i << " GeomDetEnumerator " << GeomDetEnumerators::tkDetEnum[i] << " offset "
0126         << geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) << " is it strip? "
0127         << (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) != dus.size()
0128                 ? dus[geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i])]->type().isOuterTracker()
0129                 : false);
0130     if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) != dus.size() &&
0131         dus[geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i])]->type().isOuterTracker()) {
0132       if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) < m_detectors)
0133         m_detectors = geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]);
0134     }
0135   }
0136   LogDebug("LookingForFirstStrip") << " Chosen offset: " << m_detectors;
0137 
0138   m_DetParams.resize(m_detectors);
0139   LogDebug("PixelCPEClusterRepair::fillDetParams():") << "caching " << m_detectors << " pixel detectors" << endl;
0140   //Loop through detectors and store 2D template ID
0141   bool printed_info = false;
0142   for (unsigned i = 0; i != m_detectors; ++i) {
0143     auto& p = m_DetParams[i];
0144 
0145     p.detTemplateId2D = templateDBobject2D_->getTemplateID(p.theDet->geographicalId());
0146     if (p.detTemplateId != p.detTemplateId2D && !printed_info) {
0147       edm::LogWarning("PixelCPEClusterRepair")
0148           << "different template ID between 1D and 2D " << p.detTemplateId << " " << p.detTemplateId2D << endl;
0149       printed_info = true;
0150     }
0151   }
0152 }
0153 
0154 //-----------------------------------------------------------------------------
0155 //  Clean up.
0156 //-----------------------------------------------------------------------------
0157 PixelCPEClusterRepair::~PixelCPEClusterRepair() {}
0158 
0159 std::unique_ptr<PixelCPEBase::ClusterParam> PixelCPEClusterRepair::createClusterParam(const SiPixelCluster& cl) const {
0160   return std::make_unique<ClusterParamTemplate>(cl);
0161 }
0162 
0163 //------------------------------------------------------------------
0164 //  Public methods mandated by the base class.
0165 //------------------------------------------------------------------
0166 
0167 //------------------------------------------------------------------
0168 //  The main call to the template code.
0169 //------------------------------------------------------------------
0170 LocalPoint PixelCPEClusterRepair::localPosition(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0171   ClusterParamTemplate& theClusterParam = static_cast<ClusterParamTemplate&>(theClusterParamBase);
0172   bool filled_from_2d = false;
0173 
0174   if (!GeomDetEnumerators::isTrackerPixel(theDetParam.thePart))
0175     throw cms::Exception("PixelCPEClusterRepair::localPosition :") << "A non-pixel detector type in here?";
0176 
0177   int ID1 = -9999;
0178   int ID2 = -9999;
0179   if (LoadTemplatesFromDB_) {
0180     ID1 = theDetParam.detTemplateId;
0181     ID2 = theDetParam.detTemplateId2D;
0182   } else {  // from asci file
0183     if (!GeomDetEnumerators::isEndcap(theDetParam.thePart))
0184       ID1 = ID2 = barrelTemplateID_;  // barrel
0185     else
0186       ID1 = ID2 = forwardTemplateID_;  // forward
0187   }
0188 
0189   // &&& PM, note for later: PixelCPEBase calculates minInX,Y, and maxInX,Y
0190   //     Why can't we simply use that and save time with row_offset, col_offset
0191   //     and mrow = maxInX-minInX, mcol = maxInY-minInY ... Except that we
0192   //     also need to take into account cluster_matrix_size_x,y.
0193 
0194   //--- Preparing to retrieve ADC counts from the SiPixeltheClusterParam.theCluster
0195   //    The pixels from minPixelRow() will go into clust_array_2d[0][*],
0196   //    and the pixels from minPixelCol() will go into clust_array_2d[*][0].
0197   int row_offset = theClusterParam.theCluster->minPixelRow();
0198   int col_offset = theClusterParam.theCluster->minPixelCol();
0199 
0200   //--- Store the coordinates of the center of the (0,0) pixel of the array that
0201   //    gets passed to PixelTempReco2D.  Will add these values to the output of TemplReco2D
0202   float tmp_x = float(row_offset) + 0.5f;
0203   float tmp_y = float(col_offset) + 0.5f;
0204 
0205   //--- Store these offsets (to be added later) in a LocalPoint after tranforming
0206   //    them from measurement units (pixel units) to local coordinates (cm)
0207   LocalPoint lp;
0208   if (theClusterParam.with_track_angle)
0209     //--- Update call with trk angles needed for bow/kink corrections  // Gavril
0210     lp = theDetParam.theTopol->localPosition(MeasurementPoint(tmp_x, tmp_y), theClusterParam.loc_trk_pred);
0211   else {
0212     edm::LogError("PixelCPEClusterRepair") << "@SUB = PixelCPEClusterRepair::localPosition"
0213                                            << "Should never be here. PixelCPEClusterRepair should always be called "
0214                                               "with track angles. This is a bad error !!! ";
0215     lp = theDetParam.theTopol->localPosition(MeasurementPoint(tmp_x, tmp_y));
0216   }
0217 
0218   //--- Compute the size of the matrix which will be passed to TemplateReco.
0219   //    We'll later make  clustMatrix[ mrow ][ mcol ]
0220   int mrow = 0, mcol = 0;
0221   for (int i = 0; i != theClusterParam.theCluster->size(); ++i) {
0222     auto pix = theClusterParam.theCluster->pixel(i);
0223     int irow = int(pix.x);
0224     int icol = int(pix.y);
0225     mrow = std::max(mrow, irow);
0226     mcol = std::max(mcol, icol);
0227   }
0228   mrow -= row_offset;
0229   mrow += 1;
0230   mrow = std::min(mrow, cluster_matrix_size_x);
0231   mcol -= col_offset;
0232   mcol += 1;
0233   mcol = std::min(mcol, cluster_matrix_size_y);
0234   assert(mrow > 0);
0235   assert(mcol > 0);
0236 
0237   //--- Make and fill the bool arrays flagging double pixels
0238   bool xdouble[mrow], ydouble[mcol];
0239   // x directions (shorter), rows
0240   for (int irow = 0; irow < mrow; ++irow)
0241     xdouble[irow] = theDetParam.theTopol->isItBigPixelInX(irow + row_offset);
0242   //
0243   // y directions (longer), columns
0244   for (int icol = 0; icol < mcol; ++icol)
0245     ydouble[icol] = theDetParam.theTopol->isItBigPixelInY(icol + col_offset);
0246 
0247   //--- C-style matrix.  We'll need it in either case.
0248   float clustMatrix[mrow][mcol];
0249   float clustMatrix2[mrow][mcol];
0250 
0251   //--- Prepare struct that passes pointers to TemplateReco.  It doesn't own anything.
0252   SiPixelTemplateReco::ClusMatrix clusterPayload{&clustMatrix[0][0], xdouble, ydouble, mrow, mcol};
0253   SiPixelTemplateReco2D::ClusMatrix clusterPayload2d{&clustMatrix2[0][0], xdouble, ydouble, mrow, mcol};
0254 
0255   //--- Copy clust's pixels (calibrated in electrons) into clustMatrix;
0256   memset(clustMatrix, 0, sizeof(float) * mrow * mcol);  // Wipe it clean.
0257   for (int i = 0; i != theClusterParam.theCluster->size(); ++i) {
0258     auto pix = theClusterParam.theCluster->pixel(i);
0259     int irow = int(pix.x) - row_offset;
0260     int icol = int(pix.y) - col_offset;
0261     // &&& Do we ever get pixels that are out of bounds ???  Need to check.
0262     if ((irow < mrow) & (icol < mcol))
0263       clustMatrix[irow][icol] = float(pix.adc);
0264   }
0265   // &&& Save for later: fillClustMatrix( float * clustMatrix );
0266 
0267   //--- Save a copy of clustMatrix into clustMatrix2
0268   memcpy(clustMatrix2, clustMatrix, sizeof(float) * mrow * mcol);
0269 
0270   //--- Set both return statuses, since we may call only one.
0271   theClusterParam.ierr = 0;
0272   theClusterParam.ierr2 = 0;
0273 
0274   //--- Should we run the 2D reco?
0275   checkRecommend2D(theDetParam, theClusterParam, clusterPayload, ID1);
0276   if (theClusterParam.recommended2D_) {
0277     //--- Call the Template Reco 2d with cluster repair
0278     filled_from_2d = true;
0279     callTempReco2D(theDetParam, theClusterParam, clusterPayload2d, ID2, lp);
0280   } else {
0281     //--- Call the vanilla Template Reco
0282     callTempReco1D(theDetParam, theClusterParam, clusterPayload, ID1, lp);
0283     filled_from_2d = false;
0284   }
0285 
0286   //--- Make sure cluster repair returns all info about the hit back up to caller
0287   //--- Necessary because it copied the base class so it does not modify it
0288   theClusterParamBase.isOnEdge_ = theClusterParam.isOnEdge_;
0289   theClusterParamBase.hasBadPixels_ = theClusterParam.hasBadPixels_;
0290   theClusterParamBase.spansTwoROCs_ = theClusterParam.spansTwoROCs_;
0291   theClusterParamBase.hasFilledProb_ = theClusterParam.hasFilledProb_;
0292   theClusterParamBase.qBin_ = theClusterParam.qBin_;
0293   theClusterParamBase.probabilityQ_ = theClusterParam.probabilityQ_;
0294   theClusterParamBase.filled_from_2d = filled_from_2d;
0295   if (filled_from_2d) {
0296     theClusterParamBase.probabilityX_ = theClusterParam.templProbXY_;
0297     theClusterParamBase.probabilityY_ = 0.;
0298   } else {
0299     theClusterParamBase.probabilityX_ = theClusterParam.probabilityX_;
0300     theClusterParamBase.probabilityY_ = theClusterParam.probabilityY_;
0301   }
0302 
0303   return LocalPoint(theClusterParam.templXrec_, theClusterParam.templYrec_);
0304 }
0305 
0306 //------------------------------------------------------------------
0307 //  Helper function to aggregate call & handling of Template Reco
0308 //------------------------------------------------------------------
0309 void PixelCPEClusterRepair::callTempReco1D(DetParam const& theDetParam,
0310                                            ClusterParamTemplate& theClusterParam,
0311                                            SiPixelTemplateReco::ClusMatrix& clusterPayload,
0312                                            int ID,
0313                                            LocalPoint& lp) const {
0314   SiPixelTemplate templ(*thePixelTemp_);
0315 
0316   // Output:
0317   float nonsense = -99999.9f;  // nonsense init value
0318   theClusterParam.templXrec_ = theClusterParam.templYrec_ = theClusterParam.templSigmaX_ =
0319       theClusterParam.templSigmaY_ = nonsense;
0320   // If the template recontruction fails, we want to return 1.0 for now
0321   theClusterParam.probabilityX_ = theClusterParam.probabilityY_ = theClusterParam.probabilityQ_ = 1.f;
0322   theClusterParam.qBin_ = 0;
0323   // We have a boolean denoting whether the reco failed or not
0324   theClusterParam.hasFilledProb_ = false;
0325 
0326   // In case of template reco failure, these are the lorentz drift corrections
0327   // to be applied
0328   float lorentzshiftX = 0.5f * theDetParam.lorentzShiftInCmX;
0329   float lorentzshiftY = 0.5f * theDetParam.lorentzShiftInCmY;
0330 
0331   // ******************************************************************
0332   //--- Call normal TemplateReco
0333   //
0334   float locBz = theDetParam.bz;
0335   float locBx = theDetParam.bx;
0336   //
0337   const bool deadpix = false;
0338   std::vector<std::pair<int, int>> zeropix;
0339   int nypix = 0, nxpix = 0;
0340   //
0341   theClusterParam.ierr = PixelTempReco1D(ID,
0342                                          theClusterParam.cotalpha,
0343                                          theClusterParam.cotbeta,
0344                                          locBz,
0345                                          locBx,
0346                                          clusterPayload,
0347                                          templ,
0348                                          theClusterParam.templYrec_,
0349                                          theClusterParam.templSigmaY_,
0350                                          theClusterParam.probabilityY_,
0351                                          theClusterParam.templXrec_,
0352                                          theClusterParam.templSigmaX_,
0353                                          theClusterParam.probabilityX_,
0354                                          theClusterParam.qBin_,
0355                                          speed_,
0356                                          deadpix,
0357                                          zeropix,
0358                                          theClusterParam.probabilityQ_,
0359                                          nypix,
0360                                          nxpix,
0361                                          goodEdgeAlgo_);
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(
0556           ID, theClusterParam.cotalpha, theClusterParam.cotbeta, theDetParam.bz, theDetParam.bx, goodEdgeAlgo_)) {
0557     //error setting up template, return false
0558     theClusterParam.recommended2D_ = false;
0559     return;
0560   }
0561 
0562   //length of the cluster taking into account double sized pixels
0563   float nypix = clusterPayload.mcol;
0564   for (int i = 0; i < clusterPayload.mcol; i++) {
0565     if (clusterPayload.ydouble[i])
0566       nypix += 1.;
0567   }
0568 
0569   // templ.clsleny() is the expected length of the cluster along y axis.
0570   // templ.qavg() is the expected total charge of the cluster
0571   // theClusterParam.theCluster->charge() is the total charge of this cluster
0572   float nydiff = templ.clsleny() - nypix;
0573   float qratio = theClusterParam.theCluster->charge() / templ.qavg();
0574 
0575   if (nydiff > maxSizeMismatchInY_ && qratio < minChargeRatio_) {
0576     // If the cluster is shorter than expected and has less charge, likely
0577     // due to truncated cluster, try 2D reco
0578 
0579     theClusterParam.recommended2D_ = true;
0580     theClusterParam.hasBadPixels_ = true;
0581 
0582     // if not RunDamagedClusters flag, don't try to fix any clusters
0583     if (!runDamagedClusters_) {
0584       theClusterParam.recommended2D_ = false;
0585     }
0586 
0587     // Figure out what edge flags to set for truncated cluster
0588     // Truncated clusters usually come from dead double columns
0589     //
0590     // If cluster is of even length,  either both of or neither of beginning and ending
0591     // edge are on a double column, so we cannot figure out the likely edge of
0592     // truncation, let the 2D algorithm try extending on both sides (option 3)
0593     if (theClusterParam.theCluster->sizeY() % 2 == 0)
0594       theClusterParam.edgeTypeY_ = 3;
0595     else {
0596       //If the cluster is of odd length, only one of the edges can end on
0597       //a double column, this is the likely edge of truncation
0598       //Double columns always start on even indexes
0599       int min_col = theClusterParam.theCluster->minPixelCol();
0600       if (min_col % 2 == 0) {
0601         //begining edge is at a double column (end edge cannot be,
0602         //because odd length) so likely truncated at small y (option 1)
0603         theClusterParam.edgeTypeY_ = 1;
0604       } else {
0605         //end edge is at a double column (beginning edge cannot be,
0606         //because odd length) so likely truncated at large y (option 2)
0607         theClusterParam.edgeTypeY_ = 2;
0608       }
0609     }
0610   }
0611 }
0612 
0613 //------------------------------------------------------------------
0614 //  localError() relies on localPosition() being called FIRST!!!
0615 //------------------------------------------------------------------
0616 LocalError PixelCPEClusterRepair::localError(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0617   ClusterParamTemplate& theClusterParam = static_cast<ClusterParamTemplate&>(theClusterParamBase);
0618 
0619   //--- Default is the maximum error used for edge clusters.
0620   //--- (never used, in fact: let comment it out, shut up the complains of the static analyzer, and save a few CPU cycles)
0621   float xerr = 0.0f, yerr = 0.0f;
0622 
0623   // Check if the errors were already set at the clusters splitting level
0624   if (theClusterParam.theCluster->getSplitClusterErrorX() > 0.0f &&
0625       theClusterParam.theCluster->getSplitClusterErrorX() < clusterSplitMaxError_ &&
0626       theClusterParam.theCluster->getSplitClusterErrorY() > 0.0f &&
0627       theClusterParam.theCluster->getSplitClusterErrorY() < clusterSplitMaxError_) {
0628     xerr = theClusterParam.theCluster->getSplitClusterErrorX() * micronsToCm;
0629     yerr = theClusterParam.theCluster->getSplitClusterErrorY() * micronsToCm;
0630 
0631     //cout << "Errors set at cluster splitting level : " << endl;
0632     //cout << "xerr = " << xerr << endl;
0633     //cout << "yerr = " << yerr << endl;
0634   } else {
0635     // If errors are not split at the cluster splitting level, set the errors here
0636 
0637     //--- Check status of both template calls.
0638     if UNLIKELY ((theClusterParam.ierr != 0) || (theClusterParam.ierr2 != 0)) {
0639       // If reconstruction fails the hit position is calculated from cluster center of gravity
0640       // corrected in x by average Lorentz drift. Assign huge errors.
0641       //
0642       if UNLIKELY (!GeomDetEnumerators::isTrackerPixel(theDetParam.thePart))
0643         throw cms::Exception("PixelCPEClusterRepair::localPosition :") << "A non-pixel detector type in here?";
0644 
0645       // Assign better errors based on the residuals for failed template cases
0646       if (GeomDetEnumerators::isBarrel(theDetParam.thePart)) {
0647         xerr = 55.0f * micronsToCm;  // &&& get errors from elsewhere?
0648         yerr = 36.0f * micronsToCm;
0649       } else {
0650         xerr = 42.0f * micronsToCm;
0651         yerr = 39.0f * micronsToCm;
0652       }
0653     }
0654     // Use special edge hit errors (derived from observed RMS's) for edge hits that we did not run the 2D reco on
0655     //
0656     else if (!theClusterParamBase.filled_from_2d && (theClusterParam.edgeTypeX_ || theClusterParam.edgeTypeY_)) {
0657       // for edge pixels assign errors according to observed residual RMS
0658       if (theClusterParam.edgeTypeX_ && !theClusterParam.edgeTypeY_) {
0659         xerr = xEdgeXError_ * micronsToCm;
0660         yerr = xEdgeYError_ * micronsToCm;
0661       } else if (!theClusterParam.edgeTypeX_ && theClusterParam.edgeTypeY_) {
0662         xerr = yEdgeXError_ * micronsToCm;
0663         yerr = yEdgeYError_ * micronsToCm;
0664       } else if (theClusterParam.edgeTypeX_ && theClusterParam.edgeTypeY_) {
0665         xerr = bothEdgeXError_ * micronsToCm;
0666         yerr = bothEdgeYError_ * micronsToCm;
0667       }
0668     } else {
0669       xerr = theClusterParam.templSigmaX_ * micronsToCm;
0670       yerr = theClusterParam.templSigmaY_ * micronsToCm;
0671       // &&& should also check ierr (saved as class variable) and return
0672       // &&& nonsense (another class static) if the template fit failed.
0673     }
0674   }
0675 
0676   if (theVerboseLevel > 9) {
0677     LogDebug("PixelCPEClusterRepair") << " Sizex = " << theClusterParam.theCluster->sizeX()
0678                                       << " Sizey = " << theClusterParam.theCluster->sizeY()
0679                                       << " Edgex = " << theClusterParam.edgeTypeX_
0680                                       << " Edgey = " << theClusterParam.edgeTypeY_ << " ErrX  = " << xerr
0681                                       << " ErrY  = " << yerr;
0682   }
0683 
0684   if (!(xerr > 0.0f))
0685     throw cms::Exception("PixelCPEClusterRepair::localError")
0686         << "\nERROR: Negative pixel error xerr = " << xerr << "\n\n";
0687 
0688   if (!(yerr > 0.0f))
0689     throw cms::Exception("PixelCPEClusterRepair::localError")
0690         << "\nERROR: Negative pixel error yerr = " << yerr << "\n\n";
0691 
0692   return LocalError(xerr * xerr, 0, yerr * yerr);
0693 }
0694 
0695 PixelCPEClusterRepair::Rule::Rule(const std::string& str) {
0696   static const boost::regex rule("([A-Z]+)(\\s+(\\d+))?");
0697   boost::cmatch match;
0698   // match and check it works
0699   if (!regex_match(str.c_str(), match, rule))
0700     throw cms::Exception("Configuration") << "Rule '" << str << "' not understood.\n";
0701 
0702   // subdet
0703   subdet_ = -1;
0704   if (strncmp(match[1].first, "PXB", 3) == 0)
0705     subdet_ = PixelSubdetector::PixelBarrel;
0706   else if (strncmp(match[1].first, "PXE", 3) == 0)
0707     subdet_ = PixelSubdetector::PixelEndcap;
0708   if (subdet_ == -1) {
0709     throw cms::Exception("PixelCPEClusterRepair::Configuration")
0710         << "Detector '" << match[1].first << "' not understood. Should be PXB, PXE.\n";
0711   }
0712   // layer (if present)
0713   if (match[3].first != match[3].second) {
0714     layer_ = atoi(match[3].first);
0715   } else {
0716     layer_ = 0;
0717   }
0718 }  //end Rule::Rule
0719 
0720 void PixelCPEClusterRepair::fillPSetDescription(edm::ParameterSetDescription& desc) {
0721   desc.add<int>("barrelTemplateID", 0);
0722   desc.add<int>("forwardTemplateID", 0);
0723   desc.add<int>("directoryWithTemplates", 0);
0724   desc.add<int>("speed", -2);
0725   desc.add<bool>("GoodEdgeAlgo", false);
0726   desc.add<bool>("UseClusterSplitter", false);
0727   desc.add<double>("MaxSizeMismatchInY", 0.3);
0728   desc.add<double>("MinChargeRatio", 0.8);
0729   desc.add<std::vector<std::string>>("Recommend2D", {"PXB 2", "PXB 3", "PXB 4"});
0730   desc.add<bool>("RunDamagedClusters", false);
0731 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!