CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Classes | Functions
SiPixelTemplateReco Namespace Reference

Classes

struct  ClusMatrix
 

Functions

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

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 the local B_z field for FPix (usually B_z<0 when cot(beta)>0 and B_z>0 when cot(beta)<0
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

Definition at line 125 of file SiPixelTemplateReco.cc.

References 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(), delta, SiPixelTemplate::dxone(), SiPixelTemplate::dxtwo(), SiPixelTemplate::dyone(), SiPixelTemplate::dytwo(), ENDL, Exception, f, VVIObjF::fcn(), Gamma, i, SiPixelTemplate::interpolate(), j, relval_2017::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 PixelCPETemplateReco::localPosition(), and PixelTempReco2D().

130 {
131  // Local variables
132  int i, j, k, minbin, binl, binh, binq, midpix, fypix, nypix, lypix, logypx;
133  int fxpix, nxpix, lxpix, logxpx, shifty, shiftx, nyzero[TYSIZE];
134  int nclusx, nclusy;
135  int deltaj, jmin, jmax, fxbin, lxbin, fybin, lybin, djy, djx;
136  //int fypix2D, lypix2D, fxpix2D, lxpix2D;
137  float sythr, sxthr, rnorm, delta, sigma, sigavg, pseudopix, qscale, q50;
138  float ss2, ssa, sa2, ssba, saba, sba2, rat, fq, qtotal, qpixel;
139  float originx, originy, qfy, qly, qfx, qlx, bias, maxpix, minmax;
140  double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
141  double hchi2, hndof, prvav, mpv, sigmaQ, kappa, xvav, beta2;
142  float ytemp[41][BYSIZE], xtemp[41][BXSIZE], ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
143  float chi2ybin[41], chi2xbin[41], ysig2[BYSIZE], xsig2[BXSIZE];
144  float yw2[BYSIZE], xw2[BXSIZE], ysw[BYSIZE], xsw[BXSIZE];
145  bool yd[BYSIZE], xd[BXSIZE], anyyd, anyxd, calc_probQ, use_VVIObj;
146  float ysize, xsize;
147  const float probmin={1.110223e-16};
148  const float probQmin={1.e-5};
149 
150 // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
151 
152  const double mean1pix={0.100}, chi21min={0.160};
153 
154 // First, interpolate the template needed to analyze this cluster
155 // check to see of the track direction is in the physical range of the loaded template
156 
157  if(!templ.interpolate(id, cotalpha, cotbeta, locBz)) {
158  if (theVerboseLevel > 2) {LOGDEBUG("SiPixelTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha << ", cot(beta) = " << cotbeta<< ", local B_z = " << locBz << ", template ID = " << id << ", no reconstruction performed" << ENDL;}
159  return 20;
160  }
161 
162 // Check to see if Q probability is selected
163 
164  calc_probQ = false;
165  use_VVIObj = false;
166  if(speed < 0) {
167  calc_probQ = true;
168  if(speed < -1) use_VVIObj = true;
169  speed = 0;
170  }
171 
172  if(speed > 3) {
173  calc_probQ = true;
174  if(speed < 5) use_VVIObj = true;
175  speed = 3;
176  }
177 
178 // Get pixel dimensions from the template (to allow multiple detectors in the future)
179 
180  xsize = templ.xsize();
181  ysize = templ.ysize();
182 
183 // Define size of pseudopixel
184 
185  q50 = templ.s50();
186  pseudopix = 0.2f*q50;
187 
188 // Get charge scaling factor
189 
190  qscale = templ.qscale();
191 
192 // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
193 
194  nclusx = cluster.mrow;
195  nclusy = (int)cluster.mcol;
196  auto const xdouble = cluster.xdouble;
197  auto const ydouble = cluster.ydouble;
198 
199 // First, rescale all pixel charges and compute total charge
200  qtotal = 0.f;
201  for(i=0; i<nclusx*nclusy; ++i) {
202  cluster.matrix[i] *= qscale;
203  qtotal +=cluster.matrix[i];
204  }
205 // Next, sum the total charge and "decapitate" big pixels
206  minmax = templ.pixmax();
207  for(j=0; j<nclusx; ++j)
208  for(i=0; i<nclusy; ++i) {
209  maxpix = minmax;
210  if(ydouble[i]) {maxpix *=2.f;}
211  if(cluster(j,i) > maxpix) {cluster(j,i) = maxpix;}
212  }
213 
214 
215 // Do the cluster repair here
216 
217  if(deadpix) {
218  fypix = BYM3; lypix = -1;
219  memset(nyzero, 0, TYSIZE * sizeof(int));
220  memset(ysum, 0, BYSIZE * sizeof(float));
221  for(i=0; i<nclusy; ++i) {
222 // Do preliminary cluster projection in y
223  for(j=0; j<nclusx; ++j) {
224  ysum[i] += cluster(j,i);
225  }
226  if(ysum[i] > 0.f) {
227 // identify ends of cluster to determine what the missing charge should be
228  if(i < fypix) {fypix = i;}
229  if(i > lypix) {lypix = i;}
230  }
231  }
232 
233 // Now loop over dead pixel list and "fix" everything
234 
235 //First see if the cluster ends are redefined and that we have only one dead pixel per column
236 
237  std::vector<std::pair<int, int> >::const_iterator zeroIter = zeropix.begin(), zeroEnd = zeropix.end();
238  for ( ; zeroIter != zeroEnd; ++zeroIter ) {
239  i = zeroIter->second;
240  if(i<0 || i>TYSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column y-index " << i << ", no reconstruction performed" << ENDL;
241  return 11;}
242 
243 // count the number of dead pixels in each column
244  ++nyzero[i];
245 // allow them to redefine the cluster ends
246  if(i < fypix) {fypix = i;}
247  if(i > lypix) {lypix = i;}
248  }
249 
250  nypix = lypix-fypix+1;
251 
252 // 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
253 
254  for (zeroIter = zeropix.begin(); zeroIter != zeroEnd; ++zeroIter ) {
255  i = zeroIter->second; j = zeroIter->first;
256  if(j<0 || j>TXSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column x-index " << j << ", no reconstruction performed" << ENDL;
257  return 12;}
258  if((i == fypix || i == lypix) && nypix > 1) {maxpix = templ.symax()/2.;} else {maxpix = templ.symax();}
259  if(ydouble[i]) {maxpix *=2.;}
260  if(nyzero[i] > 0 && nyzero[i] < 3) {qpixel = (maxpix - ysum[i])/(float)nyzero[i];} else {qpixel = 1.;}
261  if(qpixel < 1.) {qpixel = 1.;}
262  cluster(j,i) = qpixel;
263 // Adjust the total cluster charge to reflect the charge of the "repaired" cluster
264  qtotal += qpixel;
265  }
266 // End of cluster repair section
267  }
268 
269 // Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container
270 
271  for(i=0; i<BYSIZE; ++i) { ysum[i] = 0.f; yd[i] = false;}
272  k=0;
273  anyyd = false;
274  for(i=0; i<nclusy; ++i) {
275  for(j=0; j<nclusx; ++j) {
276  ysum[k] += cluster(j,i);
277  }
278 
279 // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
280 
281  if(ydouble[i]) {
282  ysum[k] /= 2.f;
283  ysum[k+1] = ysum[k];
284  yd[k] = true;
285  yd[k+1] = false;
286  k=k+2;
287  anyyd = true;
288  } else {
289  yd[k] = false;
290  ++k;
291  }
292  if(k > BYM1) {break;}
293  }
294 
295 // Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container
296 
297  for(i=0; i<BXSIZE; ++i) { xsum[i] = 0.f; xd[i] = false;}
298  k=0;
299  anyxd = false;
300  for(j=0; j<nclusx; ++j) {
301  for(i=0; i<nclusy; ++i) {
302  xsum[k] += cluster(j,i);
303  }
304 
305 // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
306 
307  if(xdouble[j]) {
308  xsum[k] /= 2.;
309  xsum[k+1] = xsum[k];
310  xd[k]=true;
311  xd[k+1]=false;
312  k=k+2;
313  anyxd = true;
314  } else {
315  xd[k]=false;
316  ++k;
317  }
318  if(k > BXM1) {break;}
319  }
320 
321 // next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx
322 
323  fypix=-1;
324  nypix=0;
325  lypix=0;
326  logypx=0;
327  for(i=0; i<BYSIZE; ++i) {
328  if(ysum[i] > 0.f) {
329  if(fypix == -1) {fypix = i;}
330  if(!yd[i]) {
331  ysort[logypx] = ysum[i];
332  ++logypx;
333  }
334  ++nypix;
335  lypix = i;
336  }
337  }
338 
339 // dlengthy = (float)nypix - templ.clsleny();
340 
341 // Make sure cluster is continuous
342 
343  if((lypix-fypix+1) != nypix || nypix == 0) {
344  LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
345  if (theVerboseLevel > 2) {
346  LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
347  for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}
348  LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
349  }
350 
351  return 1;
352  }
353 
354 // If cluster is longer than max template size, technique fails
355 
356  if(nypix > TYSIZE) {
357  LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster is larger than maximum template size" << ENDL;
358  if (theVerboseLevel > 2) {
359  LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
360  for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}
361  LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
362  }
363 
364  return 6;
365  }
366 
367 // Remember these numbers for later
368 
369  //fypix2D = fypix;
370  //lypix2D = lypix;
371 
372 // next, center the cluster on template center if necessary
373 
374  midpix = (fypix+lypix)/2;
375  shifty = templ.cytemp() - midpix;
376  if(shifty > 0) {
377  for(i=lypix; i>=fypix; --i) {
378  ysum[i+shifty] = ysum[i];
379  ysum[i] = 0.;
380  yd[i+shifty] = yd[i];
381  yd[i] = false;
382  }
383  } else if (shifty < 0) {
384  for(i=fypix; i<=lypix; ++i) {
385  ysum[i+shifty] = ysum[i];
386  ysum[i] = 0.;
387  yd[i+shifty] = yd[i];
388  yd[i] = false;
389  }
390  }
391  lypix +=shifty;
392  fypix +=shifty;
393 
394 // If the cluster boundaries are OK, add pesudopixels, otherwise quit
395 
396  if(fypix > 1 && fypix < BYM2) {
397  ysum[fypix-1] = pseudopix;
398  ysum[fypix-2] = pseudopix;
399  } else {return 8;}
400  if(lypix > 1 && lypix < BYM2) {
401  ysum[lypix+1] = pseudopix;
402  ysum[lypix+2] = pseudopix;
403  } else {return 8;}
404 
405 // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
406 
407  if(ydouble[0]) {
408  originy = -0.5f;
409  } else {
410  originy = 0.f;
411  }
412 
413 // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx
414 
415  fxpix=-1;
416  nxpix=0;
417  lxpix=0;
418  logxpx=0;
419  for(i=0; i<BXSIZE; ++i) {
420  if(xsum[i] > 0.) {
421  if(fxpix == -1) {fxpix = i;}
422  if(!xd[i]) {
423  xsort[logxpx] = xsum[i];
424  ++logxpx;
425  }
426  ++nxpix;
427  lxpix = i;
428  }
429  }
430 
431 // dlengthx = (float)nxpix - templ.clslenx();
432 
433 // Make sure cluster is continuous
434 
435  if((lxpix-fxpix+1) != nxpix) {
436 
437  LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
438  if (theVerboseLevel > 2) {
439  LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
440  for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}
441  LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
442  }
443 
444  return 2;
445  }
446 
447 // If cluster is longer than max template size, technique fails
448 
449  if(nxpix > TXSIZE) {
450 
451  LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
452  if (theVerboseLevel > 2) {
453  LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
454  for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}
455  LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
456  }
457 
458  return 7;
459  }
460 
461 // Remember these numbers for later
462 
463  //fxpix2D = fxpix;
464  //lxpix2D = lxpix;
465 
466 // next, center the cluster on template center if necessary
467 
468  midpix = (fxpix+lxpix)/2;
469  shiftx = templ.cxtemp() - midpix;
470  if(shiftx > 0) {
471  for(i=lxpix; i>=fxpix; --i) {
472  xsum[i+shiftx] = xsum[i];
473  xsum[i] = 0.;
474  xd[i+shiftx] = xd[i];
475  xd[i] = false;
476  }
477  } else if (shiftx < 0) {
478  for(i=fxpix; i<=lxpix; ++i) {
479  xsum[i+shiftx] = xsum[i];
480  xsum[i] = 0.;
481  xd[i+shiftx] = xd[i];
482  xd[i] = false;
483  }
484  }
485  lxpix +=shiftx;
486  fxpix +=shiftx;
487 
488 // If the cluster boundaries are OK, add pesudopixels, otherwise quit
489 
490  if(fxpix > 1 && fxpix < BXM2) {
491  xsum[fxpix-1] = pseudopix;
492  xsum[fxpix-2] = pseudopix;
493  } else {return 9;}
494  if(lxpix > 1 && lxpix < BXM2) {
495  xsum[lxpix+1] = pseudopix;
496  xsum[lxpix+2] = pseudopix;
497  } else {return 9;}
498 
499 // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
500 
501  if(xdouble[0]) {
502  originx = -0.5f;
503  } else {
504  originx = 0.f;
505  }
506 
507 // uncertainty and final corrections depend upon total charge bin
508 
509  fq = qtotal/templ.qavg();
510  if(fq > 1.5f) {
511  binq=0;
512  } else {
513  if(fq > 1.0f) {
514  binq=1;
515  } else {
516  if(fq > 0.85f) {
517  binq=2;
518  } else {
519  binq=3;
520  }
521  }
522  }
523 
524 // Return the charge bin via the parameter list unless the charge is too small (then flag it)
525 
526  qbin = binq;
527  if(!deadpix && qtotal < 0.95f*templ.qmin()) {qbin = 5;} else {
528  if(!deadpix && qtotal < 0.95f*templ.qmin(1)) {qbin = 4;}
529  }
530  if (theVerboseLevel > 9) {
531  LOGDEBUG("SiPixelTemplateReco") <<
532  "ID = " << id <<
533  " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta <<
534  " nclusx = " << nclusx << " nclusy = " << nclusy << ENDL;
535  }
536 
537 
538 // Next, copy the y- and x-templates to local arrays
539 
540 // First, decide on chi^2 min search parameters
541 
542 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
543  if(speed < 0 || speed > 3) {
544  throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::PixelTempReco2D called with illegal speed = " << speed << std::endl;
545  }
546 #else
547  assert(speed >= 0 && speed < 4);
548 #endif
549  fybin = 2; lybin = 38; fxbin = 2; lxbin = 38; djy = 1; djx = 1;
550  if(speed > 0) {
551  fybin = 8; lybin = 32;
552  if(yd[fypix]) {fybin = 4; lybin = 36;}
553  if(lypix > fypix) {
554  if(yd[lypix-1]) {fybin = 4; lybin = 36;}
555  }
556  fxbin = 8; lxbin = 32;
557  if(xd[fxpix]) {fxbin = 4; lxbin = 36;}
558  if(lxpix > fxpix) {
559  if(xd[lxpix-1]) {fxbin = 4; lxbin = 36;}
560  }
561  }
562 
563  if(speed > 1) {
564  djy = 2; djx = 2;
565  if(speed > 2) {
566  if(!anyyd) {djy = 4;}
567  if(!anyxd) {djx = 4;}
568  }
569  }
570 
571  if (theVerboseLevel > 9) {
572  LOGDEBUG("SiPixelTemplateReco") <<
573  "fypix " << fypix << " lypix = " << lypix <<
574  " fybin = " << fybin << " lybin = " << lybin <<
575  " djy = " << djy << " logypx = " << logypx << ENDL;
576  LOGDEBUG("SiPixelTemplateReco") <<
577  "fxpix " << fxpix << " lxpix = " << lxpix <<
578  " fxbin = " << fxbin << " lxbin = " << lxbin <<
579  " djx = " << djx << " logxpx = " << logxpx << ENDL;
580  }
581 
582 // Now do the copies
583 
584  templ.ytemp(fybin, lybin, ytemp);
585 
586  templ.xtemp(fxbin, lxbin, xtemp);
587 
588 // Do the y-reconstruction first
589 
590 // Apply the first-pass template algorithm to all clusters
591 
592 // Modify the template if double pixels are present
593 
594  if(nypix > logypx) {
595  i=fypix;
596  while(i < lypix) {
597  if(yd[i] && !yd[i+1]) {
598  for(j=fybin; j<=lybin; ++j) {
599 
600 // Sum the adjacent cells and put the average signal in both
601 
602  sigavg = (ytemp[j][i] + ytemp[j][i+1])/2.f;
603  ytemp[j][i] = sigavg;
604  ytemp[j][i+1] = sigavg;
605  }
606  i += 2;
607  } else {
608  ++i;
609  }
610  }
611  }
612 
613 // Define the maximum signal to allow before de-weighting a pixel
614 
615  sythr = 1.1*(templ.symax());
616 
617 // Make sure that there will be at least two pixels that are not de-weighted
618 
619  std::sort(&ysort[0], &ysort[logypx]);
620  if(logypx == 1) {sythr = 1.01f*ysort[0];} else {
621  if (ysort[1] > sythr) { sythr = 1.01f*ysort[1]; }
622  }
623 
624 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
625 
626 // for(i=0; i<BYSIZE; ++i) { ysig2[i] = 0.;}
627  templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
628 
629 // Find the template bin that minimizes the Chi^2
630 
631  chi2ymin = 1.e15;
632  for(i=fybin; i<=lybin; ++i) { chi2ybin[i] = -1.e15f;}
633  ss2 = 0.f;
634  for(i=fypix-2; i<=lypix+2; ++i) {
635  yw2[i] = 1.f/ysig2[i];
636  ysw[i] = ysum[i]*yw2[i];
637  ss2 += ysum[i]*ysw[i];
638  }
639 
640  minbin = -1;
641  deltaj = djy;
642  jmin = fybin;
643  jmax = lybin;
644  while(deltaj > 0) {
645  for(j=jmin; j<=jmax; j+=deltaj) {
646  if(chi2ybin[j] < -100.f) {
647  ssa = 0.f;
648  sa2 = 0.f;
649  for(i=fypix-2; i<=lypix+2; ++i) {
650  ssa += ysw[i]*ytemp[j][i];
651  sa2 += ytemp[j][i]*ytemp[j][i]*yw2[i];
652  }
653  rat=ssa/ss2;
654  if(rat <= 0.f) {LOGERROR("SiPixelTemplateReco") << "illegal chi2ymin normalization (1) = " << rat << ENDL; rat = 1.;}
655  chi2ybin[j]=ss2-2.f*ssa/rat+sa2/(rat*rat);
656  }
657  if(chi2ybin[j] < chi2ymin) {
658  chi2ymin = chi2ybin[j];
659  minbin = j;
660  }
661  }
662  deltaj /= 2;
663  if(minbin > fybin) {jmin = minbin - deltaj;} else {jmin = fybin;}
664  if(minbin < lybin) {jmax = minbin + deltaj;} else {jmax = lybin;}
665  }
666 
667  if (theVerboseLevel > 9) {
668  LOGDEBUG("SiPixelTemplateReco") <<
669  "minbin " << minbin << " chi2ymin = " << chi2ymin << ENDL;
670  }
671 
672 // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
673 
674  if(logypx == 1) {
675 
676  if(nypix ==1) {
677  delta = templ.dyone();
678  sigma = templ.syone();
679  } else {
680  delta = templ.dytwo();
681  sigma = templ.sytwo();
682  }
683 
684  yrec = 0.5f*(fypix+lypix-2*shifty+2.f*originy)*ysize-delta;
685  if(sigma <= 0.f) {
686  sigmay = 43.3f;
687  } else {
688  sigmay = sigma;
689  }
690 
691 // Do probability calculation for one-pixel clusters
692 
693  chi21max = fmax(chi21min, (double)templ.chi2yminone());
694  chi2ymin -=chi21max;
695  if(chi2ymin < 0.) {chi2ymin = 0.;}
696 // proby = gsl_cdf_chisq_Q(chi2ymin, mean1pix);
697  meany = fmax(mean1pix, (double)templ.chi2yavgone());
698  hchi2 = chi2ymin/2.; hndof = meany/2.;
699  proby = 1. - TMath::Gamma(hndof, hchi2);
700 
701  } else {
702 
703 // For cluster > 1 pix, make the second, interpolating pass with the templates
704 
705  binl = minbin - 1;
706  binh = binl + 2;
707  if(binl < fybin) { binl = fybin;}
708  if(binh > lybin) { binh = lybin;}
709  ssa = 0.;
710  sa2 = 0.;
711  ssba = 0.;
712  saba = 0.;
713  sba2 = 0.;
714  for(i=fypix-2; i<=lypix+2; ++i) {
715  ssa += ysw[i]*ytemp[binl][i];
716  sa2 += ytemp[binl][i]*ytemp[binl][i]*yw2[i];
717  ssba += ysw[i]*(ytemp[binh][i] - ytemp[binl][i]);
718  saba += ytemp[binl][i]*(ytemp[binh][i] - ytemp[binl][i])*yw2[i];
719  sba2 += (ytemp[binh][i] - ytemp[binl][i])*(ytemp[binh][i] - ytemp[binl][i])*yw2[i];
720  }
721 
722 // rat is the fraction of the "distance" from template a to template b
723 
724  rat=(ssba*ssa-ss2*saba)/(ss2*sba2-ssba*ssba);
725  if(rat < 0.f) {rat=0.f;}
726  if(rat > 1.f) {rat=1.0f;}
727  rnorm = (ssa+rat*ssba)/ss2;
728 
729 // Calculate the charges in the first and last pixels
730 
731  qfy = ysum[fypix];
732  if(yd[fypix]) {qfy+=ysum[fypix+1];}
733  if(logypx > 1) {
734  qly=ysum[lypix];
735  if(yd[lypix-1]) {qly+=ysum[lypix-1];}
736  } else {
737  qly = qfy;
738  }
739 
740 // Now calculate the mean bias correction and uncertainties
741 
742  float qyfrac = (qfy-qly)/(qfy+qly);
743  bias = templ.yflcorr(binq,qyfrac)+templ.yavg(binq);
744 
745 // uncertainty and final correction depend upon charge bin
746 
747  yrec = (0.125f*binl+BHY-2.5f+rat*(binh-binl)*0.125f-(float)shifty+originy)*ysize - bias;
748  sigmay = templ.yrms(binq);
749 
750 // Do goodness of fit test in y
751 
752  if(rnorm <= 0.) {LOGERROR("SiPixelTemplateReco") << "illegal chi2y normalization (2) = " << rnorm << ENDL; rnorm = 1.;}
753  chi2y=ss2-2./rnorm*ssa-2./rnorm*rat*ssba+(sa2+2.*rat*saba+rat*rat*sba2)/(rnorm*rnorm)-templ.chi2ymin(binq);
754  if(chi2y < 0.0) {chi2y = 0.0;}
755  meany = templ.chi2yavg(binq);
756  if(meany < 0.01) {meany = 0.01;}
757 // gsl function that calculates the chi^2 tail prob for non-integral dof
758 // proby = gsl_cdf_chisq_Q(chi2y, meany);
759 // proby = ROOT::Math::chisquared_cdf_c(chi2y, meany);
760  hchi2 = chi2y/2.; hndof = meany/2.;
761  proby = 1. - TMath::Gamma(hndof, hchi2);
762  }
763 
764 // Do the x-reconstruction next
765 
766 // Apply the first-pass template algorithm to all clusters
767 
768 // Modify the template if double pixels are present
769 
770  if(nxpix > logxpx) {
771  i=fxpix;
772  while(i < lxpix) {
773  if(xd[i] && !xd[i+1]) {
774  for(j=fxbin; j<=lxbin; ++j) {
775 
776 // Sum the adjacent cells and put the average signal in both
777 
778  sigavg = (xtemp[j][i] + xtemp[j][i+1])/2.f;
779  xtemp[j][i] = sigavg;
780  xtemp[j][i+1] = sigavg;
781  }
782  i += 2;
783  } else {
784  ++i;
785  }
786  }
787  }
788 
789 // Define the maximum signal to allow before de-weighting a pixel
790 
791  sxthr = 1.1f*templ.sxmax();
792 
793 // Make sure that there will be at least two pixels that are not de-weighted
794  std::sort(&xsort[0], &xsort[logxpx]);
795  if(logxpx == 1) {sxthr = 1.01f*xsort[0];} else {
796  if (xsort[1] > sxthr) { sxthr = 1.01f*xsort[1]; }
797  }
798 
799 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
800 
801 // for(i=0; i<BXSIZE; ++i) { xsig2[i] = 0.; }
802  templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
803 
804 // Find the template bin that minimizes the Chi^2
805 
806  chi2xmin = 1.e15;
807  for(i=fxbin; i<=lxbin; ++i) { chi2xbin[i] = -1.e15f;}
808  ss2 = 0.f;
809  for(i=fxpix-2; i<=lxpix+2; ++i) {
810  xw2[i] = 1.f/xsig2[i];
811  xsw[i] = xsum[i]*xw2[i];
812  ss2 += xsum[i]*xsw[i];
813  }
814  minbin = -1;
815  deltaj = djx;
816  jmin = fxbin;
817  jmax = lxbin;
818  while(deltaj > 0) {
819  for(j=jmin; j<=jmax; j+=deltaj) {
820  if(chi2xbin[j] < -100.f) {
821  ssa = 0.f;
822  sa2 = 0.f;
823  for(i=fxpix-2; i<=lxpix+2; ++i) {
824  ssa += xsw[i]*xtemp[j][i];
825  sa2 += xtemp[j][i]*xtemp[j][i]*xw2[i];
826  }
827  rat=ssa/ss2;
828  if(rat <= 0.f) {LOGERROR("SiPixelTemplateReco") << "illegal chi2xmin normalization (1) = " << rat << ENDL; rat = 1.;}
829  chi2xbin[j]=ss2-2.f*ssa/rat+sa2/(rat*rat);
830  }
831  if(chi2xbin[j] < chi2xmin) {
832  chi2xmin = chi2xbin[j];
833  minbin = j;
834  }
835  }
836  deltaj /= 2;
837  if(minbin > fxbin) {jmin = minbin - deltaj;} else {jmin = fxbin;}
838  if(minbin < lxbin) {jmax = minbin + deltaj;} else {jmax = lxbin;}
839  }
840 
841  if (theVerboseLevel > 9) {
842  LOGDEBUG("SiPixelTemplateReco") <<
843  "minbin " << minbin << " chi2xmin = " << chi2xmin << ENDL;
844  }
845 
846 // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
847 
848  if(logxpx == 1) {
849 
850  if(nxpix ==1) {
851  delta = templ.dxone();
852  sigma = templ.sxone();
853  } else {
854  delta = templ.dxtwo();
855  sigma = templ.sxtwo();
856  }
857  xrec = 0.5*(fxpix+lxpix-2*shiftx+2.*originx)*xsize-delta;
858  if(sigma <= 0.) {
859  sigmax = 28.9;
860  } else {
861  sigmax = sigma;
862  }
863 
864 // Do probability calculation for one-pixel clusters
865 
866  chi21max = fmax(chi21min, (double)templ.chi2xminone());
867  chi2xmin -=chi21max;
868  if(chi2xmin < 0.) {chi2xmin = 0.;}
869  meanx = fmax(mean1pix, (double)templ.chi2xavgone());
870  hchi2 = chi2xmin/2.; hndof = meanx/2.;
871  probx = 1. - TMath::Gamma(hndof, hchi2);
872 
873  } else {
874 
875 // Now make the second, interpolating pass with the templates
876 
877  binl = minbin - 1;
878  binh = binl + 2;
879  if(binl < fxbin) { binl = fxbin;}
880  if(binh > lxbin) { binh = lxbin;}
881  ssa = 0.;
882  sa2 = 0.;
883  ssba = 0.;
884  saba = 0.;
885  sba2 = 0.;
886  for(i=fxpix-2; i<=lxpix+2; ++i) {
887  ssa += xsw[i]*xtemp[binl][i];
888  sa2 += xtemp[binl][i]*xtemp[binl][i]*xw2[i];
889  ssba += xsw[i]*(xtemp[binh][i] - xtemp[binl][i]);
890  saba += xtemp[binl][i]*(xtemp[binh][i] - xtemp[binl][i])*xw2[i];
891  sba2 += (xtemp[binh][i] - xtemp[binl][i])*(xtemp[binh][i] - xtemp[binl][i])*xw2[i];
892  }
893 
894 // rat is the fraction of the "distance" from template a to template b
895 
896  rat=(ssba*ssa-ss2*saba)/(ss2*sba2-ssba*ssba);
897  if(rat < 0.f) {rat=0.f;}
898  if(rat > 1.f) {rat=1.0f;}
899  rnorm = (ssa+rat*ssba)/ss2;
900 
901 // Calculate the charges in the first and last pixels
902 
903  qfx = xsum[fxpix];
904  if(xd[fxpix]) {qfx+=xsum[fxpix+1];}
905  if(logxpx > 1) {
906  qlx=xsum[lxpix];
907  if(xd[lxpix-1]) {qlx+=xsum[lxpix-1];}
908  } else {
909  qlx = qfx;
910  }
911 
912 // Now calculate the mean bias correction and uncertainties
913 
914  float qxfrac = (qfx-qlx)/(qfx+qlx);
915  bias = templ.xflcorr(binq,qxfrac)+templ.xavg(binq);
916 
917 // uncertainty and final correction depend upon charge bin
918 
919  xrec = (0.125f*binl+BHX-2.5f+rat*(binh-binl)*0.125f-(float)shiftx+originx)*xsize - bias;
920  sigmax = templ.xrms(binq);
921 
922 // Do goodness of fit test in x
923 
924  if(rnorm <= 0.) {LOGERROR("SiPixelTemplateReco") << "illegal chi2x normalization (2) = " << rnorm << ENDL; rnorm = 1.;}
925  chi2x=ss2-2./rnorm*ssa-2./rnorm*rat*ssba+(sa2+2.*rat*saba+rat*rat*sba2)/(rnorm*rnorm)-templ.chi2xmin(binq);
926  if(chi2x < 0.0) {chi2x = 0.0;}
927  meanx = templ.chi2xavg(binq);
928  if(meanx < 0.01) {meanx = 0.01;}
929 // gsl function that calculates the chi^2 tail prob for non-integral dof
930 // probx = gsl_cdf_chisq_Q(chi2x, meanx);
931 // probx = ROOT::Math::chisquared_cdf_c(chi2x, meanx, trx0);
932  hchi2 = chi2x/2.; hndof = meanx/2.;
933  probx = 1. - TMath::Gamma(hndof, hchi2);
934  }
935 
936 // Don't return exact zeros for the probability
937 
938  if(proby < probmin) {proby = probmin;}
939  if(probx < probmin) {probx = probmin;}
940 
941 // Decide whether to generate a cluster charge probability
942 
943  if(calc_probQ) {
944 
945 // Calculate the Vavilov probability that the cluster charge is OK
946 
947  templ.vavilov_pars(mpv, sigmaQ, kappa);
948 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
949  if((sigmaQ <=0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
950  throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::Vavilov parameters mpv/sigmaQ/kappa = " << mpv << "/" << sigmaQ << "/" << kappa << std::endl;
951  }
952 #else
953  assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
954 #endif
955  xvav = ((double)qtotal-mpv)/sigmaQ;
956  beta2 = 1.;
957  if(use_VVIObj) {
958 // VVIObj is a private port of CERNLIB VVIDIS
959  VVIObjF vvidist(kappa, beta2, 1);
960  prvav = vvidist.fcn(xvav);
961  } else {
962 // Use faster but less accurate TMath Vavilov distribution function
963  prvav = TMath::VavilovI(xvav, kappa, beta2);
964  }
965 // Change to upper tail probability
966 // if(prvav > 0.5) prvav = 1. - prvav;
967 // probQ = (float)(2.*prvav);
968  probQ = 1. - prvav;
969  if(probQ < probQmin) {probQ = probQmin;}
970  } else {
971  probQ = -1;
972  }
973 
974  return 0;
975 } // PixelTempReco2D
dbl * delta
Definition: mlp_gen.cc:36
int i
Definition: DBlmapReader.cc:9
#define BXSIZE
float chi2xminone()
//!&lt; minimum of x chi^2 for 1 pixel clusters
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
assert(m_qm.get())
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()
//!&lt; minimum of y chi^2 for 1 pixel clusters
bool interpolate(int id, float cotalpha, float cotbeta, float locBz)
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)
int j
Definition: DBlmapReader.cc:9
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()
//!&lt; 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()
//!&lt; average x chi^2 for 1 pixel clusters
if(dp >Float(M_PI)) dp-
float pixmax()
maximum pixel charge
dbl * Gamma
Definition: mlp_gen.cc:38
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::PixelTempReco2D ( int  id,
float  cotalpha,
float  cotbeta,
float  locBz,
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 the local B_z field for FPix (usually B_z<0 when cot(beta)>0 and B_z>0 when cot(beta)<0
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 1006 of file SiPixelTemplateReco.cc.

References PixelTempReco2D().

1011 {
1012  // Local variables
1013  const bool deadpix = false;
1014  std::vector<std::pair<int, int> > zeropix;
1015 
1016  return SiPixelTemplateReco::PixelTempReco2D(id, cotalpha, cotbeta, locBz, cluster, templ,
1017  yrec, sigmay, proby, xrec, sigmax, probx, qbin, speed, deadpix, zeropix, probQ);
1018 
1019 } // PixelTempReco2D
int PixelTempReco2D(int id, float cotalpha, float cotbeta, float locBz, 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 SiPixelTemplateReco::PixelTempReco2D ( 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 1049 of file SiPixelTemplateReco.cc.

References PixelTempReco2D().

1054 {
1055  // Local variables
1056  const bool deadpix = false;
1057  std::vector<std::pair<int, int> > zeropix;
1058  float locBz = -1.f;
1059  if(cotbeta < 0.) {locBz = -locBz;}
1060 
1061  return SiPixelTemplateReco::PixelTempReco2D(id, cotalpha, cotbeta, locBz, cluster, templ,
1062  yrec, sigmay, proby, xrec, sigmax, probx, qbin, speed, deadpix, zeropix, probQ);
1063 
1064 } // PixelTempReco2D
int PixelTempReco2D(int id, float cotalpha, float cotbeta, float locBz, 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 SiPixelTemplateReco::PixelTempReco2D ( 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 1092 of file SiPixelTemplateReco.cc.

References PixelTempReco2D().

1096 {
1097  // Local variables
1098  const bool deadpix = false;
1099  std::vector<std::pair<int, int> > zeropix;
1100  float locBz = -1.f;
1101  if(cotbeta < 0.) {locBz = -locBz;}
1102  float probQ;
1103  if(speed < 0) speed = 0;
1104  if(speed > 3) speed = 3;
1105 
1106  return SiPixelTemplateReco::PixelTempReco2D(id, cotalpha, cotbeta, locBz, cluster, templ,
1107  yrec, sigmay, proby, xrec, sigmax, probx, qbin, speed, deadpix, zeropix, probQ);
1108 
1109 } // PixelTempReco2D
int PixelTempReco2D(int id, float cotalpha, float cotbeta, float locBz, 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)