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

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

References PixelTempReco1D().

1249 {
1250  // Local variables
1251  const bool deadpix = false;
1252  std::vector<std::pair<int, int> > zeropix;
1253  int nypix, nxpix;
1254 
1256  cotalpha,
1257  cotbeta,
1258  locBz,
1259  locBx,
1260  cluster,
1261  templ,
1262  yrec,
1263  sigmay,
1264  proby,
1265  xrec,
1266  sigmax,
1267  probx,
1268  qbin,
1269  speed,
1270  deadpix,
1271  zeropix,
1272  probQ,
1273  nypix,
1274  nxpix);
1275 
1276 } // 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)
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 1306 of file SiPixelTemplateReco.cc.

References f, and PixelTempReco1D().

1321 {
1322  // Local variables
1323  const bool deadpix = false;
1324  std::vector<std::pair<int, int> > zeropix;
1325  int nypix, nxpix;
1326  float locBx, locBz;
1327  locBx = 1.f;
1328  if (cotbeta < 0.f) {
1329  locBx = -1.f;
1330  }
1331  locBz = locBx;
1332  if (cotalpha < 0.f) {
1333  locBz = -locBx;
1334  }
1335 
1337  cotalpha,
1338  cotbeta,
1339  locBz,
1340  locBx,
1341  cluster,
1342  templ,
1343  yrec,
1344  sigmay,
1345  proby,
1346  xrec,
1347  sigmax,
1348  probx,
1349  qbin,
1350  speed,
1351  deadpix,
1352  zeropix,
1353  probQ,
1354  nypix,
1355  nxpix);
1356 
1357 } // 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]
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 1384 of file SiPixelTemplateReco.cc.

References f, and PixelTempReco1D().

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