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
9 
10 // Services
13 
14 #include "boost/multi_array.hpp"
15 
16 #include <iostream>
17 using namespace std;
18 
19 namespace {
20  constexpr float micronsToCm = 1.0e-4;
21  const bool MYDEBUG = false;
22 } // namespace
23 
24 //-----------------------------------------------------------------------------
26 //-----------------------------------------------------------------------------
28  const MagneticField* mag,
29  const TrackerGeometry& geom,
30  const TrackerTopology& ttopo,
31  const SiPixelLorentzAngle* lorentzAngle,
32  const SiPixelGenErrorDBObject* genErrorDBObject,
33  const SiPixelLorentzAngle* lorentzAngleWidth = nullptr)
34  : PixelCPEBase(conf, mag, geom, ttopo, lorentzAngle, genErrorDBObject, nullptr, lorentzAngleWidth, 0) {
35  if (theVerboseLevel > 0)
36  LogDebug("PixelCPEGeneric") << " constructing a generic algorithm for ideal pixel detector.\n"
37  << " CPEGeneric:: VerboseLevel = " << theVerboseLevel;
38 
39  // Externally settable cuts
40  the_eff_charge_cut_lowX = conf.getParameter<double>("eff_charge_cut_lowX");
41  the_eff_charge_cut_lowY = conf.getParameter<double>("eff_charge_cut_lowY");
42  the_eff_charge_cut_highX = conf.getParameter<double>("eff_charge_cut_highX");
43  the_eff_charge_cut_highY = conf.getParameter<double>("eff_charge_cut_highY");
44  the_size_cutX = conf.getParameter<double>("size_cutX");
45  the_size_cutY = conf.getParameter<double>("size_cutY");
46 
47  EdgeClusterErrorX_ = conf.getParameter<double>("EdgeClusterErrorX");
48  EdgeClusterErrorY_ = conf.getParameter<double>("EdgeClusterErrorY");
49 
50  // Externally settable flags to inflate errors
51  inflate_errors = conf.getParameter<bool>("inflate_errors");
52  inflate_all_errors_no_trk_angle = conf.getParameter<bool>("inflate_all_errors_no_trk_angle");
53 
54  UseErrorsFromTemplates_ = conf.getParameter<bool>("UseErrorsFromTemplates");
55  TruncatePixelCharge_ = conf.getParameter<bool>("TruncatePixelCharge");
56  IrradiationBiasCorrection_ = conf.getParameter<bool>("IrradiationBiasCorrection");
57  DoCosmics_ = conf.getParameter<bool>("DoCosmics");
58  //LoadTemplatesFromDB_ = conf.getParameter<bool>("LoadTemplatesFromDB");
59 
60  // no clear what upgrade means, is it phase1, phase2? Probably delete.
61  isUpgrade_ = false;
62  if (conf.getParameter<bool>("Upgrade"))
63  isUpgrade_ = true;
64 
65  // Select the position error source
66  // For upgrde and cosmics force the use simple errors
67  if (isUpgrade_ || (DoCosmics_))
68  UseErrorsFromTemplates_ = false;
69 
70  if (!UseErrorsFromTemplates_ && (TruncatePixelCharge_ || IrradiationBiasCorrection_ || LoadTemplatesFromDB_)) {
71  throw cms::Exception("PixelCPEGeneric::PixelCPEGeneric: ")
72  << "\nERROR: UseErrorsFromTemplates_ is set to False in PixelCPEGeneric_cfi.py. "
73  << " In this case it does not make sense to set any of the following to True: "
74  << " TruncatePixelCharge_, IrradiationBiasCorrection_, DoCosmics_, LoadTemplatesFromDB_ !!!"
75  << "\n\n";
76  }
77 
78  // Use errors from templates or from GenError
79  if (UseErrorsFromTemplates_) {
80  if (LoadTemplatesFromDB_) { // From DB
82  throw cms::Exception("InvalidCalibrationLoaded")
83  << "ERROR: GenErrors not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject version "
84  << (*genErrorDBObject_).version();
85  if (MYDEBUG)
86  cout << "Loaded genErrorDBObject v" << (*genErrorDBObject_).version() << endl;
87  } else { // From file
89  throw cms::Exception("InvalidCalibrationLoaded")
90  << "ERROR: GenErrors not loaded correctly from text file. Reconstruction will fail.";
91  } // if load from DB
92 
93  } else {
94  if (MYDEBUG)
95  cout << " Use simple parametrised errors " << endl;
96  } // if ( UseErrorsFromTemplates_ )
97 
98  // Rechit errors in case other, more correct, errors are not vailable
99  // This are constants. Maybe there is a more efficienct way to store them.
100  if (!isUpgrade_) { // normal case
101  xerr_barrel_l1_ = {0.00115, 0.00120, 0.00088};
102  xerr_barrel_l1_def_ = 0.01030;
103  yerr_barrel_l1_ = {0.00375, 0.00230, 0.00250, 0.00250, 0.00230, 0.00230, 0.00210, 0.00210, 0.00240};
104  yerr_barrel_l1_def_ = 0.00210;
105  xerr_barrel_ln_ = {0.00115, 0.00120, 0.00088};
106  xerr_barrel_ln_def_ = 0.01030;
107  yerr_barrel_ln_ = {0.00375, 0.00230, 0.00250, 0.00250, 0.00230, 0.00230, 0.00210, 0.00210, 0.00240};
108  yerr_barrel_ln_def_ = 0.00210;
109  xerr_endcap_ = {0.0020, 0.0020};
110  xerr_endcap_def_ = 0.0020;
111  yerr_endcap_ = {0.00210};
112  yerr_endcap_def_ = 0.00075;
113  } else { // isUpgrade=true
114  xerr_barrel_ln_ = {0.00025, 0.00030, 0.00035, 0.00035};
115  xerr_barrel_ln_def_ = 0.00035;
116  yerr_barrel_ln_ = {0.00210, 0.00115, 0.00125};
117  yerr_barrel_ln_def_ = 0.00125;
118  xerr_endcap_ = {0.00072, 0.00025};
119  xerr_endcap_def_ = 0.00060;
120  yerr_endcap_ = {0.00289, 0.00025};
121  yerr_endcap_def_ = 0.00180;
122 
123  if (conf.getParameter<bool>("SmallPitch")) {
124  xerr_barrel_l1_ = {0.00104, 0.000691, 0.00122};
125  xerr_barrel_l1_def_ = 0.00321;
126  yerr_barrel_l1_ = {0.00199, 0.00136, 0.0015, 0.00153, 0.00152, 0.00171, 0.00154, 0.00157, 0.00154};
127  yerr_barrel_l1_def_ = 0.00164;
128  } else {
129  xerr_barrel_l1_ = {0.00025, 0.00030, 0.00035, 0.00035};
130  xerr_barrel_l1_def_ = 0.00035;
131  yerr_barrel_l1_ = {0.00210, 0.00115, 0.00125};
132  yerr_barrel_l1_def_ = 0.00125;
133  }
134  } // if isUpgrade
135 
136  if (MYDEBUG) {
137  cout << "From PixelCPEGeneric::PixelCPEGeneric(...)" << endl;
138  cout << "(int)UseErrorsFromTemplates_ = " << (int)UseErrorsFromTemplates_ << endl;
139  cout << "TruncatePixelCharge_ = " << (int)TruncatePixelCharge_ << endl;
140  cout << "IrradiationBiasCorrection_ = " << (int)IrradiationBiasCorrection_ << endl;
141  cout << "(int)DoCosmics_ = " << (int)DoCosmics_ << endl;
142  cout << "(int)LoadTemplatesFromDB_ = " << (int)LoadTemplatesFromDB_ << endl;
143  }
144 }
145 
147  return new ClusterParamGeneric(cl);
148 }
149 
150 //-----------------------------------------------------------------------------
154 //-----------------------------------------------------------------------------
155 LocalPoint PixelCPEGeneric::localPosition(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
156  ClusterParamGeneric& theClusterParam = static_cast<ClusterParamGeneric&>(theClusterParamBase);
157 
158  //cout<<" in PixelCPEGeneric:localPosition - "<<endl; //dk
159 
160  float chargeWidthX = (theDetParam.lorentzShiftInCmX * theDetParam.widthLAFractionX);
161  float chargeWidthY = (theDetParam.lorentzShiftInCmY * theDetParam.widthLAFractionY);
162  float shiftX = 0.5f * theDetParam.lorentzShiftInCmX;
163  float shiftY = 0.5f * theDetParam.lorentzShiftInCmY;
164 
165  //cout<<" main la width "<<chargeWidthX<<" "<<chargeWidthY<<endl;
166 
168  float qclus = theClusterParam.theCluster->charge();
169  float locBz = theDetParam.bz;
170  float locBx = theDetParam.bx;
171  //cout << "PixelCPEGeneric::localPosition(...) : locBz = " << locBz << endl;
172 
173  theClusterParam.pixmx = -999; // max pixel charge for truncation of 2-D cluster
174  theClusterParam.sigmay = -999.9; // CPE Generic y-error for multi-pixel cluster
175  theClusterParam.deltay = -999.9; // CPE Generic y-bias for multi-pixel cluster
176  theClusterParam.sigmax = -999.9; // CPE Generic x-error for multi-pixel cluster
177  theClusterParam.deltax = -999.9; // CPE Generic x-bias for multi-pixel cluster
178  theClusterParam.sy1 = -999.9; // CPE Generic y-error for single single-pixel
179  theClusterParam.dy1 = -999.9; // CPE Generic y-bias for single single-pixel cluster
180  theClusterParam.sy2 = -999.9; // CPE Generic y-error for single double-pixel cluster
181  theClusterParam.dy2 = -999.9; // CPE Generic y-bias for single double-pixel cluster
182  theClusterParam.sx1 = -999.9; // CPE Generic x-error for single single-pixel cluster
183  theClusterParam.dx1 = -999.9; // CPE Generic x-bias for single single-pixel cluster
184  theClusterParam.sx2 = -999.9; // CPE Generic x-error for single double-pixel cluster
185  theClusterParam.dx2 = -999.9; // CPE Generic x-bias for single double-pixel cluster
186 
188  int gtemplID_ = theDetParam.detTemplateId;
189 
190  //int gtemplID0 = genErrorDBObject_->getGenErrorID(theDetParam.theDet->geographicalId().rawId());
191  //if(gtemplID0!=gtemplID_) cout<<" different id "<< gtemplID_<<" "<<gtemplID0<<endl;
192 
193  theClusterParam.qBin_ = gtempl.qbin(gtemplID_,
194  theClusterParam.cotalpha,
195  theClusterParam.cotbeta,
196  locBz,
197  locBx,
198  qclus,
200  theClusterParam.pixmx,
201  theClusterParam.sigmay,
202  theClusterParam.deltay,
203  theClusterParam.sigmax,
204  theClusterParam.deltax,
205  theClusterParam.sy1,
206  theClusterParam.dy1,
207  theClusterParam.sy2,
208  theClusterParam.dy2,
209  theClusterParam.sx1,
210  theClusterParam.dx1,
211  theClusterParam.sx2,
212  theClusterParam.dx2);
213 
214  // now use the charge widths stored in the new generic template headers (change to the
215  // incorrect sign convention of the base class)
216  bool useLAWidthFromGenError = false;
217  if (useLAWidthFromGenError) {
218  chargeWidthX = (-micronsToCm * gtempl.lorxwidth());
219  chargeWidthY = (-micronsToCm * gtempl.lorywidth());
220  if (MYDEBUG)
221  cout << " redefine la width (gen-error) " << chargeWidthX << " " << chargeWidthY << endl;
222  }
223  if (MYDEBUG)
224  cout << " GenError: " << gtemplID_ << endl;
225 
226  // These numbers come in microns from the qbin(...) call. Transform them to cm.
227  theClusterParam.deltax = theClusterParam.deltax * micronsToCm;
228  theClusterParam.dx1 = theClusterParam.dx1 * micronsToCm;
229  theClusterParam.dx2 = theClusterParam.dx2 * micronsToCm;
230 
231  theClusterParam.deltay = theClusterParam.deltay * micronsToCm;
232  theClusterParam.dy1 = theClusterParam.dy1 * micronsToCm;
233  theClusterParam.dy2 = theClusterParam.dy2 * micronsToCm;
234 
235  theClusterParam.sigmax = theClusterParam.sigmax * micronsToCm;
236  theClusterParam.sx1 = theClusterParam.sx1 * micronsToCm;
237  theClusterParam.sx2 = theClusterParam.sx2 * micronsToCm;
238 
239  theClusterParam.sigmay = theClusterParam.sigmay * micronsToCm;
240  theClusterParam.sy1 = theClusterParam.sy1 * micronsToCm;
241  theClusterParam.sy2 = theClusterParam.sy2 * micronsToCm;
242 
243  } // if ( UseErrorsFromTemplates_ )
244  else {
245  theClusterParam.qBin_ = 0;
246  }
247 
248  int Q_f_X;
249  int Q_l_X;
250  int Q_f_Y;
251  int Q_l_Y;
252  collect_edge_charges(theClusterParam, Q_f_X, Q_l_X, Q_f_Y, Q_l_Y);
253 
254  //--- Find the inner widths along X and Y in one shot. We
255  //--- compute the upper right corner of the inner pixels
256  //--- (== lower left corner of upper right pixel) and
257  //--- the lower left corner of the inner pixels
258  //--- (== upper right corner of lower left pixel), and then
259  //--- subtract these two points in the formula.
260 
261  //--- Upper Right corner of Lower Left pixel -- in measurement frame
262  MeasurementPoint meas_URcorn_LLpix(theClusterParam.theCluster->minPixelRow() + 1.0,
263  theClusterParam.theCluster->minPixelCol() + 1.0);
264 
265  //--- Lower Left corner of Upper Right pixel -- in measurement frame
266  MeasurementPoint meas_LLcorn_URpix(theClusterParam.theCluster->maxPixelRow(),
267  theClusterParam.theCluster->maxPixelCol());
268 
269  //--- These two now converted into the local
270  LocalPoint local_URcorn_LLpix;
271  LocalPoint local_LLcorn_URpix;
272 
273  // PixelCPEGeneric can be used with or without track angles
274  // If PixelCPEGeneric is called with track angles, use them to correct for bows/kinks:
275  if (theClusterParam.with_track_angle) {
276  local_URcorn_LLpix = theDetParam.theTopol->localPosition(meas_URcorn_LLpix, theClusterParam.loc_trk_pred);
277  local_LLcorn_URpix = theDetParam.theTopol->localPosition(meas_LLcorn_URpix, theClusterParam.loc_trk_pred);
278  } else {
279  local_URcorn_LLpix = theDetParam.theTopol->localPosition(meas_URcorn_LLpix);
280  local_LLcorn_URpix = theDetParam.theTopol->localPosition(meas_LLcorn_URpix);
281  }
282 
283 #ifdef EDM_ML_DEBUG
284  if (theVerboseLevel > 20) {
285  cout << "\n\t >>> theClusterParam.theCluster->x = " << theClusterParam.theCluster->x()
286  << "\n\t >>> theClusterParam.theCluster->y = " << theClusterParam.theCluster->y()
287  << "\n\t >>> cluster: minRow = " << theClusterParam.theCluster->minPixelRow()
288  << " minCol = " << theClusterParam.theCluster->minPixelCol()
289  << "\n\t >>> cluster: maxRow = " << theClusterParam.theCluster->maxPixelRow()
290  << " maxCol = " << theClusterParam.theCluster->maxPixelCol()
291  << "\n\t >>> meas: inner lower left = " << meas_URcorn_LLpix.x() << "," << meas_URcorn_LLpix.y()
292  << "\n\t >>> meas: inner upper right = " << meas_LLcorn_URpix.x() << "," << meas_LLcorn_URpix.y() << endl;
293  }
294 #endif
295 
296  //--- &&& Note that the cuts below should not be hardcoded (like in Orca and
297  //--- &&& CPEFromDetPosition/PixelCPEInitial), but rather be
298  //--- &&& externally settable (but tracked) parameters.
299 
300  //--- Position, including the half lorentz shift
301 
302 #ifdef EDM_ML_DEBUG
303  if (theVerboseLevel > 20)
304  cout << "\t >>> Generic:: processing X" << endl;
305 #endif
306 
307  float xPos =
308  generic_position_formula(theClusterParam.theCluster->sizeX(),
309  Q_f_X,
310  Q_l_X,
311  local_URcorn_LLpix.x(),
312  local_LLcorn_URpix.x(),
313  chargeWidthX, // lorentz shift in cm
314  theDetParam.theThickness,
315  theClusterParam.cotalpha,
316  theDetParam.thePitchX,
317  theDetParam.theRecTopol->isItBigPixelInX(theClusterParam.theCluster->minPixelRow()),
318  theDetParam.theRecTopol->isItBigPixelInX(theClusterParam.theCluster->maxPixelRow()),
321  the_size_cutX); // cut for eff charge width &&&
322 
323  // apply the lorentz offset correction
324  xPos = xPos + shiftX;
325 
326 #ifdef EDM_ML_DEBUG
327  if (theVerboseLevel > 20)
328  cout << "\t >>> Generic:: processing Y" << endl;
329 #endif
330 
331  float yPos =
332  generic_position_formula(theClusterParam.theCluster->sizeY(),
333  Q_f_Y,
334  Q_l_Y,
335  local_URcorn_LLpix.y(),
336  local_LLcorn_URpix.y(),
337  chargeWidthY, // lorentz shift in cm
338  theDetParam.theThickness,
339  theClusterParam.cotbeta,
340  theDetParam.thePitchY,
341  theDetParam.theRecTopol->isItBigPixelInY(theClusterParam.theCluster->minPixelCol()),
342  theDetParam.theRecTopol->isItBigPixelInY(theClusterParam.theCluster->maxPixelCol()),
345  the_size_cutY); // cut for eff charge width &&&
346 
347  // apply the lorentz offset correction
348  yPos = yPos + shiftY;
349 
350  // Apply irradiation corrections. NOT USED FOR NOW
352  if (theClusterParam.theCluster->sizeX() == 1) { // size=1
353  // ggiurgiu@jhu.edu, 02/03/09 : for size = 1, the Lorentz shift is already accounted by the irradiation correction
354  //float tmp1 = (0.5 * theDetParam.lorentzShiftInCmX);
355  //cout << "Apply correction correction_dx1 = " << theClusterParam.dx1 << " to xPos = " << xPos;
356  xPos = xPos - (0.5f * theDetParam.lorentzShiftInCmX);
357  // Find if pixel is double (big).
358  bool bigInX = theDetParam.theRecTopol->isItBigPixelInX(theClusterParam.theCluster->maxPixelRow());
359  if (!bigInX)
360  xPos -= theClusterParam.dx1;
361  else
362  xPos -= theClusterParam.dx2;
363  //cout<<" to "<<xPos<<" "<<(tmp1+theClusterParam.dx1)<<endl;
364  } else { // size>1
365  //cout << "Apply correction correction_deltax = " << theClusterParam.deltax << " to xPos = " << xPos;
366  xPos -= theClusterParam.deltax;
367  //cout<<" to "<<xPos<<endl;
368  }
369 
370  if (theClusterParam.theCluster->sizeY() == 1) {
371  // ggiurgiu@jhu.edu, 02/03/09 : for size = 1, the Lorentz shift is already accounted by the irradiation correction
372  yPos = yPos - (0.5f * theDetParam.lorentzShiftInCmY);
373 
374  // Find if pixel is double (big).
375  bool bigInY = theDetParam.theRecTopol->isItBigPixelInY(theClusterParam.theCluster->maxPixelCol());
376  if (!bigInY)
377  yPos -= theClusterParam.dy1;
378  else
379  yPos -= theClusterParam.dy2;
380 
381  } else {
382  //cout << "Apply correction correction_deltay = " << theClusterParam.deltay << " to yPos = " << yPos << endl;
383  yPos -= theClusterParam.deltay;
384  }
385 
386  } // if ( IrradiationBiasCorrection_ )
387 
388  //cout<<" in PixelCPEGeneric:localPosition - pos = "<<xPos<<" "<<yPos<<endl; //dk
389 
390  //--- Now put the two together
391  LocalPoint pos_in_local(xPos, yPos);
392  return pos_in_local;
393 }
394 
395 //-----------------------------------------------------------------------------
400 //-----------------------------------------------------------------------------
402  int Q_f,
403  int Q_l,
404  float upper_edge_first_pix,
405  float lower_edge_last_pix,
406  float lorentz_shift,
407  float theThickness, //detector thickness
408  float cot_angle,
409  float pitch,
410  bool first_is_big,
411  bool last_is_big,
412  float eff_charge_cut_low,
413  float eff_charge_cut_high,
414  float size_cut
415  ) const {
416  //cout<<" in PixelCPEGeneric:generic_position_formula - "<<endl; //dk
417 
418  float geom_center = 0.5f * (upper_edge_first_pix + lower_edge_last_pix);
419 
420  //--- The case of only one pixel in this projection is separate. Note that
421  //--- here first_pix == last_pix, so the average of the two is still the
422  //--- center of the pixel.
423  if (size == 1) {
424  return geom_center;
425  }
426 
427  //--- Width of the clusters minus the edge (first and last) pixels.
428  //--- In the note, they are denoted x_F and x_L (and y_F and y_L)
429  float W_inner = lower_edge_last_pix - upper_edge_first_pix; // in cm
430 
431  //--- Predicted charge width from geometry
432  float W_pred = theThickness * cot_angle // geometric correction (in cm)
433  - lorentz_shift; // (in cm) &&& check fpix!
434 
435  //cout<<" in PixelCPEGeneric:generic_position_formula - "<<W_inner<<" "<<W_pred<<endl; //dk
436 
437  //--- Total length of the two edge pixels (first+last)
438  float sum_of_edge = 2.0f;
439  if (first_is_big)
440  sum_of_edge += 1.0f;
441  if (last_is_big)
442  sum_of_edge += 1.0f;
443 
444  //--- The `effective' charge width -- particle's path in first and last pixels only
445  float W_eff = std::abs(W_pred) - W_inner;
446 
447  //--- If the observed charge width is inconsistent with the expectations
448  //--- based on the track, do *not* use W_pred-W_innner. Instead, replace
449  //--- it with an *average* effective charge width, which is the average
450  //--- length of the edge pixels.
451  //
452  // bool usedEdgeAlgo = false;
453  if ((size >= size_cut) || ((W_eff / pitch < eff_charge_cut_low) | (W_eff / pitch > eff_charge_cut_high))) {
454  W_eff = pitch * 0.5f * sum_of_edge; // ave. length of edge pixels (first+last) (cm)
455  // usedEdgeAlgo = true;
456 #ifdef EDM_ML_DEBUG
457  nRecHitsUsedEdge_++;
458 #endif
459  }
460 
461  //--- Finally, compute the position in this projection
462  float Qdiff = Q_l - Q_f;
463  float Qsum = Q_l + Q_f;
464 
465  //--- Temporary fix for clusters with both first and last pixel with charge = 0
466  if (Qsum == 0)
467  Qsum = 1.0f;
468  //float hit_pos = geom_center + 0.5f*(Qdiff/Qsum) * W_eff + half_lorentz_shift;
469  float hit_pos = geom_center + 0.5f * (Qdiff / Qsum) * W_eff;
470 
471  //cout<<" in PixelCPEGeneric:generic_position_formula - "<<hit_pos<<" "<<lorentz_shift*0.5<<endl; //dk
472 
473 #ifdef EDM_ML_DEBUG
474  //--- Debugging output
475 #warning "Debug printouts in PixelCPEGeneric.cc has been commented because they cannot be compiled"
476  /* This part is commented because some variables used here are not defined !!
477  if (theVerboseLevel > 20) {
478  if ( theDetParam.thePart == GeomDetEnumerators::PixelBarrel || theDetParam.thePart == GeomDetEnumerators::P1PXB ) {
479  cout << "\t >>> We are in the Barrel." ;
480  } else if ( theDetParam.thePart == GeomDetEnumerators::PixelEndcap ||
481  theDetParam.thePart == GeomDetEnumerators::P1PXEC ||
482  theDetParam.thePart == GeomDetEnumerators::P2PXEC ) {
483  cout << "\t >>> We are in the Forward." ;
484  } else {
485  cout << "\t >>> We are in an unexpected subdet " << theDetParam.thePart;
486  }
487  cout
488  << "\n\t >>> cot(angle) = " << cot_angle << " pitch = " << pitch << " size = " << size
489  << "\n\t >>> upper_edge_first_pix = " << upper_edge_first_pix
490  << "\n\t >>> lower_edge_last_pix = " << lower_edge_last_pix
491  << "\n\t >>> geom_center = " << geom_center
492  << "\n\t >>> half_lorentz_shift = " << half_lorentz_shift
493  << "\n\t >>> W_inner = " << W_inner
494  << "\n\t >>> W_pred = " << W_pred
495  << "\n\t >>> W_eff(orig) = " << fabs( W_pred ) - W_inner
496  << "\n\t >>> W_eff(used) = " << W_eff
497  << "\n\t >>> sum_of_edge = " << sum_of_edge
498  << "\n\t >>> Qdiff = " << Qdiff << " Qsum = " << Qsum
499  << "\n\t >>> hit_pos = " << hit_pos
500  << "\n\t >>> RecHits: total = " << nRecHitsTotal_
501  << " used edge = " << nRecHitsUsedEdge_
502  << endl;
503  if (usedEdgeAlgo)
504  cout << "\n\t >>> Used Edge algorithm." ;
505  else
506  cout << "\n\t >>> Used angle information." ;
507  cout << endl;
508  }
509  */
510 #endif
511 
512  return hit_pos;
513 }
514 
515 //-----------------------------------------------------------------------------
519 //-----------------------------------------------------------------------------
521  int& Q_f_X,
522  int& Q_l_X,
523  int& Q_f_Y,
524  int& Q_l_Y
525  ) const {
526  ClusterParamGeneric& theClusterParam = static_cast<ClusterParamGeneric&>(theClusterParamBase);
527 
528  // Initialize return variables.
529  Q_f_X = Q_l_X = 0.0;
530  Q_f_Y = Q_l_Y = 0.0;
531 
532  // Obtain boundaries in index units
533  int xmin = theClusterParam.theCluster->minPixelRow();
534  int xmax = theClusterParam.theCluster->maxPixelRow();
535  int ymin = theClusterParam.theCluster->minPixelCol();
536  int ymax = theClusterParam.theCluster->maxPixelCol();
537 
538  // Iterate over the pixels.
539  int isize = theClusterParam.theCluster->size();
540  for (int i = 0; i != isize; ++i) {
541  auto const& pixel = theClusterParam.theCluster->pixel(i);
542  // ggiurgiu@fnal.gov: add pixel charge truncation
543  int pix_adc = pixel.adc;
545  pix_adc = std::min(pix_adc, theClusterParam.pixmx);
546 
547  //
548  // X projection
549  if (pixel.x == xmin)
550  Q_f_X += pix_adc;
551  if (pixel.x == xmax)
552  Q_l_X += pix_adc;
553  //
554  // Y projection
555  if (pixel.y == ymin)
556  Q_f_Y += pix_adc;
557  if (pixel.y == ymax)
558  Q_l_Y += pix_adc;
559  }
560 
561  return;
562 }
563 
564 //============== INFLATED ERROR AND ERRORS FROM DB BELOW ================
565 
566 //-------------------------------------------------------------------------
567 // Hit error in the local frame
568 //-------------------------------------------------------------------------
569 LocalError PixelCPEGeneric::localError(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
570  ClusterParamGeneric& theClusterParam = static_cast<ClusterParamGeneric&>(theClusterParamBase);
571 
572  const bool localPrint = false;
573  // Default errors are the maximum error used for edge clusters.
574  // These are determined by looking at residuals for edge clusters
575  float xerr = EdgeClusterErrorX_ * micronsToCm;
576  float yerr = EdgeClusterErrorY_ * micronsToCm;
577 
578  // Find if cluster is at the module edge.
579  int maxPixelCol = theClusterParam.theCluster->maxPixelCol();
580  int maxPixelRow = theClusterParam.theCluster->maxPixelRow();
581  int minPixelCol = theClusterParam.theCluster->minPixelCol();
582  int minPixelRow = theClusterParam.theCluster->minPixelRow();
583 
584  bool edgex = (theDetParam.theRecTopol->isItEdgePixelInX(minPixelRow)) ||
585  (theDetParam.theRecTopol->isItEdgePixelInX(maxPixelRow));
586  bool edgey = (theDetParam.theRecTopol->isItEdgePixelInY(minPixelCol)) ||
587  (theDetParam.theRecTopol->isItEdgePixelInY(maxPixelCol));
588 
589  unsigned int sizex = theClusterParam.theCluster->sizeX();
590  unsigned int sizey = theClusterParam.theCluster->sizeY();
591  if (MYDEBUG) {
592  if (int(sizex) != (maxPixelRow - minPixelRow + 1))
593  cout << " wrong x" << endl;
594  if (int(sizey) != (maxPixelCol - minPixelCol + 1))
595  cout << " wrong y" << endl;
596  }
597 
598  // Find if cluster contains double (big) pixels.
599  bool bigInX = theDetParam.theRecTopol->containsBigPixelInX(minPixelRow, maxPixelRow);
600  bool bigInY = theDetParam.theRecTopol->containsBigPixelInY(minPixelCol, maxPixelCol);
601 
602  if (localPrint) {
603  cout << " edge clus " << xerr << " " << yerr << endl; //dk
604  if (bigInX || bigInY)
605  cout << " big " << bigInX << " " << bigInY << endl;
606  if (edgex || edgey)
607  cout << " edge " << edgex << " " << edgey << endl;
608  cout << " before if " << UseErrorsFromTemplates_ << " " << theClusterParam.qBin_ << endl;
609  if (theClusterParam.qBin_ == 0)
610  cout << " qbin 0! " << edgex << " " << edgey << " " << bigInX << " " << bigInY << " " << sizex << " " << sizey
611  << endl;
612  }
613 
614  if
616  //
617  // Use template errors
618  //cout << "Track angles are known. We can use either errors from templates or the error parameterization from DB." << endl;
619 
620  if (!edgex) { // Only use this for non-edge clusters
621  if (sizex == 1) {
622  if (!bigInX) {
623  xerr = theClusterParam.sx1;
624  } else {
625  xerr = theClusterParam.sx2;
626  }
627  } else {
628  xerr = theClusterParam.sigmax;
629  }
630  }
631 
632  if (!edgey) { // Only use for non-edge clusters
633  if (sizey == 1) {
634  if (!bigInY) {
635  yerr = theClusterParam.sy1;
636  } else {
637  yerr = theClusterParam.sy2;
638  }
639  } else {
640  yerr = theClusterParam.sigmay;
641  }
642  }
643 
644  if (localPrint) {
645  cout << " in if " << edgex << " " << edgey << " " << sizex << " " << sizey << endl;
646  cout << " errors " << xerr << " " << yerr << " " << theClusterParam.sx1 << " " << theClusterParam.sx2 << " "
647  << theClusterParam.sigmax << endl; //dk
648  }
649  }
650  else { // simple errors
651 
652  // This are the simple errors, hardcoded in the code
653  //cout << "Track angles are not known " << endl;
654  //cout << "Default angle estimation which assumes track from PV (0,0,0) does not work." << endl;
655 
656  if (GeomDetEnumerators::isTrackerPixel(theDetParam.thePart)) {
657  if (GeomDetEnumerators::isBarrel(theDetParam.thePart)) {
658  DetId id = (theDetParam.theDet->geographicalId());
659  int layer = ttopo_.layer(id);
660  if (layer == 1) {
661  if (!edgex) {
662  if (sizex <= xerr_barrel_l1_.size())
663  xerr = xerr_barrel_l1_[sizex - 1];
664  else
665  xerr = xerr_barrel_l1_def_;
666  }
667 
668  if (!edgey) {
669  if (sizey <= yerr_barrel_l1_.size())
670  yerr = yerr_barrel_l1_[sizey - 1];
671  else
672  yerr = yerr_barrel_l1_def_;
673  }
674  } else { // layer 2,3
675  if (!edgex) {
676  if (sizex <= xerr_barrel_ln_.size())
677  xerr = xerr_barrel_ln_[sizex - 1];
678  else
679  xerr = xerr_barrel_ln_def_;
680  }
681 
682  if (!edgey) {
683  if (sizey <= yerr_barrel_ln_.size())
684  yerr = yerr_barrel_ln_[sizey - 1];
685  else
686  yerr = yerr_barrel_ln_def_;
687  }
688  }
689 
690  } else { // EndCap
691 
692  if (!edgex) {
693  if (sizex <= xerr_endcap_.size())
694  xerr = xerr_endcap_[sizex - 1];
695  else
696  xerr = xerr_endcap_def_;
697  }
698 
699  if (!edgey) {
700  if (sizey <= yerr_endcap_.size())
701  yerr = yerr_endcap_[sizey - 1];
702  else
703  yerr = yerr_endcap_def_;
704  }
705  } // end endcap
706  }
707 
708  if (inflate_errors) {
709  int n_bigx = 0;
710  int n_bigy = 0;
711 
712  for (int irow = 0; irow < 7; ++irow) {
713  if (theDetParam.theRecTopol->isItBigPixelInX(irow + minPixelRow))
714  ++n_bigx;
715  }
716 
717  for (int icol = 0; icol < 21; ++icol) {
718  if (theDetParam.theRecTopol->isItBigPixelInY(icol + minPixelCol))
719  ++n_bigy;
720  }
721 
722  xerr = (float)(sizex + n_bigx) * theDetParam.thePitchX / std::sqrt(12.0f);
723  yerr = (float)(sizey + n_bigy) * theDetParam.thePitchY / std::sqrt(12.0f);
724 
725  } // if(inflate_errors)
726 
727  } // end
728 
729 #ifdef EDM_ML_DEBUG
730  if (!(xerr > 0.0))
731  throw cms::Exception("PixelCPEGeneric::localError") << "\nERROR: Negative pixel error xerr = " << xerr << "\n\n";
732 
733  if (!(yerr > 0.0))
734  throw cms::Exception("PixelCPEGeneric::localError") << "\nERROR: Negative pixel error yerr = " << yerr << "\n\n";
735 #endif
736 
737  //if(localPrint) {
738  //cout<<" errors "<<xerr<<" "<<yerr<<endl; //dk
739  //if(theClusterParam.qBin_ == 0) cout<<" qbin 0 "<<xerr<<" "<<yerr<<endl;
740  //}
741 
742  auto xerr_sq = xerr * xerr;
743  auto yerr_sq = yerr * yerr;
744 
745  return LocalError(xerr_sq, 0, yerr_sq);
746 }
747 
749  desc.add<double>("eff_charge_cut_highX", 1.0);
750  desc.add<double>("eff_charge_cut_highY", 1.0);
751  desc.add<double>("eff_charge_cut_lowX", 0.0);
752  desc.add<double>("eff_charge_cut_lowY", 0.0);
753  desc.add<double>("size_cutX", 3.0);
754  desc.add<double>("size_cutY", 3.0);
755  desc.add<double>("EdgeClusterErrorX", 50.0);
756  desc.add<double>("EdgeClusterErrorY", 85.0);
757  desc.add<bool>("inflate_errors", false);
758  desc.add<bool>("inflate_all_errors_no_trk_angle", false);
759  desc.add<bool>("UseErrorsFromTemplates", true);
760  desc.add<bool>("TruncatePixelCharge", true);
761  desc.add<bool>("IrradiationBiasCorrection", false);
762  desc.add<bool>("DoCosmics", false);
763  desc.add<bool>("Upgrade", false);
764  desc.add<bool>("SmallPitch", false);
765 }
#define LogDebug(id)
size
Write out results.
T getParameter(std::string const &) const
std::vector< float > xerr_barrel_l1_
int minPixelCol() const
T y() const
Definition: PV2DBase.h:44
static bool pushfile(int filenum, std::vector< SiPixelGenErrorStore > &pixelTemp, std::string dir="")
float the_eff_charge_cut_highY
std::vector< float > xerr_barrel_ln_
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
const SiPixelCluster * theCluster
Definition: PixelCPEBase.h:85
bool isBarrel(GeomDetEnumerators::SubDetector m)
#define nullptr
T y() const
Definition: PV3DBase.h:60
#define LIKELY(x)
Definition: Likely.h:20
std::vector< SiPixelGenErrorStore > thePixelGenError_
int charge() const
const PixelGeomDetUnit * theDet
Definition: PixelCPEBase.h:58
GeomDetType::SubDetector thePart
Definition: PixelCPEBase.h:63
LocalError localError(DetParam const &theDetParam, ClusterParam &theClusterParam) const override
bool IrradiationBiasCorrection_
const RectangularPixelTopology * theRecTopol
Definition: PixelCPEBase.h:61
std::vector< float > yerr_endcap_
int maxPixelRow() const
PixelCPEGeneric(edm::ParameterSet const &conf, const MagneticField *, const TrackerGeometry &, const TrackerTopology &, const SiPixelLorentzAngle *, const SiPixelGenErrorDBObject *, const SiPixelLorentzAngle *)
The constructor.
int qbin(int id, float cotalpha, float cotbeta, float locBz, float locBx, float qclus, bool irradiationCorrections, int &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)
bool isItEdgePixelInX(int ixbin) const override
const PixelTopology * theTopol
Definition: PixelCPEBase.h:60
int minPixelRow() const
bool LoadTemplatesFromDB_
Definition: PixelCPEBase.h:254
float the_eff_charge_cut_lowX
T sqrt(T t)
Definition: SSEVec.h:19
bool inflate_all_errors_no_trk_angle
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
DetId geographicalId() const
The label of this GeomDet.
Definition: GeomDet.h:64
double f[11][100]
T min(T a, T b)
Definition: MathUtil.h:58
float lorxwidth()
signed lorentz x-width (microns)
ParameterDescriptionBase * add(U const &iLabel, T const &value)
std::vector< float > xerr_endcap_
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
const TrackerTopology & ttopo_
Definition: PixelCPEBase.h:242
const SiPixelGenErrorDBObject * genErrorDBObject_
Definition: PixelCPEBase.h:247
float the_eff_charge_cut_lowY
ClusterParam * createClusterParam(const SiPixelCluster &cl) const override
int maxPixelCol() const
std::vector< float > yerr_barrel_l1_
virtual LocalPoint localPosition(const MeasurementPoint &) const =0
unsigned int layer(const DetId &id) const
std::vector< float > yerr_barrel_ln_
int sizeY() const
Pixel cluster – collection of neighboring pixels above threshold.
Pixel pixel(int i) const
float generic_position_formula(int size, int Q_f, int Q_l, float upper_edge_first_pix, float lower_edge_last_pix, float lorentz_shift, float theThickness, float cot_angle, float pitch, bool first_is_big, bool last_is_big, float eff_charge_cut_low, float eff_charge_cut_high, float size_cut) const
float lorywidth()
signed lorentz y-width (microns)
float y() const
bool containsBigPixelInY(int iymin, int iymax) const override
LocalPoint localPosition(DetParam const &theDetParam, ClusterParam &theClusterParam) const override
bool containsBigPixelInX(int ixmin, int ixmax) const override
T x() const
Definition: PV2DBase.h:43
bool isItEdgePixelInY(int iybin) const override
int sizeX() const
T x() const
Definition: PV3DBase.h:59
static void fillPSetDescription(edm::ParameterSetDescription &desc)
float x() const
void collect_edge_charges(ClusterParam &theClusterParam, int &Q_f_X, int &Q_l_X, int &Q_f_Y, int &Q_l_Y) const
bool isTrackerPixel(GeomDetEnumerators::SubDetector m)
#define constexpr
float the_eff_charge_cut_highX
int size() const