CMS 3D CMS Logo
Typedefs | Functions
SiPixelTemplateSplit Namespace Reference

Typedefs

typedef boost::multi_array< float, 2 > array_2d
 
typedef boost::multi_array< bool, 2 > array_2d_bool
 
typedef boost::multi_array< float, 3 > array_3d
 

Functions

int PixelTempSplit (int id, float cotalpha, float cotbeta, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec1, float &yrec2, float &sigmay, float &prob2y, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, float &prob2Q, bool resolve, int speed, float &dchisq, bool deadpix, std::vector< std::pair< int, int > > &zeropix, SiPixelTemplate2D &templ2D)
 
int PixelTempSplit (int id, float cotalpha, float cotbeta, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec1, float &yrec2, float &sigmay, float &prob2y, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, float &prob2Q, bool resolve, int speed, float &dchisq, SiPixelTemplate2D &templ2D)
 
int PixelTempSplit (int id, float cotalpha, float cotbeta, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec1, float &yrec2, float &sigmay, float &prob2y, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, float &prob2Q, bool resolve, float &dchisq, SiPixelTemplate2D &templ2D)
 
int PixelTempSplit (int id, float cotalpha, float cotbeta, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec1, float &yrec2, float &sigmay, float &prob2y, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, float &prob2Q, SiPixelTemplate2D &templ2D)
 
int PixelTempSplit (int id, float cotalpha, float cotbeta, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec1, float &yrec2, float &sigmay, float &prob2y, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, SiPixelTemplate2D &templ2D)
 

Typedef Documentation

◆ array_2d

typedef boost::multi_array<float, 2> SiPixelTemplateSplit::array_2d

Definition at line 47 of file SiPixelTemplateSplit.h.

◆ array_2d_bool

typedef boost::multi_array<bool, 2> SiPixelTemplateSplit::array_2d_bool

Definition at line 48 of file SiPixelTemplateSplit.h.

◆ array_3d

typedef boost::multi_array<float, 3> SiPixelTemplateSplit::array_3d

Definition at line 49 of file SiPixelTemplateSplit.h.

Function Documentation

◆ PixelTempSplit() [1/5]

int SiPixelTemplateSplit::PixelTempSplit ( int  id,
float  cotalpha,
float  cotbeta,
array_2d &  clust,
std::vector< bool > &  ydouble,
std::vector< bool > &  xdouble,
SiPixelTemplate &  templ,
float &  yrec1,
float &  yrec2,
float &  sigmay,
float &  prob2y,
float &  xrec1,
float &  xrec2,
float &  sigmax,
float &  prob2x,
int &  q2bin,
float &  prob2Q,
bool  resolve,
int  speed,
float &  dchisq,
bool  deadpix,
std::vector< std::pair< int, int > > &  zeropix,
SiPixelTemplate2D &  templ2D 
)

Reconstruct the best estimate of the hit positions 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
yrec1- (output) best estimate of first y-coordinate of hit in microns
yrec2- (output) best estimate of second y-coordinate of hit in microns
sigmay- (output) best estimate of uncertainty on yrec1 and yrec2 in microns
prob2y- (output) probability describing goodness-of-fit to a merged cluster hypothesis for y-reco
xrec1- (output) best estimate of first x-coordinate of hit in microns
xrec2- (output) best estimate of second x-coordinate of hit in microns
sigmax- (output) best estimate of uncertainty on xrec1 and xrec2 in microns
prob2x- (output) probability describing goodness-of-fit to a merged cluster hypothesis for x-reco
q2bin- (output) index (0-4) describing the charge of the cluster assuming a merged 2-hit cluster hypothesis [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]
prob2Q- (output) probability that the cluster charge is compatible with a 2-hit merging
resolve- (input) if true, use 2-D chisquare to resolve the 2-fold x-y association ambiguity (costs a factor of two in speed)
speed- (input) switch (-1->2) trading speed vs robustness -1 totally bombproof, searches the entire ranges of template bins, calculates Q probability w/ VVIObj 0 totally bombproof, searches the entire template bin range at full density (no Qprob) 1 faster, searches same range as 0 but at 1/2 density 2 fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
dchisq- (output) the delta chisquare estimator used to break the degeneracy (0 means no discrimination, larger than 0.1 is good)
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

Definition at line 98 of file SiPixelTemplateSplit.cc.

References cms::cuda::assert(), BHX, BHY, BXM1, BXM2, BXM3, BXSIZE, BYM1, BYM2, BYM3, BYSIZE, SiPixelTemplate::chi2xavgc2m(), SiPixelTemplate::chi2xavgone(), SiPixelTemplate::chi2xminone(), SiPixelTemplate::chi2yavgc2m(), SiPixelTemplate::chi2yavgone(), SiPixelTemplate::chi2yminone(), SiPixelTemplate::cxtemp(), SiPixelTemplate::cytemp(), dumpMFGeometry_cfg::delta, SiPixelTemplate::dxone(), SiPixelTemplate::dxtwo(), SiPixelTemplate::dyone(), SiPixelTemplate::dytwo(), MillePedeFileConverter_cfg::e, ENDL, Exception, f, VVIObj::fcn(), mps_fire::i, createfilelist::int, SiPixelTemplate::interpolate(), dqmiolumiharvest::j, dqmdumpme::k, hgcalLayerClusters_cfi::kappa, SiPixelTemplate2D::landau_par(), dqm-mbProfile::log, LOGDEBUG, LOGERROR, maxpix, SiStripPI::min, muonClassificationByHits_cfi::pixel, SiPixelTemplate::pixmax(), SiPixelTemplate::qavg(), SiPixelTemplate::qmin(), SiPixelTemplate::qscale(), SiPixelTemplate::s50(), SiPixelTemplate::simpletemplate2D(), jetUpdater_cfi::sort, HLTSiStripMonitoring_cff::speed, mathSSE::sqrt(), SiPixelTemplate::sxmax(), SiPixelTemplate::sxone(), SiPixelTemplate::sxtwo(), SiPixelTemplate::symax(), SiPixelTemplate::syone(), SiPixelTemplate::sytwo(), theVerboseLevel, TXSIZE, TYSIZE, SiPixelTemplate::vavilov2_pars(), SiPixelTemplate::xavgc2m(), SiPixelTemplate::xrmsc2m(), SiPixelTemplate::xsigma2(), xsize, SiPixelTemplate::xsize(), SiPixelTemplate::xtemp3d(), SiPixelTemplate::xtemp3d_int(), SiPixelTemplate2D::xytemp(), SiPixelTemplate::yavgc2m(), SiPixelTemplate::yrmsc2m(), SiPixelTemplate::ysigma2(), ysize, SiPixelTemplate::ysize(), SiPixelTemplate::ytemp3d(), and SiPixelTemplate::ytemp3d_int().

Referenced by PixelTempSplit().

