CMS 3D CMS Logo

Classes | Functions
SiPixelTemplateReco Namespace Reference

Classes

struct  ClusMatrix
 

Functions

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)
 
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, float &probQ)
 
int PixelTempReco1D (int id, float cotalpha, float cotbeta, ClusMatrix &cluster, SiPixelTemplate &templ, float &yrec, float &sigmay, float &proby, float &xrec, float &sigmax, float &probx, int &qbin, int speed, float &probQ)
 
int PixelTempReco1D (int id, float cotalpha, float cotbeta, ClusMatrix &cluster, SiPixelTemplate &templ, float &yrec, float &sigmay, float &proby, float &xrec, float &sigmax, float &probx, int &qbin, int speed)
 

Function Documentation

◆ PixelTempReco1D() [1/4]

int SiPixelTemplateReco::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 
)

Reconstruct the best estimate of the hit position for pixel clusters.

Parameters
id- (input) identifier of the template to use
cotalpha- (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
locBz- (input) the sign of this quantity is used to determine whether to flip cot(beta)<0 quantities from cot(beta)>0 (FPix only) for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0 for Phase 1 FPix IP-related tracks, see next comment
locBx- (input) the sign of this quantity is used to determine whether to flip cot(alpha/beta)<0 quantities from cot(alpha/beta)>0 (FPix only) for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0 for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0//!
cotbeta- (input) the cotangent of the beta track angle (see CMS IN 2004/014)
cluster- (input) boost multi_array container of 7x21 array of pixel signals, origin of local coords (0,0) at center of pixel cluster[0][0]. Set dead pixels to small non-zero values (0.1 e).
ydouble- (input) STL vector of 21 element array to flag a double-pixel
xdouble- (input) STL vector of 7 element array to flag a double-pixel
templ- (input) the template used in the reconstruction
yrec- (output) best estimate of y-coordinate of hit in microns
sigmay- (output) best estimate of uncertainty on yrec in microns
proby- (output) probability describing goodness-of-fit for y-reco
xrec- (output) best estimate of x-coordinate of hit in microns
sigmax- (output) best estimate of uncertainty on xrec in microns
probx- (output) probability describing goodness-of-fit for x-reco
qbin- (output) index (0-4) describing the charge of the cluster qbin = 0 Q/Q_avg > 1.5 [few % of all hits] 1 1.5 > Q/Q_avg > 1.0 [~30% of all hits] 2 1.0 > Q/Q_avg > 0.85 [~30% of all hits] 3 0.85 > Q/Q_avg > min1 [~30% of all hits] 4 min1 > Q/Q_avg > min2 [~0.1% of all hits] 5 min2 > Q/Q_avg [~0.1% of all hits]
speed- (input) switch (-2->5) trading speed vs robustness -2 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability w/ VVIObj (better but slower) -1 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability w/ TMath::VavilovI (poorer but faster) 0 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4) 1 faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density 2 faster yet, searches same range as 1 but at 1/2 density 3 fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster) 4 fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster), calculates Q probability w/ VVIObj (better but slower) 5 fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster), calculates Q probability w/ TMath::VavilovI (poorer but faster)
deadpix- (input) bool to indicate that there are dead pixels to be included in the analysis
zeropix- (input) vector of index pairs pointing to the dead pixels
probQ- (output) the Vavilov-distribution-based cluster charge probability
nypix- (output) the projected y-size of the cluster
nxpix- (output) the projected x-size of the cluster

Definition at line 135 of file SiPixelTemplateReco.cc.

References cms::cuda::assert(), BHX, BHY, BXM1, BXM2, BXSIZE, BYM1, BYM2, BYM3, BYSIZE, SiPixelTemplate::chi2xavg(), SiPixelTemplate::chi2xavgone(), SiPixelTemplate::chi2xmin(), SiPixelTemplate::chi2xminone(), SiPixelTemplate::chi2yavg(), SiPixelTemplate::chi2yavgone(), SiPixelTemplate::chi2ymin(), SiPixelTemplate::chi2yminone(), SiPixelTemplate::cxtemp(), SiPixelTemplate::cytemp(), dumpMFGeometry_cfg::delta, SiPixelTemplate::dxone(), SiPixelTemplate::dxtwo(), SiPixelTemplate::dyone(), SiPixelTemplate::dytwo(), ENDL, Exception, f, SiPixelTemplate::fbin(), VVIObjF::fcn(), cms::cuda::for(), mps_fire::i, createfilelist::int, SiPixelTemplate::interpolate(), dqmiolumiharvest::j, isotrackApplyRegressor::k, hltHgcalLayerClustersEE_cfi::kappa, LOGDEBUG, LOGERROR, SiPixelTemplateReco::ClusMatrix::matrix, maxpix, SiPixelTemplateReco::ClusMatrix::mcol, isotrackTrainRegressor::meany, SiPixelTemplateReco::ClusMatrix::mrow, SiPixelTemplate::pixmax(), SiPixelTemplate::qavg(), SiPixelTemplate::qmin(), SiPixelTemplate::qscale(), SiPixelTemplate::s50(), jetUpdater_cfi::sort, HLTSiStripMonitoring_cff::speed, SiPixelTemplate::sxmax(), SiPixelTemplate::sxone(), SiPixelTemplate::sxtwo(), SiPixelTemplate::symax(), SiPixelTemplate::syone(), SiPixelTemplate::sytwo(), theVerboseLevel, TXSIZE, TYSIZE, SiPixelTemplate::vavilov_pars(), SiPixelTemplate::xavg(), SiPixelTemplateReco::ClusMatrix::xdouble, SiPixelTemplate::xflcorr(), SiPixelTemplate::xrms(), SiPixelTemplate::xsigma2(), xsize, SiPixelTemplate::xsize(), SiPixelTemplate::xtemp(), SiPixelTemplate::yavg(), SiPixelTemplateReco::ClusMatrix::ydouble, SiPixelTemplate::yflcorr(), SiPixelTemplate::yrms(), SiPixelTemplate::ysigma2(), ysize, SiPixelTemplate::ysize(), and SiPixelTemplate::ytemp().

Referenced by PixelCPEClusterRepair::callTempReco1D(), PixelCPETemplateReco::localPosition(), and PixelTempReco1D().

