CMS 3D CMS Logo
Classes | Functions
SiPixelTemplateReco Namespace Reference

Classes

struct  ClusMatrix
 

Functions

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)
 
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, 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)
 

Function Documentation

◆ PixelTempReco1D() [1/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 1390 of file SiPixelTemplateReco.cc.

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

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

◆ PixelTempReco1D() [2/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 1312 of file SiPixelTemplateReco.cc.

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

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

◆ PixelTempReco1D() [3/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 134 of file SiPixelTemplateReco.cc.

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

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, 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(), 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().

◆ PixelTempReco1D() [4/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 1238 of file SiPixelTemplateReco.cc.

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

References PixelTempReco1D(), and HLTSiStripMonitoring_cff::speed.

SiPixelTemplate::sxmax
float sxmax()
average pixel signal for x-projection of cluster
Definition: SiPixelTemplate.h:417
SiPixelTemplate::dxtwo
float dxtwo()
mean offset/correction for one double-pixel x-clusters
Definition: SiPixelTemplate.h:420
mps_fire.i
i
Definition: mps_fire.py:428
SiPixelTemplate::qavg
float qavg()
average cluster charge for this set of track angles
Definition: SiPixelTemplate.h:407
BYSIZE
#define BYSIZE
Definition: SiPixelTemplateDefs.h:25
SiPixelTemplate::pixmax
float pixmax()
maximum pixel charge
Definition: SiPixelTemplate.h:408
SiPixelTemplate::sxone
float sxone()
rms for one pixel x-clusters
Definition: SiPixelTemplate.h:419
theVerboseLevel
static const int theVerboseLevel
Definition: SiPixelTemplateReco.cc:68
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
BXSIZE
#define BXSIZE
Definition: SiPixelTemplateDefs.h:33
SiPixelTemplateReco::ClusMatrix::mrow
int mrow
Definition: SiPixelTemplateReco.h:79
SiPixelTemplate::syone
float syone()
rms for one pixel y-clusters
Definition: SiPixelTemplate.h:414
SiPixelTemplate::dyone
float dyone()
mean offset/correction for one pixel y-clusters
Definition: SiPixelTemplate.h:413
cms::cuda::assert
assert(be >=bs)
ysize
const Int_t ysize
Definition: trackSplitPlot.h:43
SiPixelTemplateReco::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)
Definition: SiPixelTemplateReco.cc:134
LOGDEBUG
#define LOGDEBUG(x)
Definition: SiPixelTemplateReco.cc:74
SiPixelTemplate::cytemp
int cytemp()
Return central pixel of y template pixels above readout threshold.
Definition: SiPixelTemplate.cc:2649
SiPixelTemplate::xtemp
void xtemp(int fxbin, int lxbin, float xtemplate[41][13+4])
Definition: SiPixelTemplate.cc:2580
SiPixelTemplate::xavg
float xavg(int i)
average x-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:482
SiPixelTemplate::ysize
float ysize()
pixel y-size (microns)
Definition: SiPixelTemplate.h:676
xsize
const Int_t xsize
Definition: trackSplitPlot.h:42
LOGERROR
#define LOGERROR(x)
Definition: SiPixelTemplateReco.cc:73
SiPixelTemplate::chi2xmin
float chi2xmin(int i)
minimum y chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:552
TXSIZE
#define TXSIZE
Definition: SiPixelTemplateDefs.h:30
BXM1
#define BXM1
Definition: SiPixelTemplateDefs.h:35
SiPixelTemplate::chi2ymin
float chi2ymin(int i)
minimum y chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:532
SiPixelTemplate::s50
float s50()
1/2 of the pixel threshold signal in electrons
Definition: SiPixelTemplate.h:410
BYM2
#define BYM2
Definition: SiPixelTemplateDefs.h:28
BYM1
#define BYM1
Definition: SiPixelTemplateDefs.h:27
SiPixelTemplate::xflcorr
float xflcorr(int binq, float qflx)
Definition: SiPixelTemplate.cc:2441
SiPixelTemplate::sytwo
float sytwo()
rms for one double-pixel y-clusters
Definition: SiPixelTemplate.h:416
BXM2
#define BXM2
Definition: SiPixelTemplateDefs.h:36
BHX
#define BHX
Definition: SiPixelTemplateDefs.h:34
SiPixelTemplate::dytwo
float dytwo()
mean offset/correction for one double-pixel y-clusters
Definition: SiPixelTemplate.h:415
dqmdumpme.k
k
Definition: dqmdumpme.py:60
SiPixelTemplate::chi2yavgone
float chi2yavgone()
//!< average y chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:657
SiPixelTemplateReco::ClusMatrix::matrix
float * matrix
Definition: SiPixelTemplateReco.h:76
SiPixelTemplate::fbin
float fbin(int i)
Return lower bound of Qbin definition.
Definition: SiPixelTemplate.h:646
BYM3
#define BYM3
Definition: SiPixelTemplateDefs.h:29
SiPixelTemplate::dxone
float dxone()
mean offset/correction for one pixel x-clusters
Definition: SiPixelTemplate.h:418
SiPixelTemplate::chi2xavgone
float chi2xavgone()
//!< average x chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:659
SiPixelTemplate::symax
float symax()
average pixel signal for y-projection of cluster
Definition: SiPixelTemplate.h:412
SiPixelTemplate::chi2yminone
float chi2yminone()
//!< minimum of y chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:658
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
dumpMFGeometry_cfg.delta
delta
Definition: dumpMFGeometry_cfg.py:25
createfilelist.int
int
Definition: createfilelist.py:10
SiPixelTemplate::chi2xavg
float chi2xavg(int i)
averaage x chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:542
TYSIZE
#define TYSIZE
Definition: SiPixelTemplateDefs.h:21
maxpix
static const int maxpix
Definition: SiPixelLorentzAngle.h:48
SiPixelTemplate::cxtemp
int cxtemp()
Return central pixel of x-template pixels above readout threshold.
Definition: SiPixelTemplate.cc:2702
SiPixelTemplate::xsize
float xsize()
pixel x-size (microns)
Definition: SiPixelTemplate.h:675
HLTSiStripMonitoring_cff.speed
speed
Definition: HLTSiStripMonitoring_cff.py:22
cms::cuda::for
for(int i=first, nt=offsets[nh];i< nt;i+=gridDim.x *blockDim.x)
Definition: HistoContainer.h:27
BHY
#define BHY
Definition: SiPixelTemplateDefs.h:26
SiPixelTemplate::chi2yavg
float chi2yavg(int i)
average y chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:522
SiPixelTemplate::qscale
float qscale()
charge scaling factor
Definition: SiPixelTemplate.h:409
ENDL
#define ENDL
Definition: SiPixelTemplateReco.cc:75
SiPixelTemplate::yrms
float yrms(int i)
average y-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:452
SiPixelTemplateReco::ClusMatrix::xdouble
bool const * xdouble
Definition: SiPixelTemplateReco.h:77
SiPixelTemplateReco::ClusMatrix::mcol
int mcol
Definition: SiPixelTemplateReco.h:79
SiPixelTemplate::xrms
float xrms(int i)
average x-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:492
SiPixelTemplateReco::ClusMatrix::ydouble
bool const * ydouble
Definition: SiPixelTemplateReco.h:78
SiPixelTemplate::yflcorr
float yflcorr(int binq, float qfly)
Definition: SiPixelTemplate.cc:2385
Exception
Definition: hltDiff.cc:245
SiPixelTemplate::sxtwo
float sxtwo()
rms for one double-pixel x-clusters
Definition: SiPixelTemplate.h:421
VVIObjF
Definition: VVIObjF.h:24
SiPixelTemplate::ytemp
void ytemp(int fybin, int lybin, float ytemplate[41][21+4])
Definition: SiPixelTemplate.cc:2504
SiPixelTemplate::xsigma2
void xsigma2(int fxpix, int lxpix, float sxthr, float xsum[13+4], float xsig2[13+4])
Definition: SiPixelTemplate.cc:2267
SiPixelTemplate::ysigma2
void ysigma2(int fypix, int lypix, float sythr, float ysum[21+4], float ysig2[21+4])
kappa
static const G4double kappa
Definition: UrbanMscModel93.cc:35
SiPixelTemplate::vavilov_pars
void vavilov_pars(double &mpv, double &sigma, double &kappa)
Definition: SiPixelTemplate.cc:3995
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
SiPixelTemplate::chi2xminone
float chi2xminone()
//!< minimum of x chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:660
SiPixelTemplate::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
Definition: SiPixelTemplate.cc:1327
SiPixelTemplate::yavg
float yavg(int i)
average y-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:442
SiPixelTemplate::qmin
float qmin()
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
Definition: SiPixelTemplate.h:422
Generated for CMSSW Reference Manual by doxygen 1.8.16