122 {
123  // Local variables
124  int i, j, k, binq, midpix, fypix, nypix, lypix, logypx, lparm;
125  int fxpix, nxpix, lxpix, logxpx, shifty, shiftx;
126  int nclusx, nclusy;
127  int nybin, ycbin, nxbin, xcbin, minbinj, minbink;
128  int deltaj, jmin, jmax, kmin, kmax, km, fxbin, lxbin, fybin, lybin, djy, djx;
129  float sythr, sxthr, delta, sigma, sigavg, pseudopix, xsize, ysize, qscale, lanpar[2][5];
130  float ss2, ssa, sa2, rat, fq, qtotal, qpixel, qavg;
131  float x1p, x2p, y1p, y2p, deltachi2;
132  float originx, originy, bias, maxpix, minmax;
133  double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
134  double hchi2, hndof, sigmal1, sigmal2, mpv1, mpv2, arg1, arg2, q05, like, loglike1, loglike2;
135  double prvav, mpv, sigmaQ, kappa, xvav;
136  float ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
137  float ysig2[BYSIZE], xsig2[BXSIZE];
138  float yw2[BYSIZE], xw2[BXSIZE], ysw[BYSIZE], xsw[BXSIZE];
139  float cluster2[BXM2][BYM2], temp2d1[BXM2][BYM2], temp2d2[BXM2][BYM2];
140  bool yd[BYSIZE], xd[BXSIZE], anyyd, anyxd, any2dfail;
141  const float sqrt2x = {3.0000}, sqrt2y = {1.7000};
142  const float sqrt12 = {3.4641};
143  const float probmin = {1.110223e-16};
144  const float prob2Qmin = {1.e-5};
145  std::pair<int, int> pixel;
146 
147  // bool SiPixelTemplateSplit::SimpleTemplate2D(float cotalpha, float cotbeta, float xhit, float yhit, float thick, float lorxwidth, float lorywidth,
148  // float qavg, std::vector<bool> ydouble, std::vector<bool> xdouble, float template2d[BXM2][BYM2]);
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)) {
158  LOGDEBUG("SiPixelTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha
159  << ", cot(beta) = " << cotbeta
160  << " is not within the acceptance of template ID = " << id
161  << ", no reconstruction performed" << ENDL;
162  return 20;
163  }
164 
165  // Get pixel dimensions from the template (to allow multiple detectors in the future)
166 
167  xsize = templ.xsize();
168  ysize = templ.ysize();
169 
170  // Define size of pseudopixel
171 
172  pseudopix = templ.s50();
173  // q05 = 0.28*pseudopix;
174  q05 = 0.05f * pseudopix;
175 
176  // Get charge scaling factor
177 
178  qscale = templ.qscale();
179 
180  // Make a local copy of the cluster container so that we can muck with it
181 
182  array_2d cluster = clust;
183 
184  // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
185 
186  if (cluster.num_dimensions() != 2) {
187  LOGERROR("SiPixelTemplateReco") << "input cluster container (BOOST Multiarray) has wrong number of dimensions"
188  << ENDL;
189  return 3;
190  }
191  nclusx = (int)cluster.shape()[0];
192  nclusy = (int)cluster.shape()[1];
193  if (nclusx != (int)xdouble.size()) {
194  LOGERROR("SiPixelTemplateReco") << "input cluster container x-size is not equal to double pixel flag container size"
195  << ENDL;
196  return 4;
197  }
198  if (nclusy != (int)ydouble.size()) {
199  LOGERROR("SiPixelTemplateReco") << "input cluster container y-size is not equal to double pixel flag container size"
200  << ENDL;
201  return 5;
202  }
203 
204  // enforce maximum size
205 
206  if (nclusx > TXSIZE) {
207  nclusx = TXSIZE;
208  }
209  if (nclusy > TYSIZE) {
210  nclusy = TYSIZE;
211  }
212 
213  // First, rescale all pixel charges
214 
215  for (i = 0; i < nclusy; ++i) {
216  for (j = 0; j < nclusx; ++j) {
217  if (cluster[j][i] > 0) {
218  cluster[j][i] *= qscale;
219  }
220  }
221  }
222 
223  // Next, sum the total charge and "decapitate" big pixels
224 
225  qtotal = 0.f;
226  minmax = 2.0f * templ.pixmax();
227  for (i = 0; i < nclusy; ++i) {
228  maxpix = minmax;
229  if (ydouble[i]) {
230  maxpix *= 2.f;
231  }
232  for (j = 0; j < nclusx; ++j) {
233  qtotal += cluster[j][i];
234  if (cluster[j][i] > maxpix) {
235  cluster[j][i] = maxpix;
236  }
237  }
238  }
239 
240  // Do the cluster repair here
241 
242  if (deadpix) {
243  fypix = BYM3;
244  lypix = -1;
245  for (i = 0; i < nclusy; ++i) {
246  ysum[i] = 0.f;
247  // Do preliminary cluster projection in y
248  for (j = 0; j < nclusx; ++j) {
249  ysum[i] += cluster[j][i];
250  }
251  if (ysum[i] > 0.f) {
252  // identify ends of cluster to determine what the missing charge should be
253  if (i < fypix) {
254  fypix = i;
255  }
256  if (i > lypix) {
257  lypix = i;
258  }
259  }
260  }
261 
262  // Now loop over dead pixel list and "fix" everything
263 
264  //First see if the cluster ends are redefined and that we have only one dead pixel per column
265  int nyzero[TYSIZE]{};
266  std::vector<std::pair<int, int> >::const_iterator zeroIter = zeropix.begin(), zeroEnd = zeropix.end();
267  for (; zeroIter != zeroEnd; ++zeroIter) {
268  i = zeroIter->second;
269  if (i < 0 || i > TYSIZE - 1) {
270  LOGERROR("SiPixelTemplateReco") << "dead pixel column y-index " << i << ", no reconstruction performed" << ENDL;
271  return 11;
272  }
273 
274  // count the number of dead pixels in each column
275  ++nyzero[i];
276  // allow them to redefine the cluster ends
277  if (i < fypix) {
278  fypix = i;
279  }
280  if (i > lypix) {
281  lypix = i;
282  }
283  }
284 
285  nypix = lypix - fypix + 1;
286 
287  // 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
288 
289  for (zeroIter = zeropix.begin(); zeroIter != zeroEnd; ++zeroIter) {
290  i = zeroIter->second;
291  j = zeroIter->first;
292  if (j < 0 || j > TXSIZE - 1) {
293  LOGERROR("SiPixelTemplateReco") << "dead pixel column x-index " << j << ", no reconstruction performed" << ENDL;
294  return 12;
295  }
296  if ((i == fypix || i == lypix) && nypix > 1) {
297  maxpix = templ.symax() / 2.;
298  } else {
299  maxpix = templ.symax();
300  }
301  if (ydouble[i]) {
302  maxpix *= 2.;
303  }
304  if (nyzero[i] > 0 && nyzero[i] < 3) {
305  qpixel = (maxpix - ysum[i]) / (float)nyzero[i];
306  } else {
307  qpixel = 1.;
308  }
309  if (qpixel < 1.f) {
310  qpixel = 1.f;
311  }
312  cluster[j][i] = qpixel;
313  // Adjust the total cluster charge to reflect the charge of the "repaired" cluster
314  qtotal += qpixel;
315  }
316  // End of cluster repair section
317  }
318 
319  // Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container
320 
321  for (i = 0; i < BYSIZE; ++i) {
322  ysum[i] = 0.f;
323  yd[i] = false;
324  }
325  k = 0;
326  anyyd = false;
327  for (i = 0; i < nclusy; ++i) {
328  for (j = 0; j < nclusx; ++j) {
329  ysum[k] += cluster[j][i];
330  }
331 
332  // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
333 
334  if (ydouble[i]) {
335  ysum[k] /= 2.f;
336  ysum[k + 1] = ysum[k];
337  yd[k] = true;
338  yd[k + 1] = false;
339  k = k + 2;
340  anyyd = true;
341  } else {
342  yd[k] = false;
343  ++k;
344  }
345  if (k > BYM1) {
346  break;
347  }
348  }
349 
350  // Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container
351 
352  for (i = 0; i < BXSIZE; ++i) {
353  xsum[i] = 0.f;
354  xd[i] = false;
355  }
356  k = 0;
357  anyxd = false;
358  for (j = 0; j < nclusx; ++j) {
359  for (i = 0; i < nclusy; ++i) {
360  xsum[k] += cluster[j][i];
361  }
362 
363  // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
364 
365  if (xdouble[j]) {
366  xsum[k] /= 2.f;
367  xsum[k + 1] = xsum[k];
368  xd[k] = true;
369  xd[k + 1] = false;
370  k = k + 2;
371  anyxd = true;
372  } else {
373  xd[k] = false;
374  ++k;
375  }
376  if (k > BXM1) {
377  break;
378  }
379  }
380 
381  // next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx
382 
383  fypix = -1;
384  nypix = 0;
385  lypix = 0;
386  logypx = 0;
387  for (i = 0; i < BYSIZE; ++i) {
388  if (ysum[i] > 0.) {
389  if (fypix == -1) {
390  fypix = i;
391  }
392  if (!yd[i]) {
393  ysort[logypx] = ysum[i];
394  ++logypx;
395  }
396  ++nypix;
397  lypix = i;
398  }
399  }
400 
401  // Make sure cluster is continuous
402 
403  if ((lypix - fypix + 1) != nypix || nypix == 0) {
404  LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster doesn't agree with number of pixels above threshold"
405  << ENDL;
406  if (theVerboseLevel > 2) {
407  LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
408  for (i = 0; i < BYSIZE - 1; ++i) {
409  LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";
410  }
411  LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE - 1] << ENDL;
412  }
413 
414  return 1;
415  }
416 
417  // If cluster is longer than max template size, technique fails
418 
419  if (nypix > TYSIZE) {
420  LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster is larger than maximum template size" << ENDL;
421  if (theVerboseLevel > 2) {
422  LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
423  for (i = 0; i < BYSIZE - 1; ++i) {
424  LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";
425  }
426  LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE - 1] << ENDL;
427  }
428 
429  return 6;
430  }
431 
432  // next, center the cluster on template center if necessary
433 
434  midpix = (fypix + lypix) / 2;
435  // shifty = BHY - midpix;
436  shifty = templ.cytemp() - midpix;
437  if (shifty > 0) {
438  for (i = lypix; i >= fypix; --i) {
439  ysum[i + shifty] = ysum[i];
440  ysum[i] = 0.;
441  yd[i + shifty] = yd[i];
442  yd[i] = false;
443  }
444  } else if (shifty < 0) {
445  for (i = fypix; i <= lypix; ++i) {
446  ysum[i + shifty] = ysum[i];
447  ysum[i] = 0.;
448  yd[i + shifty] = yd[i];
449  yd[i] = false;
450  }
451  }
452  lypix += shifty;
453  fypix += shifty;
454 
455  // If the cluster boundaries are OK, add pesudopixels, otherwise quit
456 
457  if (fypix > 1 && fypix < BYM2) {
458  ysum[fypix - 1] = pseudopix;
459  ysum[fypix - 2] = 0.2f * pseudopix;
460  } else {
461  return 8;
462  }
463  if (lypix > 1 && lypix < BYM2) {
464  ysum[lypix + 1] = pseudopix;
465  ysum[lypix + 2] = 0.2f * pseudopix;
466  } else {
467  return 8;
468  }
469 
470  // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
471 
472  if (ydouble[0]) {
473  originy = -0.5f;
474  } else {
475  originy = 0.f;
476  }
477 
478  // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx
479 
480  fxpix = -1;
481  nxpix = 0;
482  lxpix = 0;
483  logxpx = 0;
484  for (i = 0; i < BXSIZE; ++i) {
485  if (xsum[i] > 0.) {
486  if (fxpix == -1) {
487  fxpix = i;
488  }
489  if (!xd[i]) {
490  xsort[logxpx] = xsum[i];
491  ++logxpx;
492  }
493  ++nxpix;
494  lxpix = i;
495  }
496  }
497 
498  // Make sure cluster is continuous
499 
500  if ((lxpix - fxpix + 1) != nxpix) {
501  LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold"
502  << ENDL;
503  if (theVerboseLevel > 2) {
504  LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
505  for (i = 0; i < BXSIZE - 1; ++i) {
506  LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";
507  }
508  LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE - 1] << ENDL;
509  }
510 
511  return 2;
512  }
513 
514  // If cluster is longer than max template size, technique fails
515 
516  if (nxpix > TXSIZE) {
517  LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
518  if (theVerboseLevel > 2) {
519  LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
520  for (i = 0; i < BXSIZE - 1; ++i) {
521  LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";
522  }
523  LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE - 1] << ENDL;
524  }
525 
526  return 7;
527  }
528 
529  // next, center the cluster on template center if necessary
530 
531  midpix = (fxpix + lxpix) / 2;
532  // shiftx = BHX - midpix;
533  shiftx = templ.cxtemp() - midpix;
534  if (shiftx > 0) {
535  for (i = lxpix; i >= fxpix; --i) {
536  xsum[i + shiftx] = xsum[i];
537  xsum[i] = 0.f;
538  xd[i + shiftx] = xd[i];
539  xd[i] = false;
540  }
541  } else if (shiftx < 0) {
542  for (i = fxpix; i <= lxpix; ++i) {
543  xsum[i + shiftx] = xsum[i];
544  xsum[i] = 0.f;
545  xd[i + shiftx] = xd[i];
546  xd[i] = false;
547  }
548  }
549  lxpix += shiftx;
550  fxpix += shiftx;
551 
552  // If the cluster boundaries are OK, add pesudopixels, otherwise quit
553 
554  if (fxpix > 1 && fxpix < BXM2) {
555  xsum[fxpix - 1] = pseudopix;
556  xsum[fxpix - 2] = 0.2f * pseudopix;
557  } else {
558  return 9;
559  }
560  if (lxpix > 1 && lxpix < BXM2) {
561  xsum[lxpix + 1] = pseudopix;
562  xsum[lxpix + 2] = 0.2f * pseudopix;
563  } else {
564  return 9;
565  }
566 
567  // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
568 
569  if (xdouble[0]) {
570  originx = -0.5f;
571  } else {
572  originx = 0.f;
573  }
574 
575  // uncertainty and final corrections depend upon total charge bin
576 
577  qavg = templ.qavg();
578  fq = qtotal / qavg;
579  if (fq > 3.0f) {
580  binq = 0;
581  } else {
582  if (fq > 2.0f) {
583  binq = 1;
584  } else {
585  if (fq > 1.70f) {
586  binq = 2;
587  } else {
588  binq = 3;
589  }
590  }
591  }
592 
593  // Calculate the Vavilov probability that the cluster charge is consistent with a merged cluster
594 
595  if (speed < 0) {
596  templ.vavilov2_pars(mpv, sigmaQ, kappa);
597 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
598  if ((sigmaQ <= 0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
599  throw cms::Exception("DataCorrupt") << "SiPixelTemplateSplit::Vavilov parameters mpv/sigmaQ/kappa = " << mpv
600  << "/" << sigmaQ << "/" << kappa << std::endl;
601  }
602 #else
603  assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
604 #endif
605  xvav = ((double)qtotal - mpv) / sigmaQ;
606  // VVIObj is a private port of CERNLIB VVIDIS
607  VVIObj vvidist(kappa);
608  prvav = vvidist.fcn(xvav);
609  prob2Q = 1. - prvav;
610  if (prob2Q < prob2Qmin) {
611  prob2Q = prob2Qmin;
612  }
613  } else {
614  prob2Q = -1.f;
615  }
616 
617  // Return the charge bin via the parameter list unless the charge is too small (then flag it)
618 
619  q2bin = binq;
620  if (!deadpix && qtotal < 1.9f * templ.qmin()) {
621  q2bin = 5;
622  } else {
623  if (!deadpix && qtotal < 1.9f * templ.qmin(1)) {
624  q2bin = 4;
625  }
626  }
627 
628  if (theVerboseLevel > 9) {
629  LOGDEBUG("SiPixelTemplateSplit") << "ID = " << id << " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta
630  << " nclusx = " << nclusx << " nclusy = " << nclusy << ENDL;
631  }
632 
633  // Next, generate the 3d y- and x-templates
634 
635  templ.ytemp3d_int(nypix, nybin);
636 
637  ycbin = nybin / 2;
638 
639  templ.xtemp3d_int(nxpix, nxbin);
640 
641  // retrieve the number of x-bins
642 
643  xcbin = nxbin / 2;
644 
645  // First, decide on chi^2 min search parameters
646 
647 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
648  if (speed < -1 || speed > 2) {
649  throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::PixelTempReco2D called with illegal speed = " << speed
650  << std::endl;
651  }
652 #else
653  assert(speed >= -1 && speed < 3);
654 #endif
655  fybin = 0;
656  lybin = nybin - 1;
657  fxbin = 0;
658  lxbin = nxbin - 1;
659  djy = 1;
660  djx = 1;
661  if (speed > 0) {
662  djy = 2;
663  djx = 2;
664  if (speed > 1) {
665  if (!anyyd) {
666  djy = 4;
667  }
668  if (!anyxd) {
669  djx = 4;
670  }
671  }
672  }
673 
674  if (theVerboseLevel > 9) {
675  LOGDEBUG("SiPixelTemplateReco") << "fypix " << fypix << " lypix = " << lypix << " fybin = " << fybin
676  << " lybin = " << lybin << " djy = " << djy << " logypx = " << logypx << ENDL;
677  LOGDEBUG("SiPixelTemplateReco") << "fxpix " << fxpix << " lxpix = " << lxpix << " fxbin = " << fxbin
678  << " lxbin = " << lxbin << " djx = " << djx << " logxpx = " << logxpx << ENDL;
679  }
680 
681  // Do the y-reconstruction first
682 
683  // Define the maximum signal to allow before de-weighting a pixel
684 
685  sythr = 2.1f * (templ.symax());
686 
687  // Make sure that there will be at least two pixels that are not de-weighted
688 
689  std::sort(&ysort[0], &ysort[logypx]);
690  if (logypx == 1) {
691  sythr = 1.01f * ysort[0];
692  } else {
693  if (ysort[1] > sythr) {
694  sythr = 1.01f * ysort[1];
695  }
696  }
697 
698  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
699 
700  // for(i=0; i<BYSIZE; ++i) { ysig2[i] = 0.;}
701  templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
702 
703  // Find the template bin that minimizes the Chi^2
704 
705  chi2ymin = 1.e15;
706  ss2 = 0.f;
707  for (i = fypix - 2; i <= lypix + 2; ++i) {
708  yw2[i] = 1.f / ysig2[i];
709  ysw[i] = ysum[i] * yw2[i];
710  ss2 += ysum[i] * ysw[i];
711  }
712  minbinj = -1;
713  minbink = -1;
714  deltaj = djy;
715  jmin = fybin;
716  jmax = lybin;
717  kmin = fybin;
718  kmax = lybin;
719  std::vector<float> ytemp(BYSIZE);
720  while (deltaj > 0) {
721  for (j = jmin; j < jmax; j += deltaj) {
722  km = std::min(kmax, j);
723  for (k = kmin; k <= km; k += deltaj) {
724  // Get the template for this set of indices
725 
726  templ.ytemp3d(j, k, ytemp);
727 
728  // Modify the template if double pixels are present
729 
730  if (nypix > logypx) {
731  i = fypix;
732  while (i < lypix) {
733  if (yd[i] && !yd[i + 1]) {
734  // Sum the adjacent cells and put the average signal in both
735 
736  sigavg = (ytemp[i] + ytemp[i + 1]) / 2.f;
737  ytemp[i] = sigavg;
738  ytemp[i + 1] = sigavg;
739  i += 2;
740  } else {
741  ++i;
742  }
743  }
744  }
745  ssa = 0.f;
746  sa2 = 0.f;
747  for (i = fypix - 2; i <= lypix + 2; ++i) {
748  ssa += ysw[i] * ytemp[i];
749  sa2 += ytemp[i] * ytemp[i] * yw2[i];
750  }
751  rat = ssa / ss2;
752  if (rat <= 0.) {
753  LOGERROR("SiPixelTemplateSplit") << "illegal chi2ymin normalization = " << rat << ENDL;
754  rat = 1.;
755  }
756  chi2y = ss2 - 2.f * ssa / rat + sa2 / (rat * rat);
757  if (chi2y < chi2ymin) {
758  chi2ymin = chi2y;
759  minbinj = j;
760  minbink = k;
761  }
762  }
763  }
764  deltaj /= 2;
765  if (minbinj > fybin) {
766  jmin = minbinj - deltaj;
767  } else {
768  jmin = fybin;
769  }
770  if (minbinj < lybin) {
771  jmax = minbinj + deltaj;
772  } else {
773  jmax = lybin;
774  }
775  if (minbink > fybin) {
776  kmin = minbink - deltaj;
777  } else {
778  kmin = fybin;
779  }
780  if (minbink < lybin) {
781  kmax = minbink + deltaj;
782  } else {
783  kmax = lybin;
784  }
785  }
786 
787  if (theVerboseLevel > 9) {
788  LOGDEBUG("SiPixelTemplateReco") << "minbins " << minbinj << "," << minbink << " chi2ymin = " << chi2ymin << ENDL;
789  }
790 
791  // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
792 
793  if (logypx == 1) {
794  if (nypix == 1) {
795  delta = templ.dyone();
796  sigma = templ.syone();
797  } else {
798  delta = templ.dytwo();
799  sigma = templ.sytwo();
800  }
801 
802  yrec1 = 0.5f * (fypix + lypix - 2 * shifty + 2.f * originy) * ysize - delta;
803  yrec2 = yrec1;
804 
805  if (sigma <= 0.) {
806  sigmay = ysize / sqrt12;
807  } else {
808  sigmay = sigma;
809  }
810 
811  // Do probability calculation for one-pixel clusters
812 
813  chi21max = fmax(chi21min, (double)templ.chi2yminone());
814  chi2ymin -= chi21max;
815  if (chi2ymin < 0.) {
816  chi2ymin = 0.;
817  }
818  // prob2y = gsl_cdf_chisq_Q(chi2ymin, mean1pix);
819  meany = fmax(mean1pix, (double)templ.chi2yavgone());
820  hchi2 = chi2ymin / 2.;
821  hndof = meany / 2.;
822  prob2y = 1. - TMath::Gamma(hndof, hchi2);
823 
824  } else {
825  // For cluster > 1 pix, use chi^2 minimm to recontruct the two y-positions
826 
827  // at small eta, the templates won't actually work on two pixel y-clusters so just return the pixel centers
828 
829  if (logypx == 2 && fabsf(cotbeta) < 0.25f) {
830  switch (nypix) {
831  case 2:
832  // Both pixels are small
833  yrec1 = (fypix - shifty + originy) * ysize;
834  yrec2 = (lypix - shifty + originy) * ysize;
835  sigmay = ysize / sqrt12;
836  break;
837  case 3:
838  // One big pixel and one small pixel
839  if (yd[fypix]) {
840  yrec1 = (fypix + 0.5f - shifty + originy) * ysize;
841  yrec2 = (lypix - shifty + originy) * ysize;
842  sigmay = ysize / sqrt12;
843  } else {
844  yrec1 = (fypix - shifty + originy) * ysize;
845  yrec2 = (lypix - 0.5f - shifty + originy) * ysize;
846  sigmay = 1.5f * ysize / sqrt12;
847  }
848  break;
849  case 4:
850  // Two big pixels
851  yrec1 = (fypix + 0.5f - shifty + originy) * ysize;
852  yrec2 = (lypix - 0.5f - shifty + originy) * ysize;
853  sigmay = 2.f * ysize / sqrt12;
854  break;
855  default:
856  // Something is screwy ...
857  LOGERROR("SiPixelTemplateReco")
858  << "weird problem: logical y-pixels = " << logypx << ", total ysize in normal pixels = " << nypix << ENDL;
859  return 10;
860  }
861  } else {
862  // uncertainty and final correction depend upon charge bin
863 
864  bias = templ.yavgc2m(binq);
865  yrec1 = (0.125f * (minbink - ycbin) + BHY - (float)shifty + originy) * ysize - bias;
866  yrec2 = (0.125f * (minbinj - ycbin) + BHY - (float)shifty + originy) * ysize - bias;
867  sigmay = sqrt2y * templ.yrmsc2m(binq);
868  }
869 
870  // Do goodness of fit test in y
871 
872  if (chi2ymin < 0.0) {
873  chi2ymin = 0.0;
874  }
875  meany = templ.chi2yavgc2m(binq);
876  if (meany < 0.01) {
877  meany = 0.01;
878  }
879  // gsl function that calculates the chi^2 tail prob for non-integral dof
880  // prob2y = gsl_cdf_chisq_Q(chi2y, meany);
881  // prob2y = ROOT::Math::chisquared_cdf_c(chi2y, meany);
882  hchi2 = chi2ymin / 2.;
883  hndof = meany / 2.;
884  prob2y = 1. - TMath::Gamma(hndof, hchi2);
885  }
886 
887  // Do the x-reconstruction next
888 
889  // Define the maximum signal to allow before de-weighting a pixel
890 
891  sxthr = 2.1f * templ.sxmax();
892 
893  // Make sure that there will be at least two pixels that are not de-weighted
894  std::sort(&xsort[0], &xsort[logxpx]);
895  if (logxpx == 1) {
896  sxthr = 1.01f * xsort[0];
897  } else {
898  if (xsort[1] > sxthr) {
899  sxthr = 1.01f * xsort[1];
900  }
901  }
902 
903  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
904 
905  // for(i=0; i<BYSIZE; ++i) { xsig2[i] = 0.; }
906  templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
907 
908  // Find the template bin that minimizes the Chi^2
909 
910  chi2xmin = 1.e15;
911  ss2 = 0.f;
912  for (i = fxpix - 2; i <= lxpix + 2; ++i) {
913  xw2[i] = 1.f / xsig2[i];
914  xsw[i] = xsum[i] * xw2[i];
915  ss2 += xsum[i] * xsw[i];
916  }
917  minbinj = -1;
918  minbink = -1;
919  deltaj = djx;
920  jmin = fxbin;
921  jmax = lxbin;
922  kmin = fxbin;
923  kmax = lxbin;
924  std::vector<float> xtemp(BXSIZE);
925  while (deltaj > 0) {
926  for (j = jmin; j < jmax; j += deltaj) {
927  km = std::min(kmax, j);
928  for (k = kmin; k <= km; k += deltaj) {
929  // Get the template for this set of indices
930 
931  templ.xtemp3d(j, k, xtemp);
932 
933  // Modify the template if double pixels are present
934 
935  if (nxpix > logxpx) {
936  i = fxpix;
937  while (i < lxpix) {
938  if (xd[i] && !xd[i + 1]) {
939  // Sum the adjacent cells and put the average signal in both
940 
941  sigavg = (xtemp[i] + xtemp[i + 1]) / 2.f;
942  xtemp[i] = sigavg;
943  xtemp[i + 1] = sigavg;
944  i += 2;
945  } else {
946  ++i;
947  }
948  }
949  }
950  ssa = 0.f;
951  sa2 = 0.f;
952  for (i = fxpix - 2; i <= lxpix + 2; ++i) {
953  ssa += xsw[i] * xtemp[i];
954  sa2 += xtemp[i] * xtemp[i] * xw2[i];
955  }
956  rat = ssa / ss2;
957  if (rat <= 0.f) {
958  LOGERROR("SiPixelTemplateSplit") << "illegal chi2xmin normalization = " << rat << ENDL;
959  rat = 1.;
960  }
961  chi2x = ss2 - 2.f * ssa / rat + sa2 / (rat * rat);
962  if (chi2x < chi2xmin) {
963  chi2xmin = chi2x;
964  minbinj = j;
965  minbink = k;
966  }
967  }
968  }
969  deltaj /= 2;
970  if (minbinj > fxbin) {
971  jmin = minbinj - deltaj;
972  } else {
973  jmin = fxbin;
974  }
975  if (minbinj < lxbin) {
976  jmax = minbinj + deltaj;
977  } else {
978  jmax = lxbin;
979  }
980  if (minbink > fxbin) {
981  kmin = minbink - deltaj;
982  } else {
983  kmin = fxbin;
984  }
985  if (minbink < lxbin) {
986  kmax = minbink + deltaj;
987  } else {
988  kmax = lxbin;
989  }
990  }
991 
992  if (theVerboseLevel > 9) {
993  LOGDEBUG("SiPixelTemplateSplit") << "minbinj/minbink " << minbinj << "/" << minbink << " chi2xmin = " << chi2xmin
994  << ENDL;
995  }
996 
997  // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
998 
999  if (logxpx == 1) {
1000  if (nxpix == 1) {
1001  delta = templ.dxone();
1002  sigma = templ.sxone();
1003  } else {
1004  delta = templ.dxtwo();
1005  sigma = templ.sxtwo();
1006  }
1007  xrec1 = 0.5f * (fxpix + lxpix - 2 * shiftx + 2.f * originx) * xsize - delta;
1008  xrec2 = xrec1;
1009  if (sigma <= 0.f) {
1010  sigmax = xsize / sqrt12;
1011  } else {
1012  sigmax = sigma;
1013  }
1014 
1015  // Do probability calculation for one-pixel clusters
1016 
1017  chi21max = fmax(chi21min, (double)templ.chi2xminone());
1018  chi2xmin -= chi21max;
1019  if (chi2xmin < 0.) {
1020  chi2xmin = 0.;
1021  }
1022  meanx = fmax(mean1pix, (double)templ.chi2xavgone());
1023  hchi2 = chi2xmin / 2.;
1024  hndof = meanx / 2.;
1025  prob2x = 1. - TMath::Gamma(hndof, hchi2);
1026 
1027  } else {
1028  // For cluster > 1 pix, use chi^2 minimm to recontruct the two x-positions
1029 
1030  // uncertainty and final correction depend upon charge bin
1031 
1032  bias = templ.xavgc2m(binq);
1033  xrec1 = (0.125f * (minbink - xcbin) + BHX - (float)shiftx + originx) * xsize - bias;
1034  xrec2 = (0.125f * (minbinj - xcbin) + BHX - (float)shiftx + originx) * xsize - bias;
1035  sigmax = sqrt2x * templ.xrmsc2m(binq);
1036 
1037  // Do goodness of fit test in y
1038 
1039  if (chi2xmin < 0.0) {
1040  chi2xmin = 0.0;
1041  }
1042  meanx = templ.chi2xavgc2m(binq);
1043  if (meanx < 0.01) {
1044  meanx = 0.01;
1045  }
1046  hchi2 = chi2xmin / 2.;
1047  hndof = meanx / 2.;
1048  prob2x = 1. - TMath::Gamma(hndof, hchi2);
1049  }
1050 
1051  // Don't return exact zeros for the probability
1052 
1053  if (prob2y < probmin) {
1054  prob2y = probmin;
1055  }
1056  if (prob2x < probmin) {
1057  prob2x = probmin;
1058  }
1059 
1060  // New code: resolve the ambiguity if resolve is set to true
1061 
1062  dchisq = 0.;
1063  if (resolve) {
1064  // First copy the unexpanded cluster into a new working buffer
1065 
1066  for (i = 0; i < BYM2; ++i) {
1067  for (j = 0; j < BXM2; ++j) {
1068  if ((i > 0 && i < BYM3) && (j > 0 && j < BXM3)) {
1069  cluster2[j][i] = qscale * clust[j - 1][i - 1];
1070  } else {
1071  cluster2[j][i] = 0.f;
1072  }
1073  }
1074  }
1075 
1076  // Next, redefine the local coordinates to start at the lower LH corner of pixel[1][1];
1077 
1078  if (xdouble[0]) {
1079  x1p = xrec1 + xsize;
1080  x2p = xrec2 + xsize;
1081  } else {
1082  x1p = xrec1 + xsize / 2.f;
1083  x2p = xrec2 + xsize / 2.f;
1084  }
1085 
1086  if (ydouble[0]) {
1087  y1p = yrec1 + ysize;
1088  y2p = yrec2 + ysize;
1089  } else {
1090  y1p = yrec1 + ysize / 2.f;
1091  y2p = yrec2 + ysize / 2.f;
1092  }
1093 
1094  // Next, calculate 2-d templates for the (x1,y1)+(x2,y2) and the (x1,y2)+(x2,y1) hypotheses
1095 
1096  // First zero the two 2-d hypotheses
1097 
1098  for (i = 0; i < BYM2; ++i) {
1099  for (j = 0; j < BXM2; ++j) {
1100  temp2d1[j][i] = 0.f;
1101  temp2d2[j][i] = 0.f;
1102  }
1103  }
1104 
1105  // Add the two hits in the first hypothesis
1106 
1107  any2dfail = templ2D.xytemp(id, cotalpha, cotbeta, x1p, y1p, xdouble, ydouble, temp2d1) &&
1108  templ2D.xytemp(id, cotalpha, cotbeta, x2p, y2p, xdouble, ydouble, temp2d1);
1109 
1110  // And then the second hypothesis
1111 
1112  any2dfail = any2dfail && templ2D.xytemp(id, cotalpha, cotbeta, x1p, y2p, xdouble, ydouble, temp2d2);
1113 
1114  any2dfail = any2dfail && templ2D.xytemp(id, cotalpha, cotbeta, x2p, y1p, xdouble, ydouble, temp2d2);
1115 
1116  // If any of these have failed, use the simple templates instead
1117 
1118  if (!any2dfail) {
1119  // Rezero the two 2-d hypotheses
1120 
1121  for (i = 0; i < BYM2; ++i) {
1122  for (j = 0; j < BXM2; ++j) {
1123  temp2d1[j][i] = 0.f;
1124  temp2d2[j][i] = 0.f;
1125  }
1126  }
1127 
1128  // Add the two hits in the first hypothesis
1129 
1130  if (!templ.simpletemplate2D(x1p, y1p, xdouble, ydouble, temp2d1)) {
1131  return 1;
1132  }
1133 
1134  if (!templ.simpletemplate2D(x2p, y2p, xdouble, ydouble, temp2d1)) {
1135  return 1;
1136  }
1137 
1138  // And then the second hypothesis
1139 
1140  if (!templ.simpletemplate2D(x1p, y2p, xdouble, ydouble, temp2d2)) {
1141  return 1;
1142  }
1143 
1144  if (!templ.simpletemplate2D(x2p, y1p, xdouble, ydouble, temp2d2)) {
1145  return 1;
1146  }
1147  }
1148 
1149  // Keep lists of pixels and nearest neighbors
1150 
1151  std::list<std::pair<int, int> > pixellst;
1152 
1153  // Loop through the array and find non-zero elements
1154 
1155  for (i = 0; i < BYM2; ++i) {
1156  for (j = 0; j < BXM2; ++j) {
1157  if (cluster2[j][i] > 0.f || temp2d1[j][i] > 0.f || temp2d2[j][i] > 0.f) {
1158  pixel.first = j;
1159  pixel.second = i;
1160  pixellst.push_back(pixel);
1161  }
1162  }
1163  }
1164 
1165  // Now calculate the product of Landau probabilities (alpha probability)
1166 
1167  templ2D.landau_par(lanpar);
1168  loglike1 = 0.;
1169  loglike2 = 0.;
1170 
1171  // Now, for each neighbor, match it again the pixel list.
1172  // If found, delete it from the neighbor list
1173 
1174  std::list<std::pair<int, int> >::const_iterator pixIter, pixEnd;
1175  pixIter = pixellst.begin();
1176  pixEnd = pixellst.end();
1177  for (; pixIter != pixEnd; ++pixIter) {
1178  j = pixIter->first;
1179  i = pixIter->second;
1180  if ((i < BHY && cotbeta > 0.) || (i >= BHY && cotbeta < 0.)) {
1181  lparm = 0;
1182  } else {
1183  lparm = 1;
1184  }
1185  mpv1 = lanpar[lparm][0] + lanpar[lparm][1] * temp2d1[j][i];
1186  sigmal1 = lanpar[lparm][2] * mpv1;
1187  sigmal1 = sqrt(sigmal1 * sigmal1 + lanpar[lparm][3] * lanpar[lparm][3]);
1188  if (sigmal1 < q05)
1189  sigmal1 = q05;
1190  arg1 = (cluster2[j][i] - mpv1) / sigmal1 - 0.22278;
1191  mpv2 = lanpar[lparm][0] + lanpar[lparm][1] * temp2d2[j][i];
1192  sigmal2 = lanpar[lparm][2] * mpv2;
1193  sigmal2 = sqrt(sigmal2 * sigmal2 + lanpar[lparm][3] * lanpar[lparm][3]);
1194  if (sigmal2 < q05)
1195  sigmal2 = q05;
1196  arg2 = (cluster2[j][i] - mpv2) / sigmal2 - 0.22278;
1197  // like = ROOT::Math::landau_pdf(arg1)/sigmal1;
1198  like = ROOT::Math::landau_pdf(arg1);
1199  if (like < 1.e-30)
1200  like = 1.e-30;
1201  loglike1 += log(like);
1202  // like = ROOT::Math::landau_pdf(arg2)/sigmal2;
1203  like = ROOT::Math::landau_pdf(arg2);
1204  if (like < 1.e-30)
1205  like = 1.e-30;
1206  loglike2 += log(like);
1207  }
1208 
1209  // Calculate chisquare difference for the two hypotheses 9don't multiply by 2 for less inconvenient scaling
1210 
1211  deltachi2 = loglike1 - loglike2;
1212 
1213  if (deltachi2 < 0.) {
1214  // Flip the x1 and x2
1215 
1216  x1p = xrec1;
1217  xrec1 = xrec2;
1218  xrec2 = x1p;
1219  }
1220 
1221  // Return a positive definite value
1222 
1223  dchisq = fabs(deltachi2);
1224  }
1225 
1226  return 0;
1227 } // PixelTempSplit
SiPixelTemplate::yrmsc2m
float yrmsc2m(int i)
1st pass chi2 min search: average y-rms of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:572
SiPixelTemplate::vavilov2_pars
void vavilov2_pars(double &mpv, double &sigma, double &kappa)
Definition: SiPixelTemplate.cc:4086
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:30
BXSIZE
#define BXSIZE
Definition: SiPixelTemplateDefs.h:33
SiPixelTemplate::chi2xminone
float chi2xminone()
//!< minimum of x chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:660
mps_fire.i
i
Definition: mps_fire.py:429
SiPixelTemplate::chi2yavgc2m
float chi2yavgc2m(int i)
1st pass chi2 min search: average y-chisq for merged clusters
Definition: SiPixelTemplate.h:582
SiPixelTemplate::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
Definition: SiPixelTemplate.cc:1328
SiPixelTemplate2D::xytemp
bool xytemp(float xhit, float yhit, bool ydouble[21+2], bool xdouble[13+2], float template2d[13+2][21+2], bool dervatives, float dpdx2d[2][13+2][21+2], float &QTemplate)
Definition: SiPixelTemplate2D.cc:992
SiPixelTemplate::symax
float symax()
average pixel signal for y-projection of cluster
Definition: SiPixelTemplate.h:412
SiPixelTemplate::cytemp
int cytemp()
Return central pixel of y template pixels above readout threshold.
Definition: SiPixelTemplate.cc:2656
BYSIZE
#define BYSIZE
Definition: SiPixelTemplateDefs.h:25
SiPixelTemplate::xsigma2
void xsigma2(int fxpix, int lxpix, float sxthr, float xsum[13+4], float xsig2[13+4])
Definition: SiPixelTemplate.cc:2274
TXSIZE
#define TXSIZE
Definition: SiPixelTemplateDefs.h:30
LOGERROR
#define LOGERROR(x)
Definition: SiPixelTemplateSplit.cc:49
BXM3
#define BXM3
Definition: SiPixelTemplateDefs.h:37
BXM1
#define BXM1
Definition: SiPixelTemplateDefs.h:35
SiPixelTemplate::simpletemplate2D
bool simpletemplate2D(float xhitp, float yhitp, std::vector< bool > &ydouble, std::vector< bool > &xdouble, float template2d[13+2][21+2])
Make simple 2-D templates from track angles set in interpolate and hit position.
Definition: SiPixelTemplate.cc:3770
SiPixelTemplate::qmin
float qmin()
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
Definition: SiPixelTemplate.h:422
SiPixelTemplate::xavgc2m
float xavgc2m(int i)
1st pass chi2 min search: average x-bias of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:604
maxpix
static const int maxpix
Definition: SiPixelLorentzAngle.h:48
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
ysize
const Int_t ysize
Definition: trackSplitPlot.h:43
SiPixelTemplate::sytwo
float sytwo()
rms for one double-pixel y-clusters
Definition: SiPixelTemplate.h:416
cms::cuda::assert
assert(be >=bs)
ENDL
#define ENDL
Definition: SiPixelTemplateSplit.cc:52
SiPixelTemplate::chi2yminone
float chi2yminone()
//!< minimum of y chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:658
SiPixelTemplate::sxone
float sxone()
rms for one pixel x-clusters
Definition: SiPixelTemplate.h:419
BYM1
#define BYM1
Definition: SiPixelTemplateDefs.h:27
BXM2
#define BXM2
Definition: SiPixelTemplateDefs.h:36
dumpMFGeometry_cfg.delta
delta
Definition: dumpMFGeometry_cfg.py:25
SiPixelTemplate::cxtemp
int cxtemp()
Return central pixel of x-template pixels above readout threshold.
Definition: SiPixelTemplate.cc:2709
SiPixelTemplate::qscale
float qscale()
charge scaling factor
Definition: SiPixelTemplate.h:409
SiPixelTemplate::chi2xavgc2m
float chi2xavgc2m(int i)
1st pass chi2 min search: average x-chisq for merged clusters
Definition: SiPixelTemplate.h:624
SiPixelTemplate::dxone
float dxone()
mean offset/correction for one pixel x-clusters
Definition: SiPixelTemplate.h:418
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
BHX
#define BHX
Definition: SiPixelTemplateDefs.h:34
SiPixelTemplate::ytemp3d
void ytemp3d(int j, int k, std::vector< float > &ytemplate)
Definition: SiPixelTemplate.cc:2842
LOGDEBUG
#define LOGDEBUG(x)
Definition: SiPixelTemplateSplit.cc:50
SiPixelTemplate::xsize
float xsize()
pixel x-size (microns)
Definition: SiPixelTemplate.h:675
SiPixelTemplate::ysigma2
void ysigma2(int fypix, int lypix, float sythr, float ysum[21+4], float ysig2[21+4])
BYM2
#define BYM2
Definition: SiPixelTemplateDefs.h:28
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
SiPixelTemplate::sxtwo
float sxtwo()
rms for one double-pixel x-clusters
Definition: SiPixelTemplate.h:421
BYM3
#define BYM3
Definition: SiPixelTemplateDefs.h:29
SiPixelTemplate::dytwo
float dytwo()
mean offset/correction for one double-pixel y-clusters
Definition: SiPixelTemplate.h:415
SiPixelTemplate::s50
float s50()
1/2 of the pixel threshold signal in electrons
Definition: SiPixelTemplate.h:410
HLTSiStripMonitoring_cff.speed
speed
Definition: HLTSiStripMonitoring_cff.py:22
SiPixelTemplate::yavgc2m
float yavgc2m(int i)
1st pass chi2 min search: average y-bias of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:562
SiPixelTemplate::xtemp3d
void xtemp3d(int j, int k, std::vector< float > &xtemplate)
Definition: SiPixelTemplate.cc:2944
TYSIZE
#define TYSIZE
Definition: SiPixelTemplateDefs.h:21
SiPixelTemplate::syone
float syone()
rms for one pixel y-clusters
Definition: SiPixelTemplate.h:414
SiPixelTemplate::ytemp3d_int
void ytemp3d_int(int nypix, int &nybins)
Definition: SiPixelTemplate.cc:2764
SiPixelTemplate::chi2yavgone
float chi2yavgone()
//!< average y chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:657
muonClassificationByHits_cfi.pixel
pixel
Definition: muonClassificationByHits_cfi.py:9
SiPixelTemplate::xtemp3d_int
void xtemp3d_int(int nxpix, int &nxbins)
Definition: SiPixelTemplate.cc:2865
SiPixelTemplate::qavg
float qavg()
average cluster charge for this set of track angles
Definition: SiPixelTemplate.h:407
BHY
#define BHY
Definition: SiPixelTemplateDefs.h:26
SiPixelTemplate::sxmax
float sxmax()
average pixel signal for x-projection of cluster
Definition: SiPixelTemplate.h:417
VVIObj
Definition: VVIObj.h:24
SiPixelTemplate::chi2xavgone
float chi2xavgone()
//!< average x chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:659
SiPixelTemplate::pixmax
float pixmax()
maximum pixel charge
Definition: SiPixelTemplate.h:408
SiPixelTemplate2D::landau_par
void landau_par(float lanpar[2][5])
Return the Landau probability parameters for this set of cot(alpha, cot(beta)
Definition: SiPixelTemplate2D.cc:1761
hgcalLayerClusters_cfi.kappa
kappa
Definition: hgcalLayerClusters_cfi.py:35
theVerboseLevel
static constexpr int theVerboseLevel
Definition: SiPixelTemplateSplit.cc:51
xsize
const Int_t xsize
Definition: trackSplitPlot.h:42
SiPixelTemplate::dyone
float dyone()
mean offset/correction for one pixel y-clusters
Definition: SiPixelTemplate.h:413
SiPixelTemplate::dxtwo
float dxtwo()
mean offset/correction for one double-pixel x-clusters
Definition: SiPixelTemplate.h:420
SiPixelTemplateSplit::array_2d
boost::multi_array< float, 2 > array_2d
Definition: SiPixelTemplateSplit.h:47
SiPixelTemplate::ysize
float ysize()
pixel y-size (microns)
Definition: SiPixelTemplate.h:676
SiPixelTemplate::xrmsc2m
float xrmsc2m(int i)
1st pass chi2 min search: average x-rms of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:614

◆ PixelTempSplit() [2/5]

int SiPixelTemplateSplit::PixelTempSplit ( int  id,
float  cotalpha,
float  cotbeta,
array_2d &  cluster,
std::vector< bool > &  ydouble,
std::vector< bool > &  xdouble,
SiPixelTemplate &  templ,
float &  yrec1,
float &  yrec2,
float &  sigmay,
float &  prob2y,
float &  xrec1,
float &  xrec2,
float &  sigmax,
float &  prob2x,
int &  q2bin,
float &  prob2Q,
bool  resolve,
int  speed,
float &  dchisq,
SiPixelTemplate2D &  templ2D 
)

Reconstruct the best estimate of the hit positions 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
yrec1- (output) best estimate of first y-coordinate of hit in microns
yrec2- (output) best estimate of second y-coordinate of hit in microns
sigmay- (output) best estimate of uncertainty on yrec in microns
prob2y- (output) probability describing goodness-of-fit for y-reco
xrec1- (output) best estimate of first x-coordinate of hit in microns
xrec2- (output) best estimate of second x-coordinate of hit in microns
sigmax- (output) best estimate of uncertainty on xrec in microns
prob2x- (output) probability describing goodness-of-fit for x-reco
q2bin- (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]
prob2Q- (output) probability that the cluster charge is compatible with a 2-hit merging
resolve- (input) use 2-D chisquare to resolve the 2-fold x-y association ambiguity (costs a factor of two in speed)
dchisq- (output) the delta chisquare estimator used to break the degeneracy (0 means no discrimination, larger than 0.1 is good)

Definition at line 1254 of file SiPixelTemplateSplit.cc.

References PixelTempSplit(), and HLTSiStripMonitoring_cff::speed.

1274  {
1275  // Local variables
1276  const bool deadpix = false;
1277  std::vector<std::pair<int, int> > zeropix;
1278 
1279  return SiPixelTemplateSplit::PixelTempSplit(id,
1280  cotalpha,
1281  cotbeta,
1282  cluster,
1283  ydouble,
1284  xdouble,
1285  templ,
1286  yrec1,
1287  yrec2,
1288  sigmay,
1289  prob2y,
1290  xrec1,
1291  xrec2,
1292  sigmax,
1293  prob2x,
1294  q2bin,
1295  prob2Q,
1296  resolve,
1297  speed,
1298  dchisq,
1299  deadpix,
1300  zeropix,
1301  templ2D);
1302 
1303 } // PixelTempSplit
SiPixelTemplateSplit::PixelTempSplit
int PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec1, float &yrec2, float &sigmay, float &prob2y, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, float &prob2Q, bool resolve, int speed, float &dchisq, bool deadpix, std::vector< std::pair< int, int > > &zeropix, SiPixelTemplate2D &templ2D)
Definition: SiPixelTemplateSplit.cc:98
HLTSiStripMonitoring_cff.speed
speed
Definition: HLTSiStripMonitoring_cff.py:22