156 {
157  // Local variables
158  int i, j, k, minbin, binl, binh, binq, midpix, fypix, lypix, logypx;
159  int fxpix, lxpix, logxpx, shifty, shiftx, nyzero[TYSIZE];
160  int nclusx, nclusy;
161  int deltaj, jmin, jmax, fxbin, lxbin, fybin, lybin, djy, djx;
162  //int fypix2D, lypix2D, fxpix2D, lxpix2D;
163  float sythr, sxthr, rnorm, delta, sigma, sigavg, pseudopix, qscale, q50;
164  float ss2, ssa, sa2, ssba, saba, sba2, rat, fq, qtotal, qpixel, fbin[3];
165  float originx, originy, qfy, qly, qfx, qlx, bias, maxpix, minmax;
166  double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
167  double hchi2, hndof, prvav, mpv, sigmaQ, kappa, xvav, beta2;
168  float ytemp[41][BYSIZE], xtemp[41][BXSIZE], ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
169  float chi2ybin[41], chi2xbin[41], ysig2[BYSIZE], xsig2[BXSIZE];
170  float yw2[BYSIZE], xw2[BXSIZE], ysw[BYSIZE], xsw[BXSIZE];
171  bool yd[BYSIZE], xd[BXSIZE], anyyd, anyxd, calc_probQ, use_VVIObj;
172  float ysize, xsize;
173  const float probmin = {1.110223e-16};
174  const float probQmin = {1.e-5};
175 
176  // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
177 
178  const double mean1pix = {0.100};
179 #ifdef SI_PIXEL_TEMPLATE_STANDALONE
180  const double chi21min = {0.};
181 #else
182  const double chi21min = {0.160};
183 #endif
184 
185  // First, interpolate the template needed to analyze this cluster
186  // check to see of the track direction is in the physical range of the loaded template
187 
188  if (id >= 0) { //if id < 0 bypass interpolation (used in calibration)
189  if (!templ.interpolate(id, cotalpha, cotbeta, locBz, locBx)) {
190  if (theVerboseLevel > 2) {
191  LOGDEBUG("SiPixelTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha
192  << ", cot(beta) = " << cotbeta << ", local B_z = " << locBz
193  << ", local B_x = " << locBx << ", template ID = " << id
194  << ", no reconstruction performed" << ENDL;
195  }
196  return 20;
197  }
198  }
199 
200  // Check to see if Q probability is selected
201 
202  calc_probQ = false;
203  use_VVIObj = false;
204  if (speed < 0) {
205  calc_probQ = true;
206  if (speed < -1)
207  use_VVIObj = true;
208  speed = 0;
209  }
210 
211  if (speed > 3) {
212  calc_probQ = true;
213  if (speed < 5)
214  use_VVIObj = true;
215  speed = 3;
216  }
217 
218  // Get pixel dimensions from the template (to allow multiple detectors in the future)
219 
220  xsize = templ.xsize();
221  ysize = templ.ysize();
222 
223  // Allow Qbin Q/Q_avg fractions to vary to optimize error estimation
224 
225  for (i = 0; i < 3; ++i) {
226  fbin[i] = templ.fbin(i);
227  }
228 
229  // Define size of pseudopixel
230 
231  q50 = templ.s50();
232  pseudopix = 0.2f * q50;
233 
234  // Get charge scaling factor
235 
236  qscale = templ.qscale();
237 
238  // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
239 
240  nclusx = cluster.mrow;
241  nclusy = (int)cluster.mcol;
242  auto const xdouble = cluster.xdouble;
243  auto const ydouble = cluster.ydouble;
244 
245  // First, rescale all pixel charges and compute total charge
246  qtotal = 0.f;
247  for (i = 0; i < nclusx * nclusy; ++i) {
248  cluster.matrix[i] *= qscale;
249  qtotal += cluster.matrix[i];
250  }
251  // Next, sum the total charge and "decapitate" big pixels
252  minmax = templ.pixmax();
253  for (j = 0; j < nclusx; ++j)
254  for (i = 0; i < nclusy; ++i) {
255  maxpix = minmax;
256  if (ydouble[i]) {
257  maxpix *= 2.f;
258  }
259  if (cluster(j, i) > maxpix) {
260  cluster(j, i) = maxpix;
261  }
262  }
263 
264  // Do the cluster repair here
265 
266  if (deadpix) {
267  fypix = BYM3;
268  lypix = -1;
269  memset(nyzero, 0, TYSIZE * sizeof(int));
270  memset(ysum, 0, BYSIZE * sizeof(float));
271  for (i = 0; i < nclusy; ++i) {
272  // Do preliminary cluster projection in y
273  for (j = 0; j < nclusx; ++j) {
274  ysum[i] += cluster(j, i);
275  }
276  if (ysum[i] > 0.f) {
277  // identify ends of cluster to determine what the missing charge should be
278  if (i < fypix) {
279  fypix = i;
280  }
281  if (i > lypix) {
282  lypix = i;
283  }
284  }
285  }
286 
287  // Now loop over dead pixel list and "fix" everything
288 
289  //First see if the cluster ends are redefined and that we have only one dead pixel per column
290 
291  std::vector<std::pair<int, int> >::const_iterator zeroIter = zeropix.begin(), zeroEnd = zeropix.end();
292  for (; zeroIter != zeroEnd; ++zeroIter) {
293  i = zeroIter->second;
294  if (i < 0 || i > TYSIZE - 1) {
295  LOGERROR("SiPixelTemplateReco") << "dead pixel column y-index " << i << ", no reconstruction performed" << ENDL;
296  return 11;
297  }
298 
299  // count the number of dead pixels in each column
300  ++nyzero[i];
301  // allow them to redefine the cluster ends
302  if (i < fypix) {
303  fypix = i;
304  }
305  if (i > lypix) {
306  lypix = i;
307  }
308  }
309 
310  nypix = lypix - fypix + 1;
311 
312  // Now adjust the charge in the dead pixels to sum to 0.5*truncation value in the end columns and the truncation value in the interior columns
313 
314  for (zeroIter = zeropix.begin(); zeroIter != zeroEnd; ++zeroIter) {
315  i = zeroIter->second;
316  j = zeroIter->first;
317  if (j < 0 || j > TXSIZE - 1) {
318  LOGERROR("SiPixelTemplateReco") << "dead pixel column x-index " << j << ", no reconstruction performed" << ENDL;
319  return 12;
320  }
321  if ((i == fypix || i == lypix) && nypix > 1) {
322  maxpix = templ.symax() / 2.;
323  } else {
324  maxpix = templ.symax();
325  }
326  if (ydouble[i]) {
327  maxpix *= 2.;
328  }
329  if (nyzero[i] > 0 && nyzero[i] < 3) {
330  qpixel = (maxpix - ysum[i]) / (float)nyzero[i];
331  } else {
332  qpixel = 1.;
333  }
334  if (qpixel < 1.) {
335  qpixel = 1.;
336  }
337  cluster(j, i) = qpixel;
338  // Adjust the total cluster charge to reflect the charge of the "repaired" cluster
339  qtotal += qpixel;
340  }
341  // End of cluster repair section
342  }
343 
344  // Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container
345 
346  for (i = 0; i < BYSIZE; ++i) {
347  ysum[i] = 0.f;
348  yd[i] = false;
349  }
350  k = 0;
351  anyyd = false;
352  for (i = 0; i < nclusy; ++i) {
353  for (j = 0; j < nclusx; ++j) {
354  ysum[k] += cluster(j, i);
355  }
356 
357  // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
358 
359  if (ydouble[i]) {
360  ysum[k] /= 2.f;
361  ysum[k + 1] = ysum[k];
362  yd[k] = true;
363  yd[k + 1] = false;
364  k = k + 2;
365  anyyd = true;
366  } else {
367  yd[k] = false;
368  ++k;
369  }
370  if (k > BYM1) {
371  break;
372  }
373  }
374 
375  // Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container
376 
377  for (i = 0; i < BXSIZE; ++i) {
378  xsum[i] = 0.f;
379  xd[i] = false;
380  }
381  k = 0;
382  anyxd = false;
383  for (j = 0; j < nclusx; ++j) {
384  for (i = 0; i < nclusy; ++i) {
385  xsum[k] += cluster(j, i);
386  }
387 
388  // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
389 
390  if (xdouble[j]) {
391  xsum[k] /= 2.;
392  xsum[k + 1] = xsum[k];
393  xd[k] = true;
394  xd[k + 1] = false;
395  k = k + 2;
396  anyxd = true;
397  } else {
398  xd[k] = false;
399  ++k;
400  }
401  if (k > BXM1) {
402  break;
403  }
404  }
405 
406  // next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx
407 
408  fypix = -1;
409  nypix = 0;
410  lypix = 0;
411  logypx = 0;
412  for (i = 0; i < BYSIZE; ++i) {
413  if (ysum[i] > 0.f) {
414  if (fypix == -1) {
415  fypix = i;
416  }
417  if (!yd[i]) {
418  ysort[logypx] = ysum[i];
419  ++logypx;
420  }
421  ++nypix;
422  lypix = i;
423  }
424  }
425 
426  // dlengthy = (float)nypix - templ.clsleny();
427 
428  // Make sure cluster is continuous
429 
430  if ((lypix - fypix + 1) != nypix || nypix == 0) {
431  LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster doesn't agree with number of pixels above threshold"
432  << ENDL;
433  if (theVerboseLevel > 2) {
434  LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
435  for (i = 0; i < BYSIZE - 1; ++i) {
436  LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";
437  }
438  LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE - 1] << ENDL;
439  }
440 
441  return 1;
442  }
443 
444  // If cluster is longer than max template size, technique fails
445 
446  if (nypix > TYSIZE) {
447  LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster is larger than maximum template size" << ENDL;
448  if (theVerboseLevel > 2) {
449  LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
450  for (i = 0; i < BYSIZE - 1; ++i) {
451  LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";
452  }
453  LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE - 1] << ENDL;
454  }
455 
456  return 6;
457  }
458 
459  // Remember these numbers for later
460 
461  //fypix2D = fypix;
462  //lypix2D = lypix;
463 
464  // next, center the cluster on template center if necessary
465 
466  midpix = (fypix + lypix) / 2;
467  shifty = templ.cytemp() - midpix;
468 
469  // calculate new cluster boundaries
470 
471  int lytmp = lypix + shifty;
472  int fytmp = fypix + shifty;
473 
474  // Check the boundaries
475 
476  if (fytmp <= 1) {
477  return 8;
478  }
479  if (lytmp >= BYM2) {
480  return 8;
481  }
482 
483  // If OK, shift everything
484 
485  if (shifty > 0) {
486  for (i = lypix; i >= fypix; --i) {
487  ysum[i + shifty] = ysum[i];
488  ysum[i] = 0.;
489  yd[i + shifty] = yd[i];
490  yd[i] = false;
491  }
492  } else if (shifty < 0) {
493  for (i = fypix; i <= lypix; ++i) {
494  ysum[i + shifty] = ysum[i];
495  ysum[i] = 0.;
496  yd[i + shifty] = yd[i];
497  yd[i] = false;
498  }
499  }
500 
501  lypix = lytmp;
502  fypix = fytmp;
503 
504  // add pseudopixels
505 
506  ysum[fypix - 1] = pseudopix;
507  ysum[fypix - 2] = pseudopix;
508  ysum[lypix + 1] = pseudopix;
509  ysum[lypix + 2] = pseudopix;
510 
511  // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
512 
513  if (ydouble[0]) {
514  originy = -0.5f;
515  } else {
516  originy = 0.f;
517  }
518 
519  // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx
520 
521  fxpix = -1;
522  nxpix = 0;
523  lxpix = 0;
524  logxpx = 0;
525  for (i = 0; i < BXSIZE; ++i) {
526  if (xsum[i] > 0.) {
527  if (fxpix == -1) {
528  fxpix = i;
529  }
530  if (!xd[i]) {
531  xsort[logxpx] = xsum[i];
532  ++logxpx;
533  }
534  ++nxpix;
535  lxpix = i;
536  }
537  }
538 
539  // dlengthx = (float)nxpix - templ.clslenx();
540 
541  // Make sure cluster is continuous
542 
543  if ((lxpix - fxpix + 1) != nxpix) {
544  LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold"
545  << ENDL;
546  if (theVerboseLevel > 2) {
547  LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
548  for (i = 0; i < BXSIZE - 1; ++i) {
549  LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";
550  }
551  LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE - 1] << ENDL;
552  }
553 
554  return 2;
555  }
556 
557  // If cluster is longer than max template size, technique fails
558 
559  if (nxpix > TXSIZE) {
560  LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
561  if (theVerboseLevel > 2) {
562  LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
563  for (i = 0; i < BXSIZE - 1; ++i) {
564  LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";
565  }
566  LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE - 1] << ENDL;
567  }
568 
569  return 7;
570  }
571 
572  // Remember these numbers for later
573 
574  //fxpix2D = fxpix;
575  //lxpix2D = lxpix;
576 
577  // next, center the cluster on template center if necessary
578 
579  midpix = (fxpix + lxpix) / 2;
580  shiftx = templ.cxtemp() - midpix;
581 
582  // calculate new cluster boundaries
583 
584  int lxtmp = lxpix + shiftx;
585  int fxtmp = fxpix + shiftx;
586 
587  // Check the boundaries
588 
589  if (fxtmp <= 1) {
590  return 9;
591  }
592  if (lxtmp >= BXM2) {
593  return 9;
594  }
595 
596  // If OK, shift everything
597 
598  if (shiftx > 0) {
599  for (i = lxpix; i >= fxpix; --i) {
600  xsum[i + shiftx] = xsum[i];
601  xsum[i] = 0.;
602  xd[i + shiftx] = xd[i];
603  xd[i] = false;
604  }
605  } else if (shiftx < 0) {
606  for (i = fxpix; i <= lxpix; ++i) {
607  xsum[i + shiftx] = xsum[i];
608  xsum[i] = 0.;
609  xd[i + shiftx] = xd[i];
610  xd[i] = false;
611  }
612  }
613 
614  lxpix = lxtmp;
615  fxpix = fxtmp;
616 
617  xsum[fxpix - 1] = pseudopix;
618  xsum[fxpix - 2] = pseudopix;
619  xsum[lxpix + 1] = pseudopix;
620  xsum[lxpix + 2] = pseudopix;
621 
622  // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
623 
624  if (xdouble[0]) {
625  originx = -0.5f;
626  } else {
627  originx = 0.f;
628  }
629 
630  // uncertainty and final corrections depend upon total charge bin
631 
632  fq = qtotal / templ.qavg();
633  if (fq > fbin[0]) {
634  binq = 0;
635  } else {
636  if (fq > fbin[1]) {
637  binq = 1;
638  } else {
639  if (fq > fbin[2]) {
640  binq = 2;
641  } else {
642  binq = 3;
643  }
644  }
645  }
646 
647  // Return the charge bin via the parameter list unless the charge is too small (then flag it)
648 
649  qbin = binq;
650  if (!deadpix && qtotal < 0.95f * templ.qmin()) {
651  qbin = 5;
652  } else {
653  if (!deadpix && qtotal < 0.95f * templ.qmin(1)) {
654  qbin = 4;
655  }
656  }
657  if (theVerboseLevel > 9) {
658  LOGDEBUG("SiPixelTemplateReco") << "ID = " << id << " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta
659  << " nclusx = " << nclusx << " nclusy = " << nclusy << ENDL;
660  }
661 
662  // Next, copy the y- and x-templates to local arrays
663 
664  // First, decide on chi^2 min search parameters
665 
666 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
667  if (speed < 0 || speed > 3) {
668  throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::PixelTempReco1D called with illegal speed = " << speed
669  << std::endl;
670  }
671 #else
672  assert(speed >= 0 && speed < 4);
673 #endif
674  fybin = 3;
675  lybin = 37;
676  fxbin = 3;
677  lxbin = 37;
678  djy = 1;
679  djx = 1;
680  if (speed > 0) {
681  fybin = 8;
682  lybin = 32;
683  if (yd[fypix]) {
684  fybin = 4;
685  lybin = 36;
686  }
687  if (lypix > fypix) {
688  if (yd[lypix - 1]) {
689  fybin = 4;
690  lybin = 36;
691  }
692  }
693  fxbin = 8;
694  lxbin = 32;
695  if (xd[fxpix]) {
696  fxbin = 4;
697  lxbin = 36;
698  }
699  if (lxpix > fxpix) {
700  if (xd[lxpix - 1]) {
701  fxbin = 4;
702  lxbin = 36;
703  }
704  }
705  }
706 
707  if (speed > 1) {
708  djy = 2;
709  djx = 2;
710  if (speed > 2) {
711  if (!anyyd) {
712  djy = 4;
713  }
714  if (!anyxd) {
715  djx = 4;
716  }
717  }
718  }
719 
720  if (theVerboseLevel > 9) {
721  LOGDEBUG("SiPixelTemplateReco") << "fypix " << fypix << " lypix = " << lypix << " fybin = " << fybin
722  << " lybin = " << lybin << " djy = " << djy << " logypx = " << logypx << ENDL;
723  LOGDEBUG("SiPixelTemplateReco") << "fxpix " << fxpix << " lxpix = " << lxpix << " fxbin = " << fxbin
724  << " lxbin = " << lxbin << " djx = " << djx << " logxpx = " << logxpx << ENDL;
725  }
726 
727  // Now do the copies
728 
729  templ.ytemp(fybin, lybin, ytemp);
730 
731  templ.xtemp(fxbin, lxbin, xtemp);
732 
733  // Do the y-reconstruction first
734 
735  // Apply the first-pass template algorithm to all clusters
736 
737  // Modify the template if double pixels are present
738 
739  if (nypix > logypx) {
740  i = fypix;
741  while (i < lypix) {
742  if (yd[i] && !yd[i + 1]) {
743  for (j = fybin; j <= lybin; ++j) {
744  // Sum the adjacent cells and put the average signal in both
745 
746  sigavg = (ytemp[j][i] + ytemp[j][i + 1]) / 2.f;
747  ytemp[j][i] = sigavg;
748  ytemp[j][i + 1] = sigavg;
749  }
750  i += 2;
751  } else {
752  ++i;
753  }
754  }
755  }
756 
757  // Define the maximum signal to allow before de-weighting a pixel
758 
759  sythr = 1.1 * (templ.symax());
760 
761  // Make sure that there will be at least two pixels that are not de-weighted
762 
763  std::sort(&ysort[0], &ysort[logypx]);
764  if (logypx == 1) {
765  sythr = 1.01f * ysort[0];
766  } else {
767  if (ysort[1] > sythr) {
768  sythr = 1.01f * ysort[1];
769  }
770  }
771 
772  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
773 
774  // for(i=0; i<BYSIZE; ++i) { ysig2[i] = 0.;}
775  templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
776 
777  // Find the template bin that minimizes the Chi^2
778 
779  chi2ymin = 1.e15;
780  for (i = fybin; i <= lybin; ++i) {
781  chi2ybin[i] = -1.e15f;
782  }
783  ss2 = 0.f;
784  for (i = fypix - 2; i <= lypix + 2; ++i) {
785  yw2[i] = 1.f / ysig2[i];
786  ysw[i] = ysum[i] * yw2[i];
787  ss2 += ysum[i] * ysw[i];
788  }
789 
790  minbin = -1;
791  deltaj = djy;
792  jmin = fybin;
793  jmax = lybin;
794  while (deltaj > 0) {
795  for (j = jmin; j <= jmax; j += deltaj) {
796  if (chi2ybin[j] < -100.f) {
797  ssa = 0.f;
798  sa2 = 0.f;
799  for (i = fypix - 2; i <= lypix + 2; ++i) {
800  ssa += ysw[i] * ytemp[j][i];
801  sa2 += ytemp[j][i] * ytemp[j][i] * yw2[i];
802  }
803  rat = ssa / ss2;
804  if (rat <= 0.f) {
805  LOGERROR("SiPixelTemplateReco") << "illegal chi2ymin normalization (1) = " << rat << ENDL;
806  rat = 1.;
807  }
808  chi2ybin[j] = ss2 - 2.f * ssa / rat + sa2 / (rat * rat);
809  }
810  if (chi2ybin[j] < chi2ymin) {
811  chi2ymin = chi2ybin[j];
812  minbin = j;
813  }
814  }
815  deltaj /= 2;
816  if (minbin > fybin) {
817  jmin = minbin - deltaj;
818  } else {
819  jmin = fybin;
820  }
821  if (minbin < lybin) {
822  jmax = minbin + deltaj;
823  } else {
824  jmax = lybin;
825  }
826  }
827 
828  if (theVerboseLevel > 9) {
829  LOGDEBUG("SiPixelTemplateReco") << "minbin " << minbin << " chi2ymin = " << chi2ymin << ENDL;
830  }
831 
832  // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
833 
834  if (logypx == 1) {
835  if (nypix == 1) {
836  delta = templ.dyone();
837  sigma = templ.syone();
838  } else {
839  delta = templ.dytwo();
840  sigma = templ.sytwo();
841  }
842 
843  yrec = 0.5f * (fypix + lypix - 2 * shifty + 2.f * originy) * ysize - delta;
844  if (sigma <= 0.f) {
845  sigmay = 43.3f;
846  } else {
847  sigmay = sigma;
848  }
849 
850  // Do probability calculation for one-pixel clusters
851 
852  chi21max = fmax(chi21min, (double)templ.chi2yminone());
853  chi2ymin -= chi21max;
854  if (chi2ymin < 0.) {
855  chi2ymin = 0.;
856  }
857  // proby = gsl_cdf_chisq_Q(chi2ymin, mean1pix);
858  meany = fmax(mean1pix, (double)templ.chi2yavgone());
859  hchi2 = chi2ymin / 2.;
860  hndof = meany / 2.;
861  proby = 1. - TMath::Gamma(hndof, hchi2);
862 
863  } else {
864  // For cluster > 1 pix, make the second, interpolating pass with the templates
865 
866  binl = minbin - 1;
867  binh = binl + 2;
868  if (binl < fybin) {
869  binl = fybin;
870  }
871  if (binh > lybin) {
872  binh = lybin;
873  }
874  ssa = 0.;
875  sa2 = 0.;
876  ssba = 0.;
877  saba = 0.;
878  sba2 = 0.;
879  for (i = fypix - 2; i <= lypix + 2; ++i) {
880  ssa += ysw[i] * ytemp[binl][i];
881  sa2 += ytemp[binl][i] * ytemp[binl][i] * yw2[i];
882  ssba += ysw[i] * (ytemp[binh][i] - ytemp[binl][i]);
883  saba += ytemp[binl][i] * (ytemp[binh][i] - ytemp[binl][i]) * yw2[i];
884  sba2 += (ytemp[binh][i] - ytemp[binl][i]) * (ytemp[binh][i] - ytemp[binl][i]) * yw2[i];
885  }
886 
887  // rat is the fraction of the "distance" from template a to template b
888 
889  rat = (ssba * ssa - ss2 * saba) / (ss2 * sba2 - ssba * ssba);
890  if (rat < 0.f) {
891  rat = 0.f;
892  }
893  if (rat > 1.f) {
894  rat = 1.0f;
895  }
896  rnorm = (ssa + rat * ssba) / ss2;
897 
898  // Calculate the charges in the first and last pixels
899 
900  qfy = ysum[fypix];
901  if (yd[fypix]) {
902  qfy += ysum[fypix + 1];
903  }
904  if (logypx > 1) {
905  qly = ysum[lypix];
906  if (yd[lypix - 1]) {
907  qly += ysum[lypix - 1];
908  }
909  } else {
910  qly = qfy;
911  }
912 
913  // Now calculate the mean bias correction and uncertainties
914 
915  float qyfrac = (qfy - qly) / (qfy + qly);
916  bias = templ.yflcorr(binq, qyfrac) + templ.yavg(binq);
917 
918  // uncertainty and final correction depend upon charge bin
919 
920  yrec = (0.125f * binl + BHY - 2.5f + rat * (binh - binl) * 0.125f - (float)shifty + originy) * ysize - bias;
921  sigmay = templ.yrms(binq);
922 
923  // Do goodness of fit test in y
924 
925  if (rnorm <= 0.) {
926  LOGERROR("SiPixelTemplateReco") << "illegal chi2y normalization (2) = " << rnorm << ENDL;
927  rnorm = 1.;
928  }
929  chi2y = ss2 - 2. / rnorm * ssa - 2. / rnorm * rat * ssba +
930  (sa2 + 2. * rat * saba + rat * rat * sba2) / (rnorm * rnorm) - templ.chi2ymin(binq);
931  if (chi2y < 0.0) {
932  chi2y = 0.0;
933  }
934  meany = templ.chi2yavg(binq);
935  if (meany < 0.01) {
936  meany = 0.01;
937  }
938  // gsl function that calculates the chi^2 tail prob for non-integral dof
939  // proby = gsl_cdf_chisq_Q(chi2y, meany);
940  // proby = ROOT::Math::chisquared_cdf_c(chi2y, meany);
941  hchi2 = chi2y / 2.;
942  hndof = meany / 2.;
943  proby = 1. - TMath::Gamma(hndof, hchi2);
944  }
945 
946  // Do the x-reconstruction next
947 
948  // Apply the first-pass template algorithm to all clusters
949 
950  // Modify the template if double pixels are present
951 
952  if (nxpix > logxpx) {
953  i = fxpix;
954  while (i < lxpix) {
955  if (xd[i] && !xd[i + 1]) {
956  for (j = fxbin; j <= lxbin; ++j) {
957  // Sum the adjacent cells and put the average signal in both
958 
959  sigavg = (xtemp[j][i] + xtemp[j][i + 1]) / 2.f;
960  xtemp[j][i] = sigavg;
961  xtemp[j][i + 1] = sigavg;
962  }
963  i += 2;
964  } else {
965  ++i;
966  }
967  }
968  }
969 
970  // Define the maximum signal to allow before de-weighting a pixel
971 
972  sxthr = 1.1f * templ.sxmax();
973 
974  // Make sure that there will be at least two pixels that are not de-weighted
975  std::sort(&xsort[0], &xsort[logxpx]);
976  if (logxpx == 1) {
977  sxthr = 1.01f * xsort[0];
978  } else {
979  if (xsort[1] > sxthr) {
980  sxthr = 1.01f * xsort[1];
981  }
982  }
983 
984  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
985 
986  // for(i=0; i<BXSIZE; ++i) { xsig2[i] = 0.; }
987  templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
988 
989  // Find the template bin that minimizes the Chi^2
990 
991  chi2xmin = 1.e15;
992  for (i = fxbin; i <= lxbin; ++i) {
993  chi2xbin[i] = -1.e15f;
994  }
995  ss2 = 0.f;
996  for (i = fxpix - 2; i <= lxpix + 2; ++i) {
997  xw2[i] = 1.f / xsig2[i];
998  xsw[i] = xsum[i] * xw2[i];
999  ss2 += xsum[i] * xsw[i];
1000  }
1001  minbin = -1;
1002  deltaj = djx;
1003  jmin = fxbin;
1004  jmax = lxbin;
1005  while (deltaj > 0) {
1006  for (j = jmin; j <= jmax; j += deltaj) {
1007  if (chi2xbin[j] < -100.f) {
1008  ssa = 0.f;
1009  sa2 = 0.f;
1010  // std::cout << j << " ";
1011  for (i = fxpix - 2; i <= lxpix + 2; ++i) {
1012  // std::cout << xtemp[j][i] << ", ";
1013  ssa += xsw[i] * xtemp[j][i];
1014  sa2 += xtemp[j][i] * xtemp[j][i] * xw2[i];
1015  }
1016  // std::cout << std::endl;
1017  rat = ssa / ss2;
1018  if (rat <= 0.f) {
1019  LOGERROR("SiPixelTemplateReco") << "illegal chi2xmin normalization (1) = " << rat << ENDL;
1020  rat = 1.;
1021  }
1022  chi2xbin[j] = ss2 - 2.f * ssa / rat + sa2 / (rat * rat);
1023  }
1024  if (chi2xbin[j] < chi2xmin) {
1025  chi2xmin = chi2xbin[j];
1026  minbin = j;
1027  }
1028  }
1029  deltaj /= 2;
1030  if (minbin > fxbin) {
1031  jmin = minbin - deltaj;
1032  } else {
1033  jmin = fxbin;
1034  }
1035  if (minbin < lxbin) {
1036  jmax = minbin + deltaj;
1037  } else {
1038  jmax = lxbin;
1039  }
1040  }
1041 
1042  if (theVerboseLevel > 9) {
1043  LOGDEBUG("SiPixelTemplateReco") << "minbin " << minbin << " chi2xmin = " << chi2xmin << ENDL;
1044  }
1045 
1046  // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
1047 
1048  if (logxpx == 1) {
1049  if (nxpix == 1) {
1050  delta = templ.dxone();
1051  sigma = templ.sxone();
1052  } else {
1053  delta = templ.dxtwo();
1054  sigma = templ.sxtwo();
1055  }
1056  xrec = 0.5 * (fxpix + lxpix - 2 * shiftx + 2. * originx) * xsize - delta;
1057  if (sigma <= 0.) {
1058  sigmax = 28.9;
1059  } else {
1060  sigmax = sigma;
1061  }
1062 
1063  // Do probability calculation for one-pixel clusters
1064 
1065  chi21max = fmax(chi21min, (double)templ.chi2xminone());
1066  chi2xmin -= chi21max;
1067  if (chi2xmin < 0.) {
1068  chi2xmin = 0.;
1069  }
1070  meanx = fmax(mean1pix, (double)templ.chi2xavgone());
1071  hchi2 = chi2xmin / 2.;
1072  hndof = meanx / 2.;
1073  probx = 1. - TMath::Gamma(hndof, hchi2);
1074 
1075  } else {
1076  // Now make the second, interpolating pass with the templates
1077 
1078  binl = minbin - 1;
1079  binh = binl + 2;
1080  if (binl < fxbin) {
1081  binl = fxbin;
1082  }
1083  if (binh > lxbin) {
1084  binh = lxbin;
1085  }
1086  ssa = 0.;
1087  sa2 = 0.;
1088  ssba = 0.;
1089  saba = 0.;
1090  sba2 = 0.;
1091  for (i = fxpix - 2; i <= lxpix + 2; ++i) {
1092  ssa += xsw[i] * xtemp[binl][i];
1093  sa2 += xtemp[binl][i] * xtemp[binl][i] * xw2[i];
1094  ssba += xsw[i] * (xtemp[binh][i] - xtemp[binl][i]);
1095  saba += xtemp[binl][i] * (xtemp[binh][i] - xtemp[binl][i]) * xw2[i];
1096  sba2 += (xtemp[binh][i] - xtemp[binl][i]) * (xtemp[binh][i] - xtemp[binl][i]) * xw2[i];
1097  }
1098 
1099  // rat is the fraction of the "distance" from template a to template b
1100 
1101  rat = (ssba * ssa - ss2 * saba) / (ss2 * sba2 - ssba * ssba);
1102  if (rat < 0.f) {
1103  rat = 0.f;
1104  }
1105  if (rat > 1.f) {
1106  rat = 1.0f;
1107  }
1108  rnorm = (ssa + rat * ssba) / ss2;
1109 
1110  // Calculate the charges in the first and last pixels
1111 
1112  qfx = xsum[fxpix];
1113  if (xd[fxpix]) {
1114  qfx += xsum[fxpix + 1];
1115  }
1116  if (logxpx > 1) {
1117  qlx = xsum[lxpix];
1118  if (xd[lxpix - 1]) {
1119  qlx += xsum[lxpix - 1];
1120  }
1121  } else {
1122  qlx = qfx;
1123  }
1124 
1125  // Now calculate the mean bias correction and uncertainties
1126 
1127  float qxfrac = (qfx - qlx) / (qfx + qlx);
1128  bias = templ.xflcorr(binq, qxfrac) + templ.xavg(binq);
1129 
1130  // uncertainty and final correction depend upon charge bin
1131 
1132  xrec = (0.125f * binl + BHX - 2.5f + rat * (binh - binl) * 0.125f - (float)shiftx + originx) * xsize - bias;
1133  sigmax = templ.xrms(binq);
1134 
1135  // Do goodness of fit test in x
1136 
1137  if (rnorm <= 0.) {
1138  LOGERROR("SiPixelTemplateReco") << "illegal chi2x normalization (2) = " << rnorm << ENDL;
1139  rnorm = 1.;
1140  }
1141  chi2x = ss2 - 2. / rnorm * ssa - 2. / rnorm * rat * ssba +
1142  (sa2 + 2. * rat * saba + rat * rat * sba2) / (rnorm * rnorm) - templ.chi2xmin(binq);
1143  if (chi2x < 0.0) {
1144  chi2x = 0.0;
1145  }
1146  meanx = templ.chi2xavg(binq);
1147  if (meanx < 0.01) {
1148  meanx = 0.01;
1149  }
1150  // gsl function that calculates the chi^2 tail prob for non-integral dof
1151  // probx = gsl_cdf_chisq_Q(chi2x, meanx);
1152  // probx = ROOT::Math::chisquared_cdf_c(chi2x, meanx, trx0);
1153  hchi2 = chi2x / 2.;
1154  hndof = meanx / 2.;
1155  probx = 1. - TMath::Gamma(hndof, hchi2);
1156  }
1157 
1158  // Don't return exact zeros for the probability
1159 
1160  if (proby < probmin) {
1161  proby = probmin;
1162  }
1163  if (probx < probmin) {
1164  probx = probmin;
1165  }
1166 
1167  // Decide whether to generate a cluster charge probability
1168 
1169  if (calc_probQ) {
1170  // Calculate the Vavilov probability that the cluster charge is OK
1171 
1172  templ.vavilov_pars(mpv, sigmaQ, kappa);
1173 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1174  if ((sigmaQ <= 0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
1175  throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::Vavilov parameters mpv/sigmaQ/kappa = " << mpv << "/"
1176  << sigmaQ << "/" << kappa << std::endl;
1177  }
1178 #else
1179  assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
1180 #endif
1181  xvav = ((double)qtotal - mpv) / sigmaQ;
1182  beta2 = 1.;
1183  if (use_VVIObj) {
1184  // VVIObj is a private port of CERNLIB VVIDIS
1185  VVIObjF vvidist(kappa);
1186  prvav = vvidist.fcn(xvav);
1187  } else {
1188  // Use faster but less accurate TMath Vavilov distribution function
1189  prvav = TMath::VavilovI(xvav, kappa, beta2);
1190  }
1191  // Change to upper tail probability
1192  // if(prvav > 0.5) prvav = 1. - prvav;
1193  // probQ = (float)(2.*prvav);
1194  probQ = 1. - prvav;
1195  if (probQ < probQmin) {
1196  probQ = probQmin;
1197  }
1198  } else {
1199  probQ = -1;
1200  }
1201 
1202  return 0;
1203 } // PixelTempReco1D
float fbin(int i)
Return lower bound of Qbin definition.
#define BXSIZE
float chi2xminone()
//!< minimum of x chi^2 for 1 pixel clusters
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
float symax()
average pixel signal for y-projection of cluster
float yavg(int i)
average y-bias of reconstruction binned in 4 charge bins
int cytemp()
Return central pixel of y template pixels above readout threshold.
#define BYSIZE
void xsigma2(int fxpix, int lxpix, float sxthr, float xsum[13+4], float xsig2[13+4])
#define TXSIZE
for(int i=first, nt=offsets[nh];i< nt;i+=gridDim.x *blockDim.x)
float chi2xmin(int i)
minimum y chi^2 in 4 charge bins
#define BXM1
float chi2ymin(int i)
minimum y chi^2 in 4 charge bins
float xrms(int i)
average x-rms of reconstruction binned in 4 charge bins
float qmin()
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
static const int maxpix
const Int_t ysize
float xflcorr(int binq, float qflx)
float sytwo()
rms for one double-pixel y-clusters
assert(be >=bs)
float chi2yminone()
//!< minimum of y chi^2 for 1 pixel clusters
float sxone()
rms for one pixel x-clusters
#define BYM1
#define BXM2
void ytemp(int fybin, int lybin, float ytemplate[41][21+4])
int cxtemp()
Return central pixel of x-template pixels above readout threshold.
float qscale()
charge scaling factor
float dxone()
mean offset/correction for one pixel x-clusters
float chi2yavg(int i)
average y chi^2 in 4 charge bins
void vavilov_pars(double &mpv, double &sigma, double &kappa)
float yrms(int i)
average y-rms of reconstruction binned in 4 charge bins
#define BHX
float yflcorr(int binq, float qfly)
float xsize()
pixel x-size (microns)
void ysigma2(int fypix, int lypix, float sythr, float ysum[21+4], float ysig2[21+4])
#define BYM2
double f[11][100]
float sxtwo()
rms for one double-pixel x-clusters
#define BYM3
float dytwo()
mean offset/correction for one double-pixel y-clusters
#define LOGDEBUG(x)
float s50()
1/2 of the pixel threshold signal in electrons
#define TYSIZE
void xtemp(int fxbin, int lxbin, float xtemplate[41][13+4])
#define LOGERROR(x)
float syone()
rms for one pixel y-clusters
static const int theVerboseLevel
float chi2yavgone()
//!< average y chi^2 for 1 pixel clusters
float qavg()
average cluster charge for this set of track angles
#define BHY
float sxmax()
average pixel signal for x-projection of cluster
float chi2xavgone()
//!< average x chi^2 for 1 pixel clusters
float pixmax()
maximum pixel charge
#define ENDL
float chi2xavg(int i)
averaage x chi^2 in 4 charge bins
const Int_t xsize
float dyone()
mean offset/correction for one pixel y-clusters
float xavg(int i)
average x-bias of reconstruction binned in 4 charge bins
float dxtwo()
mean offset/correction for one double-pixel x-clusters
float ysize()
pixel y-size (microns)

