Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:00:39 

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