◆ PixelTempSplit() [3/5]

int SiPixelTemplateSplit::PixelTempSplit ( int  id,
float  cotalpha,
float  cotbeta,
array_2d &  cluster,
std::vector< bool > &  ydouble,
std::vector< bool > &  xdouble,
SiPixelTemplate &  templ,
float &  yrec1,
float &  yrec2,
float &  sigmay,
float &  prob2y,
float &  xrec1,
float &  xrec2,
float &  sigmax,
float &  prob2x,
int &  q2bin,
float &  prob2Q,
bool  resolve,
float &  dchisq,
SiPixelTemplate2D &  templ2D 
)

Reconstruct the best estimate of the hit positions 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
yrec1- (output) best estimate of first y-coordinate of hit in microns
yrec2- (output) best estimate of second y-coordinate of hit in microns
sigmay- (output) best estimate of uncertainty on yrec in microns
prob2y- (output) probability describing goodness-of-fit for y-reco
xrec1- (output) best estimate of first x-coordinate of hit in microns
xrec2- (output) best estimate of second x-coordinate of hit in microns
sigmax- (output) best estimate of uncertainty on xrec in microns
prob2x- (output) probability describing goodness-of-fit for x-reco
q2bin- (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]
prob2Q- (output) probability that the cluster charge is compatible with a 2-hit merging
resolve- (input) use 2-D chisquare to resolve the 2-fold x-y association ambiguity (costs a factor of two in speed)
dchisq- (output) the delta chisquare estimator used to break the degeneracy (0 means no discrimination, larger than 0.1 is good)