◆ PixelTempReco1D() [2/4]

int SiPixelTemplateReco::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,
float &  probQ 
)

Reconstruct the best estimate of the hit position for pixel clusters.

Parameters
id- (input) identifier of the template to use
cotalpha- (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
cotbeta- (input) the cotangent of the beta track angle (see CMS IN 2004/014)
locBz- (input) the sign of this quantity is used to determine whether to flip cot(beta)<0 quantities from cot(beta)>0 (FPix only) for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0 for Phase 1 FPix IP-related tracks, see next comment
locBx- (input) the sign of this quantity is used to determine whether to flip cot(alpha/beta)<0 quantities from cot(alpha/beta)>0 (FPix only) for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0 for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0//!
cotbeta- (input) the cotangent of the beta track angle (see CMS IN 2004/014)
cluster- (input) boost multi_array container of 7x21 array of pixel signals, origin of local coords (0,0) at center of pixel cluster[0][0].
ydouble- (input) STL vector of 21 element array to flag a double-pixel
xdouble- (input) STL vector of 7 element array to flag a double-pixel
templ- (input) the template used in the reconstruction
yrec- (output) best estimate of y-coordinate of hit in microns
sigmay- (output) best estimate of uncertainty on yrec in microns
proby- (output) probability describing goodness-of-fit for y-reco
xrec- (output) best estimate of x-coordinate of hit in microns
sigmax- (output) best estimate of uncertainty on xrec in microns
probx- (output) probability describing goodness-of-fit for x-reco
qbin- (output) index (0-4) describing the charge of the cluster [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin]
speed- (input) switch (-1-4) trading speed vs robustness -1 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability 0 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4) 1 faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density 2 faster yet, searches same range as 1 but at 1/2 density 3 fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster) 4 fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster), calculates Q probability
probQ- (output) the Vavilov-distribution-based cluster charge probability

