CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PixelCPEClusterRepair.cc
Go to the documentation of this file.
1 // Include our own header first
3 
4 // Geometry services
7 
8 // MessageLogger
10 
11 // Magnetic field
13 
14 // Commented for now (3/10/17) until we figure out how to resuscitate 2D template splitter
16 
17 #include <vector>
18 #include "boost/multi_array.hpp"
19 #include <boost/regex.hpp>
20 #include <map>
21 
22 #include <iostream>
23 
24 using namespace SiPixelTemplateReco;
25 //using namespace SiPixelTemplateSplit;
26 using namespace std;
27 
28 namespace {
29  constexpr float micronsToCm = 1.0e-4;
30  constexpr int cluster_matrix_size_x = 13;
31  constexpr int cluster_matrix_size_y = 21;
32 } // namespace
33 
34 //-----------------------------------------------------------------------------
35 // Constructor.
36 //
37 //-----------------------------------------------------------------------------
39  const MagneticField* mag,
40  const TrackerGeometry& geom,
41  const TrackerTopology& ttopo,
42  const SiPixelLorentzAngle* lorentzAngle,
43  const SiPixelTemplateDBObject* templateDBobject,
44  const SiPixel2DTemplateDBObject* templateDBobject2D)
45  : PixelCPEBase(conf, mag, geom, ttopo, lorentzAngle, nullptr, templateDBobject, nullptr, 1) {
46  LogDebug("PixelCPEClusterRepair::(constructor)") << endl;
47 
48  //--- Parameter to decide between DB or text file template access
50  // Initialize template store to the selected ID [Morris, 6/25/08]
52  throw cms::Exception("PixelCPEClusterRepair")
53  << "\nERROR: Templates not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject version "
54  << (*templateDBobject_).version() << "\n\n";
55 
56  // Initialize template store to the selected ID [Morris, 6/25/08]
57  if (!SiPixelTemplate2D::pushfile(*templateDBobject2D, thePixelTemp2D_))
58  throw cms::Exception("PixelCPEClusterRepair")
59  << "\nERROR: Templates not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject2D version "
60  << (*templateDBobject2D).version() << "\n\n";
61  } else {
62  LogDebug("PixelCPEClusterRepair") << "Loading templates for barrel and forward from ASCII files." << endl;
63  //--- (Archaic) Get configurable template IDs. This code executes only if we loading pixels from ASCII
64  // files, and then they are mandatory.
65  barrelTemplateID_ = conf.getParameter<int>("barrelTemplateID");
66  forwardTemplateID_ = conf.getParameter<int>("forwardTemplateID");
67  templateDir_ = conf.getParameter<int>("directoryWithTemplates");
68 
70  throw cms::Exception("PixelCPEClusterRepair")
71  << "\nERROR: Template ID " << barrelTemplateID_
72  << " not loaded correctly from text file. Reconstruction will fail.\n\n";
73 
74  if (!SiPixelTemplate::pushfile(forwardTemplateID_, thePixelTemp_, templateDir_))
75  throw cms::Exception("PixelCPEClusterRepair")
76  << "\nERROR: Template ID " << forwardTemplateID_
77  << " not loaded correctly from text file. Reconstruction will fail.\n\n";
78  }
79 
80  speed_ = conf.getParameter<int>("speed");
81  LogDebug("PixelCPEClusterRepair::PixelCPEClusterRepair:") << "Template speed = " << speed_ << "\n";
82 
83  // this returns the magnetic field value in kgauss (1T = 10 kgauss)
84  int theMagField = mag->nominalValue();
85 
86  if (theMagField >= 36 && theMagField < 39) {
87  LogDebug("PixelCPEClusterRepair::PixelCPEClusterRepair:")
88  << "Magnetic field value is: " << theMagField << " kgauss. Algorithm is being run \n";
89 
90  templateDBobject2D_ = templateDBobject2D;
92  }
93 
94  UseClusterSplitter_ = conf.getParameter<bool>("UseClusterSplitter");
95 
96  maxSizeMismatchInY_ = conf.getParameter<double>("MaxSizeMismatchInY");
97  minChargeRatio_ = conf.getParameter<double>("MinChargeRatio");
98 
99  // read sub-detectors and apply rule to recommend 2D
100  // can be:
101  // XYX (XYZ = PXB, PXE)
102  // XYZ n (XYZ as above, n = layer, wheel or disk = 1 .. 6 ;)
103  std::vector<std::string> str_recommend2D = conf.getParameter<std::vector<std::string>>("Recommend2D");
104  recommend2D_.reserve(str_recommend2D.size());
105  for (auto& str : str_recommend2D) {
106  recommend2D_.push_back(str);
107  }
108 
109  // do not recommend 2D if theMagField!=3.8T
110  if (theMagField < 36 || theMagField > 39) {
111  recommend2D_.clear();
112  }
113 
114  // run CR on damaged clusters (and not only on edge hits)
115  runDamagedClusters_ = conf.getParameter<bool>("RunDamagedClusters");
116 }
117 
118 //-----------------------------------------------------------------------------
119 // Fill all 2D template IDs
120 //-----------------------------------------------------------------------------
122  auto const& dus = geom_.detUnits();
123  unsigned m_detectors = dus.size();
124  for (unsigned int i = 1; i < 7; ++i) {
125  LogDebug("PixelCPEClusterRepair:: LookingForFirstStrip")
126  << "Subdetector " << i << " GeomDetEnumerator " << GeomDetEnumerators::tkDetEnum[i] << " offset "
127  << geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) << " is it strip? "
128  << (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) != dus.size()
129  ? dus[geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i])]->type().isOuterTracker()
130  : false);
131  if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) != dus.size() &&
132  dus[geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i])]->type().isOuterTracker()) {
133  if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) < m_detectors)
134  m_detectors = geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]);
135  }
136  }
137  LogDebug("LookingForFirstStrip") << " Chosen offset: " << m_detectors;
138 
139  m_DetParams.resize(m_detectors);
140  LogDebug("PixelCPEClusterRepair::fillDetParams():") << "caching " << m_detectors << " pixel detectors" << endl;
141  //Loop through detectors and store 2D template ID
142  bool printed_info = false;
143  for (unsigned i = 0; i != m_detectors; ++i) {
144  auto& p = m_DetParams[i];
145 
146  p.detTemplateId2D = templateDBobject2D_->getTemplateID(p.theDet->geographicalId());
147  if (p.detTemplateId != p.detTemplateId2D && !printed_info) {
148  edm::LogWarning("PixelCPEClusterRepair")
149  << "different template ID between 1D and 2D " << p.detTemplateId << " " << p.detTemplateId2D << endl;
150  printed_info = true;
151  }
152  }
153 }
154 
155 //-----------------------------------------------------------------------------
156 // Clean up.
157 //-----------------------------------------------------------------------------
159  for (auto x : thePixelTemp_)
160  x.destroy();
161  for (auto x : thePixelTemp2D_)
162  x.destroy();
163 }
164 
166  return new ClusterParamTemplate(cl);
167 }
168 
169 //------------------------------------------------------------------
170 // Public methods mandated by the base class.
171 //------------------------------------------------------------------
172 
173 //------------------------------------------------------------------
174 // The main call to the template code.
175 //------------------------------------------------------------------
176 LocalPoint PixelCPEClusterRepair::localPosition(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
177  ClusterParamTemplate& theClusterParam = static_cast<ClusterParamTemplate&>(theClusterParamBase);
178  bool filled_from_2d = false;
179 
181  throw cms::Exception("PixelCPEClusterRepair::localPosition :") << "A non-pixel detector type in here?";
182 
183  int ID1 = -9999;
184  int ID2 = -9999;
185  if (LoadTemplatesFromDB_) {
186  ID1 = theDetParam.detTemplateId;
187  ID2 = theDetParam.detTemplateId2D;
188  } else { // from asci file
189  if (!GeomDetEnumerators::isEndcap(theDetParam.thePart))
190  ID1 = ID2 = barrelTemplateID_; // barrel
191  else
192  ID1 = ID2 = forwardTemplateID_; // forward
193  }
194 
195  // &&& PM, note for later: PixelCPEBase calculates minInX,Y, and maxInX,Y
196  // Why can't we simply use that and save time with row_offset, col_offset
197  // and mrow = maxInX-minInX, mcol = maxInY-minInY ... Except that we
198  // also need to take into account cluster_matrix_size_x,y.
199 
200  //--- Preparing to retrieve ADC counts from the SiPixeltheClusterParam.theCluster
201  // The pixels from minPixelRow() will go into clust_array_2d[0][*],
202  // and the pixels from minPixelCol() will go into clust_array_2d[*][0].
203  int row_offset = theClusterParam.theCluster->minPixelRow();
204  int col_offset = theClusterParam.theCluster->minPixelCol();
205 
206  //--- Store the coordinates of the center of the (0,0) pixel of the array that
207  // gets passed to PixelTempReco2D. Will add these values to the output of TemplReco2D
208  float tmp_x = float(row_offset) + 0.5f;
209  float tmp_y = float(col_offset) + 0.5f;
210 
211  //--- Store these offsets (to be added later) in a LocalPoint after tranforming
212  // them from measurement units (pixel units) to local coordinates (cm)
213  LocalPoint lp;
214  if (theClusterParam.with_track_angle)
215  //--- Update call with trk angles needed for bow/kink corrections // Gavril
216  lp = theDetParam.theTopol->localPosition(MeasurementPoint(tmp_x, tmp_y), theClusterParam.loc_trk_pred);
217  else {
218  edm::LogError("PixelCPEClusterRepair") << "@SUB = PixelCPEClusterRepair::localPosition"
219  << "Should never be here. PixelCPEClusterRepair should always be called "
220  "with track angles. This is a bad error !!! ";
221  lp = theDetParam.theTopol->localPosition(MeasurementPoint(tmp_x, tmp_y));
222  }
223 
224  //--- Compute the size of the matrix which will be passed to TemplateReco.
225  // We'll later make clustMatrix[ mrow ][ mcol ]
226  int mrow = 0, mcol = 0;
227  for (int i = 0; i != theClusterParam.theCluster->size(); ++i) {
228  auto pix = theClusterParam.theCluster->pixel(i);
229  int irow = int(pix.x);
230  int icol = int(pix.y);
231  mrow = std::max(mrow, irow);
232  mcol = std::max(mcol, icol);
233  }
234  mrow -= row_offset;
235  mrow += 1;
236  mrow = std::min(mrow, cluster_matrix_size_x);
237  mcol -= col_offset;
238  mcol += 1;
239  mcol = std::min(mcol, cluster_matrix_size_y);
240  assert(mrow > 0);
241  assert(mcol > 0);
242 
243  //--- Make and fill the bool arrays flagging double pixels
244  bool xdouble[mrow], ydouble[mcol];
245  // x directions (shorter), rows
246  for (int irow = 0; irow < mrow; ++irow)
247  xdouble[irow] = theDetParam.theRecTopol->isItBigPixelInX(irow + row_offset);
248  //
249  // y directions (longer), columns
250  for (int icol = 0; icol < mcol; ++icol)
251  ydouble[icol] = theDetParam.theRecTopol->isItBigPixelInY(icol + col_offset);
252 
253  //--- C-style matrix. We'll need it in either case.
254  float clustMatrix[mrow][mcol];
255  float clustMatrix2[mrow][mcol];
256 
257  //--- Prepare struct that passes pointers to TemplateReco. It doesn't own anything.
258  SiPixelTemplateReco::ClusMatrix clusterPayload{&clustMatrix[0][0], xdouble, ydouble, mrow, mcol};
259  SiPixelTemplateReco2D::ClusMatrix clusterPayload2d{&clustMatrix2[0][0], xdouble, ydouble, mrow, mcol};
260 
261  //--- Copy clust's pixels (calibrated in electrons) into clustMatrix;
262  memset(clustMatrix, 0, sizeof(float) * mrow * mcol); // Wipe it clean.
263  for (int i = 0; i != theClusterParam.theCluster->size(); ++i) {
264  auto pix = theClusterParam.theCluster->pixel(i);
265  int irow = int(pix.x) - row_offset;
266  int icol = int(pix.y) - col_offset;
267  // &&& Do we ever get pixels that are out of bounds ??? Need to check.
268  if ((irow < mrow) & (icol < mcol))
269  clustMatrix[irow][icol] = float(pix.adc);
270  }
271  // &&& Save for later: fillClustMatrix( float * clustMatrix );
272 
273  //--- Save a copy of clustMatrix into clustMatrix2
274  memcpy(clustMatrix2, clustMatrix, sizeof(float) * mrow * mcol);
275 
276  //--- Set both return statuses, since we may call only one.
277  theClusterParam.ierr = 0;
278  theClusterParam.ierr2 = 0;
279 
280  //--- Should we run the 2D reco?
281  checkRecommend2D(theDetParam, theClusterParam, clusterPayload, ID1);
282  if (theClusterParam.recommended2D_) {
283  //--- Call the Template Reco 2d with cluster repair
284  filled_from_2d = true;
285  callTempReco2D(theDetParam, theClusterParam, clusterPayload2d, ID2, lp);
286  } else {
287  //--- Call the vanilla Template Reco
288  callTempReco1D(theDetParam, theClusterParam, clusterPayload, ID1, lp);
289  filled_from_2d = false;
290  }
291 
292  //--- Make sure cluster repair returns all info about the hit back up to caller
293  //--- Necessary because it copied the base class so it does not modify it
294  theClusterParamBase.isOnEdge_ = theClusterParam.isOnEdge_;
295  theClusterParamBase.hasBadPixels_ = theClusterParam.hasBadPixels_;
296  theClusterParamBase.spansTwoROCs_ = theClusterParam.spansTwoROCs_;
297  theClusterParamBase.hasFilledProb_ = theClusterParam.hasFilledProb_;
298  theClusterParamBase.qBin_ = theClusterParam.qBin_;
299  theClusterParamBase.probabilityQ_ = theClusterParam.probabilityQ_;
300  theClusterParamBase.filled_from_2d = filled_from_2d;
301  if (filled_from_2d) {
302  theClusterParamBase.probabilityX_ = theClusterParam.templProbXY_;
303  theClusterParamBase.probabilityY_ = 0.;
304  } else {
305  theClusterParamBase.probabilityX_ = theClusterParam.probabilityX_;
306  theClusterParamBase.probabilityY_ = theClusterParam.probabilityY_;
307  }
308 
309  return LocalPoint(theClusterParam.templXrec_, theClusterParam.templYrec_);
310 }
311 
312 //------------------------------------------------------------------
313 // Helper function to aggregate call & handling of Template Reco
314 //------------------------------------------------------------------
316  ClusterParamTemplate& theClusterParam,
317  SiPixelTemplateReco::ClusMatrix& clusterPayload,
318  int ID,
319  LocalPoint& lp) const {
321 
322  // Output:
323  float nonsense = -99999.9f; // nonsense init value
324  theClusterParam.templXrec_ = theClusterParam.templYrec_ = theClusterParam.templSigmaX_ =
325  theClusterParam.templSigmaY_ = nonsense;
326  // If the template recontruction fails, we want to return 1.0 for now
327  theClusterParam.probabilityX_ = theClusterParam.probabilityY_ = theClusterParam.probabilityQ_ = 1.f;
328  theClusterParam.qBin_ = 0;
329  // We have a boolean denoting whether the reco failed or not
330  theClusterParam.hasFilledProb_ = false;
331 
332  // ******************************************************************
333  //--- Call normal TemplateReco
334  //
335  float locBz = theDetParam.bz;
336  float locBx = theDetParam.bx;
337  //
338  const bool deadpix = false;
339  std::vector<std::pair<int, int>> zeropix;
340  int nypix = 0, nxpix = 0;
341  //
342  theClusterParam.ierr = PixelTempReco1D(ID,
343  theClusterParam.cotalpha,
344  theClusterParam.cotbeta,
345  locBz,
346  locBx,
347  clusterPayload,
348  templ,
349  theClusterParam.templYrec_,
350  theClusterParam.templSigmaY_,
351  theClusterParam.probabilityY_,
352  theClusterParam.templXrec_,
353  theClusterParam.templSigmaX_,
354  theClusterParam.probabilityX_,
355  theClusterParam.qBin_,
356  speed_,
357  deadpix,
358  zeropix,
359  theClusterParam.probabilityQ_,
360  nypix,
361  nxpix);
362  // ******************************************************************
363 
364  //--- Check exit status
365  if
366  UNLIKELY(theClusterParam.ierr != 0) {
367  LogDebug("PixelCPEClusterRepair::localPosition")
368  << "reconstruction failed with error " << theClusterParam.ierr << "\n";
369 
370  theClusterParam.probabilityX_ = theClusterParam.probabilityY_ = theClusterParam.probabilityQ_ = 0.f;
371  theClusterParam.qBin_ = 0;
372 
373  // Gavril: what do we do in this case ? For now, just return the cluster center of gravity in microns
374  // In the x case, apply a rough Lorentz drift average correction
375  // To do: call PixelCPEGeneric whenever PixelTempReco1D fails
376  float lorentz_drift = -999.9;
377  if (!GeomDetEnumerators::isEndcap(theDetParam.thePart))
378  lorentz_drift = 60.0f; // in microns
379  else
380  lorentz_drift = 10.0f; // in microns
381  // GG: trk angles needed to correct for bows/kinks
382  if (theClusterParam.with_track_angle) {
383  theClusterParam.templXrec_ =
384  theDetParam.theTopol->localX(theClusterParam.theCluster->x(), theClusterParam.loc_trk_pred) -
385  lorentz_drift * micronsToCm; // rough Lorentz drift correction
386  theClusterParam.templYrec_ =
387  theDetParam.theTopol->localY(theClusterParam.theCluster->y(), theClusterParam.loc_trk_pred);
388  } else {
389  edm::LogError("PixelCPEClusterRepair") << "@SUB = PixelCPEClusterRepair::localPosition"
390  << "Should never be here. PixelCPEClusterRepair should always be called "
391  "with track angles. This is a bad error !!! ";
392 
393  theClusterParam.templXrec_ = theDetParam.theTopol->localX(theClusterParam.theCluster->x()) -
394  lorentz_drift * micronsToCm; // rough Lorentz drift correction
395  theClusterParam.templYrec_ = theDetParam.theTopol->localY(theClusterParam.theCluster->y());
396  }
397  }
398  else {
399  //--- Template Reco succeeded. The probabilities are filled.
400  theClusterParam.hasFilledProb_ = true;
401 
402  //--- Go from microns to centimeters
403  theClusterParam.templXrec_ *= micronsToCm;
404  theClusterParam.templYrec_ *= micronsToCm;
405 
406  //--- Go back to the module coordinate system
407  theClusterParam.templXrec_ += lp.x();
408  theClusterParam.templYrec_ += lp.y();
409  }
410  return;
411 }
412 
413 //------------------------------------------------------------------
414 // Helper function to aggregate call & handling of Template 2D fit
415 //------------------------------------------------------------------
417  ClusterParamTemplate& theClusterParam,
418  SiPixelTemplateReco2D::ClusMatrix& clusterPayload,
419  int ID,
420  LocalPoint& lp) const {
422 
423  // Output:
424  float nonsense = -99999.9f; // nonsense init value
425  theClusterParam.templXrec_ = theClusterParam.templYrec_ = theClusterParam.templSigmaX_ =
426  theClusterParam.templSigmaY_ = nonsense;
427  // If the template recontruction fails, we want to return 1.0 for now
428  theClusterParam.probabilityX_ = theClusterParam.probabilityY_ = theClusterParam.probabilityQ_ = 1.f;
429  theClusterParam.qBin_ = 0;
430  // We have a boolean denoting whether the reco failed or not
431  theClusterParam.hasFilledProb_ = false;
432 
433  // ******************************************************************
434  //--- Call 2D TemplateReco
435  //
436  float locBz = theDetParam.bz;
437  float locBx = theDetParam.bx;
438 
439  //--- Input:
440  // edgeflagy - (input) flag to indicate the present of edges in y:
441  // 0=none (or interior gap),1=edge at small y, 2=edge at large y, 3=edge at either end
442  //
443  // edgeflagx - (input) flag to indicate the present of edges in x:
444  // 0=none, 1=edge at small x, 2=edge at large x
445  //
446  // These two variables are calculated in setTheClu() and stored in edgeTypeX_ and edgeTypeY_
447  //
448  //--- Output:
449  // deltay - (output) template y-length - cluster length [when > 0, possibly missing end]
450  // npixels - ??? &&& Ask Morris
451 
452  float deltay = 0; // return param
453  int npixels = 0; // return param
454 
455  //For now, turn off edgeX_ flags
456  theClusterParam.edgeTypeX_ = 0;
457 
458  if (clusterPayload.mrow > 4) {
459  // The cluster is too big, the 2D reco will perform horribly.
460  // Better to return immediately in error
461  theClusterParam.ierr2 = 8;
462 
463  } else {
464  theClusterParam.ierr2 = PixelTempReco2D(ID,
465  theClusterParam.cotalpha,
466  theClusterParam.cotbeta,
467  locBz,
468  locBx,
469  theClusterParam.edgeTypeY_,
470  theClusterParam.edgeTypeX_,
471  clusterPayload,
472  templ2d,
473  theClusterParam.templYrec_,
474  theClusterParam.templSigmaY_,
475  theClusterParam.templXrec_,
476  theClusterParam.templSigmaX_,
477  theClusterParam.templProbXY_,
478  theClusterParam.probabilityQ_,
479  theClusterParam.qBin_,
480  deltay,
481  npixels);
482  }
483  // ******************************************************************
484 
485  //--- Check exit status
486  if
487  UNLIKELY(theClusterParam.ierr2 != 0) {
488  LogDebug("PixelCPEClusterRepair::localPosition")
489  << "2D reconstruction failed with error " << theClusterParam.ierr2 << "\n";
490 
491  theClusterParam.probabilityX_ = theClusterParam.probabilityY_ = theClusterParam.probabilityQ_ = 0.f;
492  theClusterParam.qBin_ = 0;
493  // GG: what do we do in this case? For now, just return the cluster center of gravity in microns
494  // In the x case, apply a rough Lorentz drift average correction
495  float lorentz_drift = -999.9;
496  if (!GeomDetEnumerators::isEndcap(theDetParam.thePart))
497  lorentz_drift = 60.0f; // in microns // &&& replace with a constant (globally)
498  else
499  lorentz_drift = 10.0f; // in microns
500  // GG: trk angles needed to correct for bows/kinks
501  if (theClusterParam.with_track_angle) {
502  theClusterParam.templXrec_ =
503  theDetParam.theTopol->localX(theClusterParam.theCluster->x(), theClusterParam.loc_trk_pred) -
504  lorentz_drift * micronsToCm; // rough Lorentz drift correction
505  theClusterParam.templYrec_ =
506  theDetParam.theTopol->localY(theClusterParam.theCluster->y(), theClusterParam.loc_trk_pred);
507  } else {
508  edm::LogError("PixelCPEClusterRepair") << "@SUB = PixelCPEClusterRepair::localPosition"
509  << "Should never be here. PixelCPEClusterRepair should always be called "
510  "with track angles. This is a bad error !!! ";
511 
512  theClusterParam.templXrec_ = theDetParam.theTopol->localX(theClusterParam.theCluster->x()) -
513  lorentz_drift * micronsToCm; // rough Lorentz drift correction
514  theClusterParam.templYrec_ = theDetParam.theTopol->localY(theClusterParam.theCluster->y());
515  }
516  }
517  else {
518  //--- Template Reco succeeded.
519  theClusterParam.hasFilledProb_ = true;
520 
521  //--- Go from microns to centimeters
522  theClusterParam.templXrec_ *= micronsToCm;
523  theClusterParam.templYrec_ *= micronsToCm;
524 
525  //--- Go back to the module coordinate system
526  theClusterParam.templXrec_ += lp.x();
527  theClusterParam.templYrec_ += lp.y();
528  }
529  return;
530 }
531 //---------------------------------------------------------------------------------
532 // Helper function to see if we should recommend 2D reco to be run on this
533 // cluster
534 // ---------------------------------------------------------------------------------
536  ClusterParamTemplate& theClusterParam,
537  SiPixelTemplateReco::ClusMatrix& clusterPayload,
538  int ID) const
539 
540 {
541  // recommend2D is false by default
542  theClusterParam.recommended2D_ = false;
543 
544  DetId id = (theDetParam.theDet->geographicalId());
545 
546  bool recommend = false;
547  for (auto& rec : recommend2D_) {
548  recommend = rec.recommend(id, ttopo_);
549  if (recommend)
550  break;
551  }
552 
553  // only run on those layers recommended by configuration
554  if (!recommend) {
555  theClusterParam.recommended2D_ = false;
556  return;
557  }
558 
559  if (theClusterParam.edgeTypeY_) {
560  // For clusters that have edges in Y, run 2d reco.
561  // Flags already set by CPEBase
562  theClusterParam.recommended2D_ = true;
563  return;
564  }
565  // The 1d pixel template
567  if (!templ.interpolate(ID, theClusterParam.cotalpha, theClusterParam.cotbeta, theDetParam.bz, theDetParam.bx)) {
568  //error setting up template, return false
569  theClusterParam.recommended2D_ = false;
570  return;
571  }
572 
573  //length of the cluster taking into account double sized pixels
574  float nypix = clusterPayload.mcol;
575  for (int i = 0; i < clusterPayload.mcol; i++) {
576  if (clusterPayload.ydouble[i])
577  nypix += 1.;
578  }
579 
580  // templ.clsleny() is the expected length of the cluster along y axis.
581  // templ.qavg() is the expected total charge of the cluster
582  // theClusterParam.theCluster->charge() is the total charge of this cluster
583  float nydiff = templ.clsleny() - nypix;
584  float qratio = theClusterParam.theCluster->charge() / templ.qavg();
585 
586  if (nydiff > maxSizeMismatchInY_ && qratio < minChargeRatio_) {
587  // If the cluster is shorter than expected and has less charge, likely
588  // due to truncated cluster, try 2D reco
589 
590  theClusterParam.recommended2D_ = true;
591  theClusterParam.hasBadPixels_ = true;
592 
593  // if not RunDamagedClusters flag, don't try to fix any clusters
594  if (!runDamagedClusters_) {
595  theClusterParam.recommended2D_ = false;
596  }
597 
598  // Figure out what edge flags to set for truncated cluster
599  // Truncated clusters usually come from dead double columns
600  //
601  // If cluster is of even length, either both of or neither of beginning and ending
602  // edge are on a double column, so we cannot figure out the likely edge of
603  // truncation, let the 2D algorithm try extending on both sides (option 3)
604  if (theClusterParam.theCluster->sizeY() % 2 == 0)
605  theClusterParam.edgeTypeY_ = 3;
606  else {
607  //If the cluster is of odd length, only one of the edges can end on
608  //a double column, this is the likely edge of truncation
609  //Double columns always start on even indexes
610  int min_col = theClusterParam.theCluster->minPixelCol();
611  if (min_col % 2 == 0) {
612  //begining edge is at a double column (end edge cannot be,
613  //because odd length) so likely truncated at small y (option 1)
614  theClusterParam.edgeTypeY_ = 1;
615  } else {
616  //end edge is at a double column (beginning edge cannot be,
617  //because odd length) so likely truncated at large y (option 2)
618  theClusterParam.edgeTypeY_ = 2;
619  }
620  }
621  }
622 }
623 
624 //------------------------------------------------------------------
625 // localError() relies on localPosition() being called FIRST!!!
626 //------------------------------------------------------------------
627 LocalError PixelCPEClusterRepair::localError(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
628  ClusterParamTemplate& theClusterParam = static_cast<ClusterParamTemplate&>(theClusterParamBase);
629 
630  //--- Default is the maximum error used for edge clusters.
631  //--- (never used, in fact: let comment it out, shut up the complains of the static analyzer, and save a few CPU cycles)
632  float xerr = 0.0f, yerr = 0.0f;
633 
634  //--- Check status of both template calls.
635  if
636  UNLIKELY((theClusterParam.ierr != 0) || (theClusterParam.ierr2 != 0)) {
637  // If reconstruction fails the hit position is calculated from cluster center of gravity
638  // corrected in x by average Lorentz drift. Assign huge errors.
639  //
640  if
642  throw cms::Exception("PixelCPEClusterRepair::localPosition :") << "A non-pixel detector type in here?";
643 
644  // Assign better errors based on the residuals for failed template cases
645  if (GeomDetEnumerators::isBarrel(theDetParam.thePart)) {
646  xerr = 55.0f * micronsToCm; // &&& get errors from elsewhere?
647  yerr = 36.0f * micronsToCm;
648  } else {
649  xerr = 42.0f * micronsToCm;
650  yerr = 39.0f * micronsToCm;
651  }
652  }
653  // Leave commented for now, until we study the interplay of failure modes
654  // of 1D template reco and edges. For edge hits we run 2D reco by default!
655  //
656  // else if ( theClusterParam.edgeTypeX_ || theClusterParam.edgeTypeY_ ) {
657  // // for edge pixels assign errors according to observed residual RMS
658  // if ( theClusterParam.edgeTypeX_ && !theClusterParam.edgeTypeY_ ) {
659  // xerr = 23.0f * micronsToCm;
660  // yerr = 39.0f * micronsToCm;
661  // }
662  // else if ( !theClusterParam.edgeTypeX_ && theClusterParam.edgeTypeY_ ) {
663  // xerr = 24.0f * micronsToCm;
664  // yerr = 96.0f * micronsToCm;
665  // }
666  // else if ( theClusterParam.edgeTypeX_ && theClusterParam.edgeTypeY_ ) {
667  // xerr = 31.0f * micronsToCm;
668  // yerr = 90.0f * micronsToCm;
669  // }
670  // }
671  else {
672  xerr = theClusterParam.templSigmaX_ * micronsToCm;
673  yerr = theClusterParam.templSigmaY_ * micronsToCm;
674  // &&& should also check ierr (saved as class variable) and return
675  // &&& nonsense (another class static) if the template fit failed.
676  }
677 
678  if (theVerboseLevel > 9) {
679  LogDebug("PixelCPEClusterRepair") << " Sizex = " << theClusterParam.theCluster->sizeX()
680  << " Sizey = " << theClusterParam.theCluster->sizeY()
681  << " Edgex = " << theClusterParam.edgeTypeX_
682  << " Edgey = " << theClusterParam.edgeTypeY_ << " ErrX = " << xerr
683  << " ErrY = " << yerr;
684  }
685 
686  if (!(xerr > 0.0f))
687  throw cms::Exception("PixelCPEClusterRepair::localError")
688  << "\nERROR: Negative pixel error xerr = " << xerr << "\n\n";
689 
690  if (!(yerr > 0.0f))
691  throw cms::Exception("PixelCPEClusterRepair::localError")
692  << "\nERROR: Negative pixel error yerr = " << yerr << "\n\n";
693 
694  return LocalError(xerr * xerr, 0, yerr * yerr);
695 }
696 
698  static const boost::regex rule("([A-Z]+)(\\s+(\\d+))?");
699  boost::cmatch match;
700  // match and check it works
701  if (!regex_match(str.c_str(), match, rule))
702  throw cms::Exception("Configuration") << "Rule '" << str << "' not understood.\n";
703 
704  // subdet
705  subdet_ = -1;
706  if (strncmp(match[1].first, "PXB", 3) == 0)
708  else if (strncmp(match[1].first, "PXE", 3) == 0)
710  if (subdet_ == -1) {
711  throw cms::Exception("PixelCPEClusterRepair::Configuration")
712  << "Detector '" << match[1].first << "' not understood. Should be PXB, PXE.\n";
713  }
714  // layer (if present)
715  if (match[3].first != match[3].second) {
716  layer_ = atoi(match[3].first);
717  } else {
718  layer_ = 0;
719  }
720 } //end Rule::Rule
721 
723  desc.add<int>("barrelTemplateID", 0);
724  desc.add<int>("forwardTemplateID", 0);
725  desc.add<int>("directoryWithTemplates", 0);
726  desc.add<int>("speed", -2);
727  desc.add<bool>("UseClusterSplitter", false);
728  desc.add<double>("MaxSizeMismatchInY", 0.3);
729  desc.add<double>("MinChargeRatio", 0.8);
730  desc.add<std::vector<std::string>>("Recommend2D", {"PXB 2", "PXB 3", "PXB 4"});
731  desc.add<bool>("RunDamagedClusters", false);
732 }
#define LogDebug(id)
LocalError localError(DetParam const &theDetParam, ClusterParam &theClusterParam) const override
T getParameter(std::string const &) const
static const char layer_[]
int minPixelCol() const
Point3DBase< Scalar, LocalTag > LocalPoint
Definition: Definitions.h:30
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
const DetContainer & detUnits() const override
Returm a vector of all GeomDet.
int PixelTempReco2D(int id, float cotalpha, float cotbeta, float locBz, float locBx, int edgeflagy, int edgeflagx, ClusMatrix &cluster, SiPixelTemplate2D &templ, float &yrec, float &sigmay, float &xrec, float &sigmax, float &probxy, float &probQ, int &qbin, float &deltay, int &npixel)
int nominalValue() const
The nominal field value for this map in kGauss.
Definition: MagneticField.h:49
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
const SiPixelCluster * theCluster
Definition: PixelCPEBase.h:85
uint32_t ID
Definition: Definitions.h:24
bool isBarrel(GeomDetEnumerators::SubDetector m)
#define nullptr
std::vector< SiPixelTemplateStore2D > thePixelTemp2D_
T y() const
Definition: PV3DBase.h:60
DetParams m_DetParams
Definition: PixelCPEBase.h:280
int PixelTempReco1D(int id, float cotalpha, float cotbeta, float locBz, float locBx, ClusMatrix &cluster, SiPixelTemplate &templ, float &yrec, float &sigmay, float &proby, float &xrec, float &sigmax, float &probx, int &qbin, int speed, bool deadpix, std::vector< std::pair< int, int > > &zeropix, float &probQ, int &nypix, int &nxpix)
ClusterParam * createClusterParam(const SiPixelCluster &cl) const override
int charge() const
const PixelGeomDetUnit * theDet
Definition: PixelCPEBase.h:58
GeomDetType::SubDetector thePart
Definition: PixelCPEBase.h:63
std::vector< SiPixelTemplateStore > thePixelTemp_
const RectangularPixelTopology * theRecTopol
Definition: PixelCPEBase.h:61
std::vector< Rule > recommend2D_
U second(std::pair< T, U > const &p)
PixelCPEClusterRepair(edm::ParameterSet const &conf, const MagneticField *, const TrackerGeometry &, const TrackerTopology &, const SiPixelLorentzAngle *, const SiPixelTemplateDBObject *, const SiPixel2DTemplateDBObject *)
Measurement2DPoint MeasurementPoint
Measurement points are two-dimensional by default.
const PixelTopology * theTopol
Definition: PixelCPEBase.h:60
int minPixelRow() const
bool LoadTemplatesFromDB_
Definition: PixelCPEBase.h:254
const SiPixel2DTemplateDBObject * templateDBobject2D_
float clsleny()
y-size of smaller interpolated template in pixels
unsigned int offsetDU(SubDetector sid) const
const SiPixelTemplateDBObject * templateDBobject_
Definition: PixelCPEBase.h:250
virtual float localX(float mpX) const =0
DetId geographicalId() const
The label of this GeomDet.
Definition: GeomDet.h:64
double f[11][100]
static void fillPSetDescription(edm::ParameterSetDescription &desc)
T min(T a, T b)
Definition: MathUtil.h:58
void callTempReco2D(DetParam const &theDetParam, ClusterParamTemplate &theClusterParam, SiPixelTemplateReco2D::ClusMatrix &clusterPayload, int ID, LocalPoint &lp) const
ParameterDescriptionBase * add(U const &iLabel, T const &value)
SubDetector tkDetEnum[8]
bool isEndcap(GeomDetEnumerators::SubDetector m)
bool isItBigPixelInY(const int iybin) const override
bool isItBigPixelInX(const int ixbin) const override
Definition: DetId.h:17
Topology::LocalTrackPred loc_trk_pred
Definition: PixelCPEBase.h:99
short getTemplateID(const uint32_t &detid) const
const TrackerTopology & ttopo_
Definition: PixelCPEBase.h:242
void checkRecommend2D(DetParam const &theDetParam, ClusterParamTemplate &theClusterParam, SiPixelTemplateReco::ClusMatrix &clusterPayload, int ID) const
static bool pushfile(int filenum, std::vector< SiPixelTemplateStore2D > &pixelTemp, std::string dir="")
virtual LocalPoint localPosition(const MeasurementPoint &) const =0
float qavg()
average cluster charge for this set of track angles
void callTempReco1D(DetParam const &theDetParam, ClusterParamTemplate &theClusterParam, SiPixelTemplateReco::ClusMatrix &clusterPayload, int ID, LocalPoint &lp) const
static bool pushfile(int filenum, std::vector< SiPixelTemplateStore > &pixelTemp, std::string dir="CalibTracker/SiPixelESProducers/data/")
int sizeY() const
Pixel cluster – collection of neighboring pixels above threshold.
const TrackerGeometry & geom_
Definition: PixelCPEBase.h:241
Pixel pixel(int i) const
float y() const
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:10
virtual float localY(float mpY) const =0
#define UNLIKELY(x)
Definition: Likely.h:21
#define str(s)
int sizeX() const
T x() const
Definition: PV3DBase.h:59
float x() const
bool isTrackerPixel(GeomDetEnumerators::SubDetector m)
#define constexpr
LocalPoint localPosition(DetParam const &theDetParam, ClusterParam &theClusterParam) const override
int size() const