Definition at line 1330 of file SiPixelTemplateSplit.cc.

References PixelTempSplit(), and HLTSiStripMonitoring_cff::speed.

1349  {
1350  // Local variables
1351  const bool deadpix = false;
1352  std::vector<std::pair<int, int> > zeropix;
1353  const int speed = 1;
1354 
1355  return SiPixelTemplateSplit::PixelTempSplit(id,
1356  cotalpha,
1357  cotbeta,
1358  cluster,
1359  ydouble,
1360  xdouble,
1361  templ,
1362  yrec1,
1363  yrec2,
1364  sigmay,
1365  prob2y,
1366  xrec1,
1367  xrec2,
1368  sigmax,
1369  prob2x,
1370  q2bin,
1371  prob2Q,
1372  resolve,
1373  speed,
1374  dchisq,
1375  deadpix,
1376  zeropix,
1377  templ2D);
1378 
1379 } // PixelTempSplit
SiPixelTemplateSplit::PixelTempSplit
int PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec1, float &yrec2, float &sigmay, float &prob2y, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, float &prob2Q, bool resolve, int speed, float &dchisq, bool deadpix, std::vector< std::pair< int, int > > &zeropix, SiPixelTemplate2D &templ2D)
Definition: SiPixelTemplateSplit.cc:98
HLTSiStripMonitoring_cff.speed
speed
Definition: HLTSiStripMonitoring_cff.py:22