Definition at line 1239 of file SiPixelTemplateReco.cc.

References PixelTempReco1D(), and HLTSiStripMonitoring_cff::speed.

1256 {
1257  // Local variables
1258  const bool deadpix = false;
1259  std::vector<std::pair<int, int> > zeropix;
1260  int nypix, nxpix;
1261 
1263  cotalpha,
1264  cotbeta,
1265  locBz,
1266  locBx,
1267  cluster,
1268  templ,
1269  yrec,
1270  sigmay,
1271  proby,
1272  xrec,
1273  sigmax,
1274  probx,
1275  qbin,
1276  speed,
1277  deadpix,
1278  zeropix,
1279  probQ,
1280  nypix,
1281  nxpix);
1282 
1283 } // PixelTempReco1D
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)

◆ PixelTempReco1D() [3/4]

int SiPixelTemplateReco::PixelTempReco1D ( int  id,
float  cotalpha,
float  cotbeta,
ClusMatrix cluster,
SiPixelTemplate templ,
float &  yrec,
float &  sigmay,
float &  proby,
float &  xrec,
float &  sigmax,
float &  probx,
int &  qbin,
int  speed,
float &  probQ 
)

Reconstruct the best estimate of the hit position for pixel clusters.

Parameters
id- (input) identifier of the template to use
cotalpha- (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
cotbeta- (input) the cotangent of the beta track angle (see CMS IN 2004/014)
cluster- (input) boost multi_array container of 7x21 array of pixel signals, origin of local coords (0,0) at center of pixel cluster[0][0].
ydouble- (input) STL vector of 21 element array to flag a double-pixel
xdouble- (input) STL vector of 7 element array to flag a double-pixel
templ- (input) the template used in the reconstruction
yrec- (output) best estimate of y-coordinate of hit in microns
sigmay- (output) best estimate of uncertainty on yrec in microns
proby- (output) probability describing goodness-of-fit for y-reco
xrec- (output) best estimate of x-coordinate of hit in microns
sigmax- (output) best estimate of uncertainty on xrec in microns
probx- (output) probability describing goodness-of-fit for x-reco
qbin- (output) index (0-4) describing the charge of the cluster [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin]
speed- (input) switch (-1-4) trading speed vs robustness -1 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability 0 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4) 1 faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density 2 faster yet, searches same range as 1 but at 1/2 density 3 fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster) 4 fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster), calculates Q probability
probQ- (output) the Vavilov-distribution-based cluster charge probability

