CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PixelCPEGeneric.cc
Go to the documentation of this file.
2 
5 
6 // this is needed to get errors from templates
10 
11 
12 // Services
15 
16 #include "boost/multi_array.hpp"
17 
18 #include <iostream>
19 using namespace std;
20 
21 const double HALF_PI = 1.57079632679489656;
22 
23 //-----------------------------------------------------------------------------
25 //-----------------------------------------------------------------------------
27  const MagneticField * mag, const SiPixelLorentzAngle * lorentzAngle, const SiPixelCPEGenericErrorParm * genErrorParm, const SiPixelTemplateDBObject * templateDBobject)
28  : PixelCPEBase(conf, mag, lorentzAngle, genErrorParm, templateDBobject)
29 {
30 
31  if (theVerboseLevel > 0)
32  LogDebug("PixelCPEGeneric")
33  << " constructing a generic algorithm for ideal pixel detector.\n"
34  << " CPEGeneric:: VerboseLevel = " << theVerboseLevel;
35 
36  // Externally settable cuts
37  the_eff_charge_cut_lowX = conf.getParameter<double>("eff_charge_cut_lowX");
38  the_eff_charge_cut_lowY = conf.getParameter<double>("eff_charge_cut_lowY");
39  the_eff_charge_cut_highX = conf.getParameter<double>("eff_charge_cut_highX");
40  the_eff_charge_cut_highY = conf.getParameter<double>("eff_charge_cut_highY");
41  the_size_cutX = conf.getParameter<double>("size_cutX");
42  the_size_cutY = conf.getParameter<double>("size_cutY");
43 
44  EdgeClusterErrorX_ = conf.getParameter<double>("EdgeClusterErrorX");
45  EdgeClusterErrorY_ = conf.getParameter<double>("EdgeClusterErrorY");
46 
47  // Externally settable flags to inflate errors
48  inflate_errors = conf.getParameter<bool>("inflate_errors");
49  inflate_all_errors_no_trk_angle = conf.getParameter<bool>("inflate_all_errors_no_trk_angle");
50 
51  UseErrorsFromTemplates_ = conf.getParameter<bool>("UseErrorsFromTemplates");
52  TruncatePixelCharge_ = conf.getParameter<bool>("TruncatePixelCharge");
53  IrradiationBiasCorrection_ = conf.getParameter<bool>("IrradiationBiasCorrection");
54  DoCosmics_ = conf.getParameter<bool>("DoCosmics");
55  LoadTemplatesFromDB_ = conf.getParameter<bool>("LoadTemplatesFromDB");
56 
57  if ( !UseErrorsFromTemplates_ && ( TruncatePixelCharge_ ||
58  IrradiationBiasCorrection_ ||
59  DoCosmics_ ||
60  LoadTemplatesFromDB_ ) )
61  {
62  throw cms::Exception("PixelCPEGeneric::PixelCPEGeneric: ")
63  << "\nERROR: UseErrorsFromTemplates_ is set to False in PixelCPEGeneric_cfi.py. "
64  << " In this case it does not make sense to set any of the following to True: "
65  << " TruncatePixelCharge_, IrradiationBiasCorrection_, DoCosmics_, LoadTemplatesFromDB_ !!!"
66  << "\n\n";
67  }
68 
69  if ( UseErrorsFromTemplates_ )
70  {
71  templID_ = -999;
72  if ( LoadTemplatesFromDB_ )
73  {
74  // Initialize template store to the selected ID [Morris, 6/25/08]
76  throw cms::Exception("InvalidCalibrationLoaded")
77  << "ERROR: Templates not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject version "
78  << ( *templateDBobject_ ).version() << ". Template ID is " << templID_;
79  }
80  else
81  {
82  if ( !templ_.pushfile( templID_ ) )
83  throw cms::Exception("InvalidCalibrationLoaded")
84  << "ERROR: Templates not loaded correctly from text file. Reconstruction will fail." << " Template ID is " << templID_;
85  }
86 
87  } // if ( UseErrorsFromTemplates_ )
88 
89  //cout << endl;
90  //cout << "From PixelCPEGeneric::PixelCPEGeneric(...)" << endl;
91  //cout << "(int)UseErrorsFromTemplates_ = " << (int)UseErrorsFromTemplates_ << endl;
92  //cout << "TruncatePixelCharge_ = " << (int)TruncatePixelCharge_ << endl;
93  //cout << "IrradiationBiasCorrection_ = " << (int)IrradiationBiasCorrection_ << endl;
94  //cout << "(int)DoCosmics_ = " << (int)DoCosmics_ << endl;
95  //cout << "(int)LoadTemplatesFromDB_ = " << (int)LoadTemplatesFromDB_ << endl;
96  //cout << endl;
97 
98  //yes, these should be config parameters!
99  //default case...
100  xerr_barrel_l1_= {0.00115, 0.00120, 0.00088};
101  xerr_barrel_l1_def_=0.01030;
102  yerr_barrel_l1_= {0.00375,0.00230,0.00250,0.00250,0.00230,0.00230,0.00210,0.00210,0.00240};
103  yerr_barrel_l1_def_=0.00210;
104  xerr_barrel_ln_= {0.00115, 0.00120, 0.00088};
105  xerr_barrel_ln_def_=0.01030;
106  yerr_barrel_ln_= {0.00375,0.00230,0.00250,0.00250,0.00230,0.00230,0.00210,0.00210,0.00240};
107  yerr_barrel_ln_def_=0.00210;
108  xerr_endcap_= {0.0020, 0.0020};
109  xerr_endcap_def_=0.0020;
110  yerr_endcap_= {0.00210};
111  yerr_endcap_def_=0.00075;
112 
113  bool isUpgrade=false;
114  if ( conf.exists("Upgrade") && conf.getParameter<bool>("Upgrade")) {
115  isUpgrade=true;
116  xerr_barrel_ln_= {0.00114,0.00104,0.00214};
117  xerr_barrel_ln_def_=0.00425;
118  yerr_barrel_ln_= {0.00299,0.00203,0.0023,0.00237,0.00233,0.00243,0.00232,0.00259,0.00176};
119  yerr_barrel_ln_def_=0.00245;
120  xerr_endcap_= {0.00151,0.000813,0.00221};
121  xerr_endcap_def_=0.00218;
122  yerr_endcap_= {0.00261,0.00107,0.00264};
123  yerr_endcap_def_=0.00357;
124 
125  if ( conf.exists("SmallPitch") && conf.getParameter<bool>("SmallPitch")) {
126  xerr_barrel_l1_= {0.00104, 0.000691, 0.00122};
127  xerr_barrel_l1_def_=0.00321;
128  yerr_barrel_l1_= {0.00199,0.00136,0.0015,0.00153,0.00152,0.00171,0.00154,0.00157,0.00154};
129  yerr_barrel_l1_def_=0.00164;
130  }
131  else{
132  xerr_barrel_l1_= {0.00114,0.00104,0.00214};
133  xerr_barrel_l1_def_=0.00425;
134  yerr_barrel_l1_= {0.00299,0.00203,0.0023,0.00237,0.00233,0.00243,0.00232,0.00259,0.00176};
135  yerr_barrel_l1_def_=0.00245;
136  }
137  }
138 
139  isUpgrade_=isUpgrade;
140 
141 }
142 
143 
144 
145 //-----------------------------------------------------------------------------
149 //-----------------------------------------------------------------------------
152  const GeomDetUnit & det) const
153 {
154  setTheDet( det, cluster );
157  {
159  /*bool fpix; //!< barrel(false) or forward(true)
160  if ( thePart == GeomDetEnumerators::PixelBarrel )
161  fpix = false; // no, it's not forward -- it's barrel
162  else
163  fpix = true; // yes, it's forward
164  */
165  float qclus = cluster.charge();
166 
167 
168  float locBz = (*theParam).bz;
169  //cout << "PixelCPEGeneric::localPosition(...) : locBz = " << locBz << endl;
170 
171  pixmx = -999.9; // max pixel charge for truncation of 2-D cluster
172  sigmay = -999.9; // CPE Generic y-error for multi-pixel cluster
173  deltay = -999.9; // CPE Generic y-bias for multi-pixel cluster
174  sigmax = -999.9; // CPE Generic x-error for multi-pixel cluster
175  deltax = -999.9; // CPE Generic x-bias for multi-pixel cluster
176  sy1 = -999.9; // CPE Generic y-error for single single-pixel
177  dy1 = -999.9; // CPE Generic y-bias for single single-pixel cluster
178  sy2 = -999.9; // CPE Generic y-error for single double-pixel cluster
179  dy2 = -999.9; // CPE Generic y-bias for single double-pixel cluster
180  sx1 = -999.9; // CPE Generic x-error for single single-pixel cluster
181  dx1 = -999.9; // CPE Generic x-bias for single single-pixel cluster
182  sx2 = -999.9; // CPE Generic x-error for single double-pixel cluster
183  dx2 = -999.9; // CPE Generic x-bias for single double-pixel cluster
184 
185  qBin_ = templ_.qbin( templID_, cotalpha_, cotbeta_, locBz, qclus, // inputs
186  pixmx, // returned by reference
187  sigmay, deltay, sigmax, deltax, // returned by reference
188  sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2 ); // returned by reference
189 
190  // These numbers come in microns from the qbin(...) call. Transform them to cm.
191  const float micronsToCm = 1.0e-4;
192 
194  dx1 = dx1 * micronsToCm;
195  dx2 = dx2 * micronsToCm;
196 
198  dy1 = dy1 * micronsToCm;
199  dy2 = dy2 * micronsToCm;
200 
202  sx1 = sx1 * micronsToCm;
203  sx2 = sx2 * micronsToCm;
204 
206  sy1 = sy1 * micronsToCm;
207  sy2 = sy2 * micronsToCm;
208 
209  } // if ( UseErrorsFromTemplates_ )
210 
211 
212  float Q_f_X = 0.0;
213  float Q_l_X = 0.0;
214  float Q_m_X = 0.0;
215  float Q_f_Y = 0.0;
216  float Q_l_Y = 0.0;
217  float Q_m_Y = 0.0;
218  collect_edge_charges( cluster,
219  Q_f_X, Q_l_X, Q_m_X,
220  Q_f_Y, Q_l_Y, Q_m_Y );
221 
222  //--- Find the inner widths along X and Y in one shot. We
223  //--- compute the upper right corner of the inner pixels
224  //--- (== lower left corner of upper right pixel) and
225  //--- the lower left corner of the inner pixels
226  //--- (== upper right corner of lower left pixel), and then
227  //--- subtract these two points in the formula.
228 
229  //--- Upper Right corner of Lower Left pixel -- in measurement frame
230  MeasurementPoint meas_URcorn_LLpix( cluster.minPixelRow()+1.0,
231  cluster.minPixelCol()+1.0 );
232 
233  //--- Lower Left corner of Upper Right pixel -- in measurement frame
234  MeasurementPoint meas_LLcorn_URpix( cluster.maxPixelRow(),
235  cluster.maxPixelCol() );
236 
237  //--- These two now converted into the local
238 
239  LocalPoint local_URcorn_LLpix;
240  LocalPoint local_LLcorn_URpix;
241 
242 
243  // PixelCPEGeneric can be used with or without track angles
244  // If PixelCPEGeneric is called with track angles, use them to correct for bows/kinks:
245  if ( with_track_angle )
246  {
247  local_URcorn_LLpix = theTopol->localPosition(meas_URcorn_LLpix, loc_trk_pred_);
248  local_LLcorn_URpix = theTopol->localPosition(meas_LLcorn_URpix, loc_trk_pred_);
249  }
250  else
251  {
252  local_URcorn_LLpix = theTopol->localPosition(meas_URcorn_LLpix);
253  local_LLcorn_URpix = theTopol->localPosition(meas_LLcorn_URpix);
254  }
255 
256  if (theVerboseLevel > 20) {
257  cout
258  << "\n\t >>> cluster.x = " << cluster.x()
259  << "\n\t >>> cluster.y = " << cluster.y()
260  << "\n\t >>> cluster: minRow = " << cluster.minPixelRow()
261  << " minCol = " << cluster.minPixelCol()
262  << "\n\t >>> cluster: maxRow = " << cluster.maxPixelRow()
263  << " maxCol = " << cluster.maxPixelCol()
264  << "\n\t >>> meas: inner lower left = " << meas_URcorn_LLpix.x()
265  << "," << meas_URcorn_LLpix.y()
266  << "\n\t >>> meas: inner upper right = " << meas_LLcorn_URpix.x()
267  << "," << meas_LLcorn_URpix.y()
268  << endl;
269  }
270 
271  //--- &&& Note that the cuts below should not be hardcoded (like in Orca and
272  //--- &&& CPEFromDetPosition/PixelCPEInitial), but rather be
273  //--- &&& externally settable (but tracked) parameters.
274 
275 
276  //--- Position, including the half lorentz shift
277  if (theVerboseLevel > 20)
278  cout << "\t >>> Generic:: processing X" << endl;
279 
280  float xPos =
281  generic_position_formula( cluster.sizeX(),
282  Q_f_X, Q_l_X,
283  local_URcorn_LLpix.x(), local_LLcorn_URpix.x(),
284  0.5*lorentzShiftInCmX_, // 0.5 * lorentz shift in
285  cotalpha_,
286  thePitchX,
287  theRecTopol->isItBigPixelInX( cluster.minPixelRow() ),
288  theRecTopol->isItBigPixelInX( cluster.maxPixelRow() ),
291  the_size_cutX); // cut for eff charge width &&&
292 
293 
294  if (theVerboseLevel > 20)
295  cout << "\t >>> Generic:: processing Y" << endl;
296  float yPos =
297  generic_position_formula( cluster.sizeY(),
298  Q_f_Y, Q_l_Y,
299  local_URcorn_LLpix.y(), local_LLcorn_URpix.y(),
300  0.5*lorentzShiftInCmY_, // 0.5 * lorentz shift in cm
301  cotbeta_,
302  thePitchY, // 0.5 * lorentz shift (may be 0)
303  theRecTopol->isItBigPixelInY( cluster.minPixelCol() ),
304  theRecTopol->isItBigPixelInY( cluster.maxPixelCol() ),
307  the_size_cutY); // cut for eff charge width &&&
308 
309  // Apply irradiation corrections.
311  {
312  if ( cluster.sizeX() == 1 )
313  {
314 
315  // Find if pixel is double (big).
316  bool bigInX = theRecTopol->isItBigPixelInX( cluster.maxPixelRow() );
317 
318  if ( !bigInX )
319  {
320  //cout << "Apply correction dx1 = " << dx1 << " to xPos = " << xPos << endl;
321  xPos -= dx1;
322  }
323  else
324  {
325  //cout << "Apply correction dx2 = " << dx2 << " to xPos = " << xPos << endl;
326  xPos -= dx2;
327  }
328  }
329  else
330  {
331  //cout << "Apply correction correction_deltax = " << deltax << " to xPos = " << xPos << endl;
332  xPos -= deltax;
333  }
334 
335  if ( cluster.sizeY() == 1 )
336  {
337 
338  // Find if pixel is double (big).
339  bool bigInY = theRecTopol->isItBigPixelInY( cluster.maxPixelCol() );
340 
341  if ( !bigInY )
342  {
343  //cout << "Apply correction dy1 = " << dy1 << " to yPos = " << yPos << endl;
344  yPos -= dy1;
345  }
346  else
347  {
348  //cout << "Apply correction dy2 = " << dy2 << " to yPos = " << yPos << endl;
349  yPos -= dy2;
350  }
351  }
352  else
353  {
354  //cout << "Apply correction deltay = " << deltay << " to yPos = " << yPos << endl;
355  yPos -= deltay;
356  }
357 
358  } // if ( IrradiationBiasCorrection_ )
359 
360  //--- Now put the two together
361  LocalPoint pos_in_local( xPos, yPos );
362  return pos_in_local;
363 }
364 
365 
366 
367 //-----------------------------------------------------------------------------
372 //-----------------------------------------------------------------------------
373 double
376  double Q_f,
377  double Q_l,
378  double upper_edge_first_pix,
379  double lower_edge_last_pix,
380  double half_lorentz_shift,
381  double cot_angle,
382  double pitch,
383  bool first_is_big,
384  bool last_is_big,
385  double eff_charge_cut_low,
386  double eff_charge_cut_high,
387  double size_cut
388  ) const
389 {
390  double geom_center = 0.5 * ( upper_edge_first_pix + lower_edge_last_pix );
391 
392  //--- The case of only one pixel in this projection is separate. Note that
393  //--- here first_pix == last_pix, so the average of the two is still the
394  //--- center of the pixel.
395  if ( size == 1 )
396  {
397  // ggiurgiu@jhu.edu, 02/03/09 : for size = 1, the Lorentz shift is already accounted by the irradiation correction
399  return geom_center;
400  else
401  return geom_center + half_lorentz_shift;
402  }
403 
404 
405  //--- Width of the clusters minus the edge (first and last) pixels.
406  //--- In the note, they are denoted x_F and x_L (and y_F and y_L)
407  double W_inner = lower_edge_last_pix - upper_edge_first_pix; // in cm
408 
409 
410  //--- Predicted charge width from geometry
411  double W_pred =
412  theThickness * cot_angle // geometric correction (in cm)
413  - 2 * half_lorentz_shift; // (in cm) &&& check fpix!
414 
415 
416  //--- Total length of the two edge pixels (first+last)
417  double sum_of_edge = 0.0;
418  if (first_is_big) sum_of_edge += 2.0;
419  else sum_of_edge += 1.0;
420 
421  if (last_is_big) sum_of_edge += 2.0;
422  else sum_of_edge += 1.0;
423 
424 
425  //--- The `effective' charge width -- particle's path in first and last pixels only
426  double W_eff = fabs( W_pred ) - W_inner;
427 
428 
429  //--- If the observed charge width is inconsistent with the expectations
430  //--- based on the track, do *not* use W_pred-W_innner. Instead, replace
431  //--- it with an *average* effective charge width, which is the average
432  //--- length of the edge pixels.
433  //
434  bool usedEdgeAlgo = false;
435  if (( W_eff/pitch < eff_charge_cut_low ) ||
436  ( W_eff/pitch > eff_charge_cut_high ) || (size >= size_cut))
437  {
438  W_eff = pitch * 0.5 * sum_of_edge; // ave. length of edge pixels (first+last) (cm)
439  usedEdgeAlgo = true;
441  }
442 
443 
444  //--- Finally, compute the position in this projection
445  double Qdiff = Q_l - Q_f;
446  double Qsum = Q_l + Q_f;
447 
448  //--- Temporary fix for clusters with both first and last pixel with charge = 0
449  if(Qsum==0) Qsum=1.0;
450  double hit_pos = geom_center + 0.5*(Qdiff/Qsum) * W_eff + half_lorentz_shift;
451 
452  //--- Debugging output
453  if (theVerboseLevel > 20) {
455  cout << "\t >>> We are in the Barrel." ;
456  } else {
457  cout << "\t >>> We are in the Forward." ;
458  }
459  cout
460  << "\n\t >>> cot(angle) = " << cot_angle << " pitch = " << pitch << " size = " << size
461  << "\n\t >>> upper_edge_first_pix = " << upper_edge_first_pix
462  << "\n\t >>> lower_edge_last_pix = " << lower_edge_last_pix
463  << "\n\t >>> geom_center = " << geom_center
464  << "\n\t >>> half_lorentz_shift = " << half_lorentz_shift
465  << "\n\t >>> W_inner = " << W_inner
466  << "\n\t >>> W_pred = " << W_pred
467  << "\n\t >>> W_eff(orig) = " << fabs( W_pred ) - W_inner
468  << "\n\t >>> W_eff(used) = " << W_eff
469  << "\n\t >>> sum_of_edge = " << sum_of_edge
470  << "\n\t >>> Qdiff = " << Qdiff << " Qsum = " << Qsum
471  << "\n\t >>> hit_pos = " << hit_pos
472  << "\n\t >>> RecHits: total = " << nRecHitsTotal_
473  << " used edge = " << nRecHitsUsedEdge_
474  << endl;
475  if (usedEdgeAlgo)
476  cout << "\n\t >>> Used Edge algorithm." ;
477  else
478  cout << "\n\t >>> Used angle information." ;
479  cout << endl;
480  }
481 
482 
483  return hit_pos;
484 }
485 
486 
487 
488 
489 
490 //-----------------------------------------------------------------------------
494 //-----------------------------------------------------------------------------
495 void
498  float & Q_f_X,
499  float & Q_l_X,
500  float & Q_m_X,
501  float & Q_f_Y,
502  float & Q_l_Y,
503  float & Q_m_Y
504  ) const
505 {
506  // Initialize return variables.
507  Q_f_X = Q_l_X = Q_m_X = 0.0;
508  Q_f_Y = Q_l_Y = Q_m_Y = 0.0;
509 
510 
511  // Obtain boundaries in index units
512  int xmin = cluster.minPixelRow();
513  int xmax = cluster.maxPixelRow();
514  int ymin = cluster.minPixelCol();
515  int ymax = cluster.maxPixelCol();
516 
517 
518 
519  // Iterate over the pixels.
520  int isize = cluster.size();
521 
522  for (int i = 0; i != isize; ++i)
523  {
524  auto const & pixel = cluster.pixel(i);
525  // ggiurgiu@fnal.gov: add pixel charge truncation
526  float pix_adc = pixel.adc;
528  pix_adc = std::min(pix_adc, pixmx );
529 
530  //
531  // X projection
532  if ( pixel.x == xmin ) // need to match with tolerance!!! &&&
533  Q_f_X += pix_adc;
534  else if ( pixel.x == xmax )
535  Q_l_X += pix_adc;
536  else
537  Q_m_X += pix_adc;
538  //
539  // Y projection
540  if ( pixel.y == ymin )
541  Q_f_Y += pix_adc;
542  else if ( pixel.y == ymax )
543  Q_l_Y += pix_adc;
544  else
545  Q_m_Y += pix_adc;
546  }
547 
548  return;
549 }
550 
551 
552 //============== INFLATED ERROR AND ERRORS FROM DB BELOW ================
553 
554 //-------------------------------------------------------------------------
555 // Hit error in the local frame
556 //-------------------------------------------------------------------------
557 LocalError
559  const GeomDetUnit & det) const
560 {
561  setTheDet( det, cluster );
562 
563  // The squared errors
564  float xerr_sq = -99999.9f;
565  float yerr_sq = -99999.9f;
566 
567 
568  // Default errors are the maximum error used for edge clusters.
569  /*
570  int row_offset = cluster.minPixelRow();
571  int col_offset = cluster.minPixelCol();
572  int n_bigInX = 0;
573  for (int irow = 0; irow < sizex; ++irow)
574  if ( theTopol->isItBigPixelInX( irow+row_offset ) )
575  ++n_bigInX;
576  int n_bigInY = 0;
577  for (int icol = 0; icol < sizey; ++icol)
578  if ( theTopol->isItBigPixelInY( icol+col_offset ) )
579  ++n_bigInX;
580  float xerr = (float)(sizex + n_bigInX) * thePitchX / sqrt(12.0);
581  float yerr = (float)(sizey + n_bigInY) * thePitchY / sqrt(12.0);
582  */
583 
584  // These are determined by looking at residuals for edge clusters
585  const float micronsToCm = 1.0e-4f;
586  float xerr = EdgeClusterErrorX_ * micronsToCm;
587  float yerr = EdgeClusterErrorY_ * micronsToCm;
588 
589  // Find if cluster is at the module edge.
590  int maxPixelCol = cluster.maxPixelCol();
591  int maxPixelRow = cluster.maxPixelRow();
592  int minPixelCol = cluster.minPixelCol();
593  int minPixelRow = cluster.minPixelRow();
594  unsigned int sizex = maxPixelRow - minPixelRow+1;
595  unsigned int sizey = maxPixelCol - minPixelCol+1;
596 
597  bool edgex = ( theRecTopol->isItEdgePixelInX( minPixelRow ) ) || ( theRecTopol->isItEdgePixelInX( maxPixelRow ) );
598  bool edgey = ( theRecTopol->isItEdgePixelInY( minPixelCol ) ) || ( theRecTopol->isItEdgePixelInY( maxPixelCol ) );
599 
600  // Find if cluster contains double (big) pixels.
601  bool bigInX = theRecTopol->containsBigPixelInX( minPixelRow, maxPixelRow );
602  bool bigInY = theRecTopol->containsBigPixelInY( minPixelCol, maxPixelCol );
603  if ( isUpgrade_ ||(!with_track_angle && DoCosmics_) )
604  {
605  //cout << "Track angles are not known and we are processing cosmics." << endl;
606  //cout << "Default angle estimation which assumes track from PV (0,0,0) does not work." << endl;
607  //cout << "Use an error parameterization which only depends on cluster size (by Vincenzo Chiochia)." << endl;
608 
610  {
611  DetId id = (det.geographicalId());
612  int layer=PXBDetId(id).layer();
613  if ( layer==1 ) {
614  if ( !edgex )
615  {
616  if ( sizex<=xerr_barrel_l1_.size() ) xerr=xerr_barrel_l1_[sizex-1];
617  else xerr=xerr_barrel_l1_def_;
618  }
619 
620  if ( !edgey )
621  {
622  if ( sizey<=yerr_barrel_l1_.size() ) yerr=yerr_barrel_l1_[sizey-1];
623  else yerr=yerr_barrel_l1_def_;
624  }
625  }
626  else{
627  if ( !edgex )
628  {
629  if ( sizex<=xerr_barrel_ln_.size() ) xerr=xerr_barrel_ln_[sizex-1];
630  else xerr=xerr_barrel_ln_def_;
631  }
632 
633  if ( !edgey )
634  {
635  if ( sizey<=yerr_barrel_ln_.size() ) yerr=yerr_barrel_ln_[sizey-1];
636  else yerr=yerr_barrel_ln_def_;
637  }
638  }
639  }
640  else // EndCap
641  {
642  if ( !edgex )
643  {
644  if ( sizex<=xerr_endcap_.size() ) xerr=xerr_endcap_[sizex-1];
645  else xerr=xerr_endcap_def_;
646  }
647 
648  if ( !edgey )
649  {
650  if ( sizey<=yerr_endcap_.size() ) yerr=yerr_endcap_[sizey-1];
651  else yerr=yerr_endcap_def_;
652  }
653  }
654 
655  } // if ( !with_track_angle )
656  else
657  {
658  //cout << "Track angles are known. We can use either errors from templates or the error parameterization from DB." << endl;
659 
661  {
662  if (qBin_ == 0 && inflate_errors )
663  {
664  int n_bigx = 0;
665  int n_bigy = 0;
666 
667  int row_offset = minPixelRow;
668  int col_offset = minPixelCol;
669 
670  for (int irow = 0; irow < 7; ++irow)
671  {
672  if ( theRecTopol->isItBigPixelInX( irow+row_offset ) )
673  ++n_bigx;
674  }
675 
676  for (int icol = 0; icol < 21; ++icol)
677  {
678  if ( theRecTopol->isItBigPixelInY( icol+col_offset ) )
679  ++n_bigy;
680  }
681 
682  xerr = (float)(sizex + n_bigx) * thePitchX / sqrt( 12.0f );
683  yerr = (float)(sizey + n_bigy) * thePitchY / sqrt( 12.0f );
684 
685  } // if ( qbin == 0 && inflate_errors )
686  else
687  {
688  // Default errors
689 
690  if ( !edgex )
691  {
692  if ( sizex == 1 )
693  {
694  if ( !bigInX )
695  xerr = sx1;
696  else
697  xerr = sx2;
698  }
699  else
700  xerr = sigmax;
701 
702  }
703 
704  if ( !edgey )
705  {
706  if ( sizey == 1 )
707  {
708  if ( !bigInY )
709  yerr = sy1;
710  else
711  yerr = sy2;
712  }
713  else
714  yerr = sigmay;
715 
716  }
717  } // if ( qbin == 0 && inflate_errors ) else
718 
719  } //if ( UseErrorsFromTemplates_ )
720  else
721  {
722  //cout << endl << "Use errors from DB:" << endl;
723 
724  if ( edgex && edgey )
725  {
726  //--- Both axes on the edge, no point in calling PixelErrorParameterization,
727  //--- just return the max errors on both.
728  }
729  else
730  {
731  pair<float,float> errPair =
733  alpha_, beta_, bigInX, bigInY );
734  if ( !edgex )
735  xerr = errPair.first;
736  if ( !edgey )
737  yerr = errPair.second;
738  }
739 
740  if (theVerboseLevel > 9)
741  {
742  LogDebug("PixelCPEGeneric") <<
743  " Sizex = " << cluster.sizeX() << " Sizey = " << cluster.sizeY() <<
744  " Edgex = " << edgex << " Edgey = " << edgey <<
745  " ErrX = " << xerr << " ErrY = " << yerr;
746  }
747 
748  } //if ( UseErrorsFromTemplates_ ) else
749 
750  } // if ( !with_track_angle ) else
751 
752  if ( !(xerr > 0.0) )
753  throw cms::Exception("PixelCPEGeneric::localError")
754  << "\nERROR: Negative pixel error xerr = " << xerr << "\n\n";
755 
756  if ( !(yerr > 0.0) )
757  throw cms::Exception("PixelCPEGeneric::localError")
758  << "\nERROR: Negative pixel error yerr = " << yerr << "\n\n";
759 
760  xerr_sq = xerr*xerr;
761  yerr_sq = yerr*yerr;
762 
763  return LocalError( xerr_sq, 0, yerr_sq );
764 
765 }
766 
767 
768 
769 
770 
771 
#define LogDebug(id)
T getParameter(std::string const &) const
std::vector< float > xerr_barrel_l1_
int i
Definition: DBlmapReader.cc:9
Topology::LocalTrackPred loc_trk_pred_
Definition: PixelCPEBase.h:276
float charge() const
int minPixelCol() const
int qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float &pixmx, float &sigmay, float &deltay, float &sigmax, float &deltax, float &sy1, float &dy1, float &sy2, float &dy2, float &sx1, float &dx1, float &sx2, float &dx2, float &lorywidth, float &lorxwidth)
void collect_edge_charges(const SiPixelCluster &cluster, float &Q_f_X, float &Q_l_X, float &Q_m_X, float &Q_f_Y, float &Q_l_Y, float &Q_m_Y) const
bool with_track_angle
Definition: PixelCPEBase.h:226
const float micronsToCm
double generic_position_formula(int size, double Q_f, double Q_l, double upper_edge_first_pix, double lower_edge_last_pix, double half_lorentz_shift, double cot_angle, double pitch, bool first_is_big, bool last_is_big, double eff_charge_cut_low, double eff_charge_cut_high, double size_cut) const
virtual LocalPoint localPosition(const MeasurementPoint &) const =0
SiPixelCPEGenericDBErrorParametrization * genErrorsFromDB_
std::vector< float > xerr_barrel_ln_
short getTemplateID(const uint32_t &detid) const
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
float theThickness
Definition: PixelCPEBase.h:196
T y() const
Definition: PV3DBase.h:63
bool exists(std::string const &parameterName) const
checks if a parameter exists
#define min(a, b)
Definition: mlp_lapack.h:161
std::pair< float, float > getError(const SiPixelCPEGenericErrorParm *parmErrors, GeomDetType::SubDetector pixelPart, int sizex, int sizey, float alpha, float beta, bool bigInX=false, bool bigInY=false)
GeomDetType::SubDetector thePart
Definition: PixelCPEBase.h:194
bool IrradiationBiasCorrection_
unsigned int layer() const
layer id
Definition: PXBDetId.h:35
std::vector< float > yerr_endcap_
double lorentzShiftInCmY_
Definition: PixelCPEBase.h:256
int maxPixelRow() const
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
double xPos
double the_eff_charge_cut_highX
int minPixelRow() const
LocalPoint localPosition(const SiPixelCluster &cluster, const GeomDetUnit &det) const
void computeLorentzShifts() const
double the_eff_charge_cut_lowY
T sqrt(T t)
Definition: SSEVec.h:48
PixelCPEGeneric(edm::ParameterSet const &conf, const MagneticField *, const SiPixelLorentzAngle *, const SiPixelCPEGenericErrorParm *, const SiPixelTemplateDBObject *)
The constructor.
bool inflate_all_errors_no_trk_angle
double lorentzShiftInCmX_
Definition: PixelCPEBase.h:255
const SiPixelTemplateDBObject * templateDBobject_
Definition: PixelCPEBase.h:269
double yPos
DetId geographicalId() const
The label of this GeomDet.
Definition: GeomDet.h:72
double f[11][100]
std::vector< float > xerr_endcap_
LocalError localError(const SiPixelCluster &cl, const GeomDetUnit &det) const
SiPixelTemplate templ_
tuple conf
Definition: dbtoconf.py:185
Definition: DetId.h:20
bool pushfile(int filenum)
const RectangularPixelTopology * theRecTopol
Definition: PixelCPEBase.h:190
const PixelTopology * theTopol
Definition: PixelCPEBase.h:189
int maxPixelCol() const
std::vector< float > yerr_barrel_l1_
std::vector< float > yerr_barrel_ln_
int sizeY() const
void setTheDet(const GeomDetUnit &det, const SiPixelCluster &cluster) const
Definition: PixelCPEBase.cc:89
Pixel cluster – collection of neighboring pixels above threshold.
Pixel pixel(int i) const
double the_eff_charge_cut_lowX
float y() const
tuple cout
Definition: gather_cfg.py:121
const PixelGeomDetUnit * theDet
Definition: PixelCPEBase.h:185
double the_eff_charge_cut_highY
const double HALF_PI
int sizeX() const
T x() const
Definition: PV3DBase.h:62
tuple size
Write out results.
float x() const
const SiPixelCPEGenericErrorParm * genErrorParm_
Definition: PixelCPEBase.h:267
int nRecHitsUsedEdge_
Definition: PixelCPEBase.h:223
int size() const