◆ PixelTempSplit() [4/5]

int SiPixelTemplateSplit::PixelTempSplit ( int  id,
float  cotalpha,
float  cotbeta,
array_2d &  cluster,
std::vector< bool > &  ydouble,
std::vector< bool > &  xdouble,
SiPixelTemplate &  templ,
float &  yrec1,
float &  yrec2,
float &  sigmay,
float &  prob2y,
float &  xrec1,
float &  xrec2,
float &  sigmax,
float &  prob2x,
int &  q2bin,
float &  prob2Q,
SiPixelTemplate2D &  templ2D 
)

Reconstruct the best estimate of the hit positions for pixel clusters. Legacy interface to older code.

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
yrec1- (output) best estimate of first y-coordinate of hit in microns
yrec2- (output) best estimate of second y-coordinate of hit in microns
sigmay- (output) best estimate of uncertainty on yrec in microns
prob2y- (output) probability describing goodness-of-fit for y-reco
xrec1- (output) best estimate of first x-coordinate of hit in microns
xrec2- (output) best estimate of second x-coordinate of hit in microns
sigmax- (output) best estimate of uncertainty on xrec in microns
prob2x- (output) probability describing goodness-of-fit for x-reco
q2bin- (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]
prob2Q- (output) probability that the cluster charge is compatible with a 2-hit merging