Definition at line 1313 of file SiPixelTemplateReco.cc.

References f, PixelTempReco1D(), and HLTSiStripMonitoring_cff::speed.

1328 {
1329  // Local variables
1330  const bool deadpix = false;
1331  std::vector<std::pair<int, int> > zeropix;
1332  int nypix, nxpix;
1333  float locBx, locBz;
1334  locBx = 1.f;
1335  if (cotbeta < 0.f) {
1336  locBx = -1.f;
1337  }
1338  locBz = locBx;
1339  if (cotalpha < 0.f) {
1340  locBz = -locBx;
1341  }
1342 
1344  cotalpha,
1345  cotbeta,
1346  locBz,
1347  locBx,
1348  cluster,
1349  templ,
1350  yrec,
1351  sigmay,
1352  proby,
1353  xrec,
1354  sigmax,
1355  probx,
1356  qbin,
1357  speed,
1358  deadpix,
1359  zeropix,
1360  probQ,
1361  nypix,
1362  nxpix);
1363 
1364 } // PixelTempReco1D
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)
double f[11][100]

◆ PixelTempReco1D() [4/4]

int SiPixelTemplateReco::PixelTempReco1D ( int  id,
float  cotalpha,
float  cotbeta,
ClusMatrix cluster,
SiPixelTemplate templ,
float &  yrec,
float &  sigmay,
float &  proby,
float &  xrec,
float &  sigmax,
float &  probx,
int &  qbin,
int  speed 
)

Reconstruct the best estimate of the hit position for pixel clusters.

Parameters
id- (input) identifier of the template to use
cotalpha- (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
cotbeta- (input) the cotangent of the beta track angle (see CMS IN 2004/014)
cluster- (input) boost multi_array container of 7x21 array of pixel signals, origin of local coords (0,0) at center of pixel cluster[0][0].
ydouble- (input) STL vector of 21 element array to flag a double-pixel
xdouble- (input) STL vector of 7 element array to flag a double-pixel
templ- (input) the template used in the reconstruction
yrec- (output) best estimate of y-coordinate of hit in microns
sigmay- (output) best estimate of uncertainty on yrec in microns
proby- (output) probability describing goodness-of-fit for y-reco
xrec- (output) best estimate of x-coordinate of hit in microns
sigmax- (output) best estimate of uncertainty on xrec in microns
probx- (output) probability describing goodness-of-fit for x-reco
qbin- (output) index (0-4) describing the charge of the cluster [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin]
speed- (input) switch (0-3) trading speed vs robustness 0 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4) 1 faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density 2 faster yet, searches same range as 1 but at 1/2 density 3 fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)

Definition at line 1391 of file SiPixelTemplateReco.cc.

References f, PixelTempReco1D(), and HLTSiStripMonitoring_cff::speed.

1405 {
1406  // Local variables
1407  const bool deadpix = false;
1408  std::vector<std::pair<int, int> > zeropix;
1409  int nypix, nxpix;
1410  float locBx, locBz;
1411  locBx = 1.f;
1412  if (cotbeta < 0.f) {
1413  locBx = -1.f;
1414  }
1415  locBz = locBx;
1416  if (cotalpha < 0.f) {
1417  locBz = -locBx;
1418  }
1419  float probQ;
1420  if (speed < 0)
1421  speed = 0;
1422  if (speed > 3)
1423  speed = 3;
1424 
1426  cotalpha,
1427  cotbeta,
1428  locBz,
1429  locBx,
1430  cluster,
1431  templ,
1432  yrec,
1433  sigmay,
1434  proby,
1435  xrec,
1436  sigmax,
1437  probx,
1438  qbin,
1439  speed,
1440  deadpix,
1441  zeropix,
1442  probQ,
1443  nypix,
1444  nxpix);
1445 
1446 } // PixelTempReco1D
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)
double f[11][100]