Definition at line 1404 of file SiPixelTemplateSplit.cc.

References PixelTempSplit(), and HLTSiStripMonitoring_cff::speed.

1421  {
1422  // Local variables
1423  const bool deadpix = false;
1424  const bool resolve = true;
1425  float dchisq;
1426  std::vector<std::pair<int, int> > zeropix;
1427  const int speed = 1;
1428 
1429  return SiPixelTemplateSplit::PixelTempSplit(id,
1430  cotalpha,
1431  cotbeta,
1432  cluster,
1433  ydouble,
1434  xdouble,
1435  templ,
1436  yrec1,
1437  yrec2,
1438  sigmay,
1439  prob2y,
1440  xrec1,
1441  xrec2,
1442  sigmax,
1443  prob2x,
1444  q2bin,
1445  prob2Q,
1446  resolve,
1447  speed,
1448  dchisq,
1449  deadpix,
1450  zeropix,
1451  templ2D);
1452 
1453 } // PixelTempSplit
SiPixelTemplateSplit::PixelTempSplit
int PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec1, float &yrec2, float &sigmay, float &prob2y, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, float &prob2Q, bool resolve, int speed, float &dchisq, bool deadpix, std::vector< std::pair< int, int > > &zeropix, SiPixelTemplate2D &templ2D)
Definition: SiPixelTemplateSplit.cc:98
HLTSiStripMonitoring_cff.speed
speed
Definition: HLTSiStripMonitoring_cff.py:22

◆ PixelTempSplit() [5/5]

int SiPixelTemplateSplit::PixelTempSplit ( int  id,
float  cotalpha,
float  cotbeta,
array_2d &  cluster,
std::vector< bool > &  ydouble,
std::vector< bool > &  xdouble,
SiPixelTemplate &  templ,
float &  yrec1,
float &  yrec2,
float &  sigmay,
float &  prob2y,
float &  xrec1,
float &  xrec2,
float &  sigmax,
float &  prob2x,
int &  q2bin,
SiPixelTemplate2D &  templ2D 
)

Reconstruct the best estimate of the hit positions for pixel clusters. Legacy interface to older code.

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
yrec1- (output) best estimate of first y-coordinate of hit in microns
yrec2- (output) best estimate of second y-coordinate of hit in microns
sigmay- (output) best estimate of uncertainty on yrec in microns
prob2y- (output) probability describing goodness-of-fit for y-reco
xrec1- (output) best estimate of first x-coordinate of hit in microns
xrec2- (output) best estimate of second x-coordinate of hit in microns
sigmax- (output) best estimate of uncertainty on xrec in microns
prob2x- (output) probability describing goodness-of-fit for x-reco
q2bin- (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]

Definition at line 1477 of file SiPixelTemplateSplit.cc.

References PixelTempSplit(), and HLTSiStripMonitoring_cff::speed.

1493  {
1494  // Local variables
1495  const bool deadpix = false;
1496  const bool resolve = true;
1497  float dchisq, prob2Q;
1498  std::vector<std::pair<int, int> > zeropix;
1499  const int speed = 1;
1500 
1501  return SiPixelTemplateSplit::PixelTempSplit(id,
1502  cotalpha,
1503  cotbeta,
1504  cluster,
1505  ydouble,
1506  xdouble,
1507  templ,
1508  yrec1,
1509  yrec2,
1510  sigmay,
1511  prob2y,
1512  xrec1,
1513  xrec2,
1514  sigmax,
1515  prob2x,
1516  q2bin,
1517  prob2Q,
1518  resolve,
1519  speed,
1520  dchisq,
1521  deadpix,
1522  zeropix,
1523  templ2D);
1524 
1525 } // PixelTempSplit
SiPixelTemplateSplit::PixelTempSplit
int PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec1, float &yrec2, float &sigmay, float &prob2y, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, float &prob2Q, bool resolve, int speed, float &dchisq, bool deadpix, std::vector< std::pair< int, int > > &zeropix, SiPixelTemplate2D &templ2D)
Definition: SiPixelTemplateSplit.cc:98
HLTSiStripMonitoring_cff.speed
speed
Definition: HLTSiStripMonitoring_cff.py:22
Generated for CMSSW Reference Manual by doxygen 1.8.14