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, 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, 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, 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 &  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.

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

References PixelTempSplit(), and HLTSiStripMonitoring_cff::speed.

◆ PixelTempSplit() [2/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 97 of file SiPixelTemplateSplit.cc.

121 {
122  // Local variables
123  int i, j, k, binq, midpix, fypix, nypix, lypix, logypx, lparm;
124  int fxpix, nxpix, lxpix, logxpx, shifty, shiftx;
125  int nclusx, nclusy;
126  int nybin, ycbin, nxbin, xcbin, minbinj, minbink;
127  int deltaj, jmin, jmax, kmin, kmax, km, fxbin, lxbin, fybin, lybin, djy, djx;
128  float sythr, sxthr, delta, sigma, sigavg, pseudopix, xsize, ysize, qscale, lanpar[2][5];
129  float ss2, ssa, sa2, rat, fq, qtotal, qpixel, qavg;
130  float x1p, x2p, y1p, y2p, deltachi2;
131  float originx, originy, bias, maxpix, minmax;
132  double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
133  double hchi2, hndof, sigmal1, sigmal2, mpv1, mpv2, arg1, arg2, q05, like, loglike1, loglike2;
134  double prvav, mpv, sigmaQ, kappa, xvav, beta2;
135  float ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
136  float ysig2[BYSIZE], xsig2[BXSIZE];
137  float yw2[BYSIZE], xw2[BXSIZE], ysw[BYSIZE], xsw[BXSIZE];
138  float cluster2[BXM2][BYM2], temp2d1[BXM2][BYM2], temp2d2[BXM2][BYM2];
139  bool yd[BYSIZE], xd[BXSIZE], anyyd, anyxd, any2dfail;
140  const float sqrt2x = {3.0000}, sqrt2y = {1.7000};
141  const float sqrt12 = {3.4641};
142  const float probmin = {1.110223e-16};
143  const float prob2Qmin = {1.e-5};
144  std::pair<int, int> pixel;
145 
146  // bool SiPixelTemplateSplit::SimpleTemplate2D(float cotalpha, float cotbeta, float xhit, float yhit, float thick, float lorxwidth, float lorywidth,
147  // float qavg, std::vector<bool> ydouble, std::vector<bool> xdouble, float template2d[BXM2][BYM2]);
148 
149  // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
150 
151  const double mean1pix = {0.100}, chi21min = {0.160};
152 
153  // First, interpolate the template needed to analyze this cluster
154  // check to see of the track direction is in the physical range of the loaded template
155 
156  if (!templ.interpolate(id, cotalpha, cotbeta)) {
157  LOGDEBUG("SiPixelTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha
158  << ", cot(beta) = " << cotbeta
159  << " is not within the acceptance of template ID = " << id
160  << ", no reconstruction performed" << ENDL;
161  return 20;
162  }
163 
164  // Get pixel dimensions from the template (to allow multiple detectors in the future)
165 
166  xsize = templ.xsize();
167  ysize = templ.ysize();
168 
169  // Define size of pseudopixel
170 
171  pseudopix = templ.s50();
172  // q05 = 0.28*pseudopix;
173  q05 = 0.05f * pseudopix;
174 
175  // Get charge scaling factor
176 
177  qscale = templ.qscale();
178 
179  // Make a local copy of the cluster container so that we can muck with it
180 
181  array_2d cluster = clust;
182 
183  // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
184 
185  if (cluster.num_dimensions() != 2) {
186  LOGERROR("SiPixelTemplateReco") << "input cluster container (BOOST Multiarray) has wrong number of dimensions"
187  << ENDL;
188  return 3;
189  }
190  nclusx = (int)cluster.shape()[0];
191  nclusy = (int)cluster.shape()[1];
192  if (nclusx != (int)xdouble.size()) {
193  LOGERROR("SiPixelTemplateReco") << "input cluster container x-size is not equal to double pixel flag container size"
194  << ENDL;
195  return 4;
196  }
197  if (nclusy != (int)ydouble.size()) {
198  LOGERROR("SiPixelTemplateReco") << "input cluster container y-size is not equal to double pixel flag container size"
199  << ENDL;
200  return 5;
201  }
202 
203  // enforce maximum size
204 
205  if (nclusx > TXSIZE) {
206  nclusx = TXSIZE;
207  }
208  if (nclusy > TYSIZE) {
209  nclusy = TYSIZE;
210  }
211 
212  // First, rescale all pixel charges
213 
214  for (i = 0; i < nclusy; ++i) {
215  for (j = 0; j < nclusx; ++j) {
216  if (cluster[j][i] > 0) {
217  cluster[j][i] *= qscale;
218  }
219  }
220  }
221 
222  // Next, sum the total charge and "decapitate" big pixels
223 
224  qtotal = 0.f;
225  minmax = 2.0f * templ.pixmax();
226  for (i = 0; i < nclusy; ++i) {
227  maxpix = minmax;
228  if (ydouble[i]) {
229  maxpix *= 2.f;
230  }
231  for (j = 0; j < nclusx; ++j) {
232  qtotal += cluster[j][i];
233  if (cluster[j][i] > maxpix) {
234  cluster[j][i] = maxpix;
235  }
236  }
237  }
238 
239  // Do the cluster repair here
240 
241  if (deadpix) {
242  fypix = BYM3;
243  lypix = -1;
244  for (i = 0; i < nclusy; ++i) {
245  ysum[i] = 0.f;
246  // Do preliminary cluster projection in y
247  for (j = 0; j < nclusx; ++j) {
248  ysum[i] += cluster[j][i];
249  }
250  if (ysum[i] > 0.f) {
251  // identify ends of cluster to determine what the missing charge should be
252  if (i < fypix) {
253  fypix = i;
254  }
255  if (i > lypix) {
256  lypix = i;
257  }
258  }
259  }
260 
261  // Now loop over dead pixel list and "fix" everything
262 
263  //First see if the cluster ends are redefined and that we have only one dead pixel per column
264  int nyzero[TYSIZE]{};
265  std::vector<std::pair<int, int> >::const_iterator zeroIter = zeropix.begin(), zeroEnd = zeropix.end();
266  for (; zeroIter != zeroEnd; ++zeroIter) {
267  i = zeroIter->second;
268  if (i < 0 || i > TYSIZE - 1) {
269  LOGERROR("SiPixelTemplateReco") << "dead pixel column y-index " << i << ", no reconstruction performed" << ENDL;
270  return 11;
271  }
272 
273  // count the number of dead pixels in each column
274  ++nyzero[i];
275  // allow them to redefine the cluster ends
276  if (i < fypix) {
277  fypix = i;
278  }
279  if (i > lypix) {
280  lypix = i;
281  }
282  }
283 
284  nypix = lypix - fypix + 1;
285 
286  // 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
287 
288  for (zeroIter = zeropix.begin(); zeroIter != zeroEnd; ++zeroIter) {
289  i = zeroIter->second;
290  j = zeroIter->first;
291  if (j < 0 || j > TXSIZE - 1) {
292  LOGERROR("SiPixelTemplateReco") << "dead pixel column x-index " << j << ", no reconstruction performed" << ENDL;
293  return 12;
294  }
295  if ((i == fypix || i == lypix) && nypix > 1) {
296  maxpix = templ.symax() / 2.;
297  } else {
298  maxpix = templ.symax();
299  }
300  if (ydouble[i]) {
301  maxpix *= 2.;
302  }
303  if (nyzero[i] > 0 && nyzero[i] < 3) {
304  qpixel = (maxpix - ysum[i]) / (float)nyzero[i];
305  } else {
306  qpixel = 1.;
307  }
308  if (qpixel < 1.f) {
309  qpixel = 1.f;
310  }
311  cluster[j][i] = qpixel;
312  // Adjust the total cluster charge to reflect the charge of the "repaired" cluster
313  qtotal += qpixel;
314  }
315  // End of cluster repair section
316  }
317 
318  // Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container
319 
320  for (i = 0; i < BYSIZE; ++i) {
321  ysum[i] = 0.f;
322  yd[i] = false;
323  }
324  k = 0;
325  anyyd = false;
326  for (i = 0; i < nclusy; ++i) {
327  for (j = 0; j < nclusx; ++j) {
328  ysum[k] += cluster[j][i];
329  }
330 
331  // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
332 
333  if (ydouble[i]) {
334  ysum[k] /= 2.f;
335  ysum[k + 1] = ysum[k];
336  yd[k] = true;
337  yd[k + 1] = false;
338  k = k + 2;
339  anyyd = true;
340  } else {
341  yd[k] = false;
342  ++k;
343  }
344  if (k > BYM1) {
345  break;
346  }
347  }
348 
349  // Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container
350 
351  for (i = 0; i < BXSIZE; ++i) {
352  xsum[i] = 0.f;
353  xd[i] = false;
354  }
355  k = 0;
356  anyxd = false;
357  for (j = 0; j < nclusx; ++j) {
358  for (i = 0; i < nclusy; ++i) {
359  xsum[k] += cluster[j][i];
360  }
361 
362  // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
363 
364  if (xdouble[j]) {
365  xsum[k] /= 2.f;
366  xsum[k + 1] = xsum[k];
367  xd[k] = true;
368  xd[k + 1] = false;
369  k = k + 2;
370  anyxd = true;
371  } else {
372  xd[k] = false;
373  ++k;
374  }
375  if (k > BXM1) {
376  break;
377  }
378  }
379 
380  // next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx
381 
382  fypix = -1;
383  nypix = 0;
384  lypix = 0;
385  logypx = 0;
386  for (i = 0; i < BYSIZE; ++i) {
387  if (ysum[i] > 0.) {
388  if (fypix == -1) {
389  fypix = i;
390  }
391  if (!yd[i]) {
392  ysort[logypx] = ysum[i];
393  ++logypx;
394  }
395  ++nypix;
396  lypix = i;
397  }
398  }
399 
400  // Make sure cluster is continuous
401 
402  if ((lypix - fypix + 1) != nypix || nypix == 0) {
403  LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster doesn't agree with number of pixels above threshold"
404  << ENDL;
405  if (theVerboseLevel > 2) {
406  LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
407  for (i = 0; i < BYSIZE - 1; ++i) {
408  LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";
409  }
410  LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE - 1] << ENDL;
411  }
412 
413  return 1;
414  }
415 
416  // If cluster is longer than max template size, technique fails
417 
418  if (nypix > TYSIZE) {
419  LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster is larger than maximum template size" << ENDL;
420  if (theVerboseLevel > 2) {
421  LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
422  for (i = 0; i < BYSIZE - 1; ++i) {
423  LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";
424  }
425  LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE - 1] << ENDL;
426  }
427 
428  return 6;
429  }
430 
431  // next, center the cluster on template center if necessary
432 
433  midpix = (fypix + lypix) / 2;
434  // shifty = BHY - midpix;
435  shifty = templ.cytemp() - midpix;
436  if (shifty > 0) {
437  for (i = lypix; i >= fypix; --i) {
438  ysum[i + shifty] = ysum[i];
439  ysum[i] = 0.;
440  yd[i + shifty] = yd[i];
441  yd[i] = false;
442  }
443  } else if (shifty < 0) {
444  for (i = fypix; i <= lypix; ++i) {
445  ysum[i + shifty] = ysum[i];
446  ysum[i] = 0.;
447  yd[i + shifty] = yd[i];
448  yd[i] = false;
449  }
450  }
451  lypix += shifty;
452  fypix += shifty;
453 
454  // If the cluster boundaries are OK, add pesudopixels, otherwise quit
455 
456  if (fypix > 1 && fypix < BYM2) {
457  ysum[fypix - 1] = pseudopix;
458  ysum[fypix - 2] = 0.2f * pseudopix;
459  } else {
460  return 8;
461  }
462  if (lypix > 1 && lypix < BYM2) {
463  ysum[lypix + 1] = pseudopix;
464  ysum[lypix + 2] = 0.2f * pseudopix;
465  } else {
466  return 8;
467  }
468 
469  // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
470 
471  if (ydouble[0]) {
472  originy = -0.5f;
473  } else {
474  originy = 0.f;
475  }
476 
477  // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx
478 
479  fxpix = -1;
480  nxpix = 0;
481  lxpix = 0;
482  logxpx = 0;
483  for (i = 0; i < BXSIZE; ++i) {
484  if (xsum[i] > 0.) {
485  if (fxpix == -1) {
486  fxpix = i;
487  }
488  if (!xd[i]) {
489  xsort[logxpx] = xsum[i];
490  ++logxpx;
491  }
492  ++nxpix;
493  lxpix = i;
494  }
495  }
496 
497  // Make sure cluster is continuous
498 
499  if ((lxpix - fxpix + 1) != nxpix) {
500  LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold"
501  << ENDL;
502  if (theVerboseLevel > 2) {
503  LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
504  for (i = 0; i < BXSIZE - 1; ++i) {
505  LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";
506  }
507  LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE - 1] << ENDL;
508  }
509 
510  return 2;
511  }
512 
513  // If cluster is longer than max template size, technique fails
514 
515  if (nxpix > TXSIZE) {
516  LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
517  if (theVerboseLevel > 2) {
518  LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
519  for (i = 0; i < BXSIZE - 1; ++i) {
520  LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";
521  }
522  LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE - 1] << ENDL;
523  }
524 
525  return 7;
526  }
527 
528  // next, center the cluster on template center if necessary
529 
530  midpix = (fxpix + lxpix) / 2;
531  // shiftx = BHX - midpix;
532  shiftx = templ.cxtemp() - midpix;
533  if (shiftx > 0) {
534  for (i = lxpix; i >= fxpix; --i) {
535  xsum[i + shiftx] = xsum[i];
536  xsum[i] = 0.f;
537  xd[i + shiftx] = xd[i];
538  xd[i] = false;
539  }
540  } else if (shiftx < 0) {
541  for (i = fxpix; i <= lxpix; ++i) {
542  xsum[i + shiftx] = xsum[i];
543  xsum[i] = 0.f;
544  xd[i + shiftx] = xd[i];
545  xd[i] = false;
546  }
547  }
548  lxpix += shiftx;
549  fxpix += shiftx;
550 
551  // If the cluster boundaries are OK, add pesudopixels, otherwise quit
552 
553  if (fxpix > 1 && fxpix < BXM2) {
554  xsum[fxpix - 1] = pseudopix;
555  xsum[fxpix - 2] = 0.2f * pseudopix;
556  } else {
557  return 9;
558  }
559  if (lxpix > 1 && lxpix < BXM2) {
560  xsum[lxpix + 1] = pseudopix;
561  xsum[lxpix + 2] = 0.2f * pseudopix;
562  } else {
563  return 9;
564  }
565 
566  // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
567 
568  if (xdouble[0]) {
569  originx = -0.5f;
570  } else {
571  originx = 0.f;
572  }
573 
574  // uncertainty and final corrections depend upon total charge bin
575 
576  qavg = templ.qavg();
577  fq = qtotal / qavg;
578  if (fq > 3.0f) {
579  binq = 0;
580  } else {
581  if (fq > 2.0f) {
582  binq = 1;
583  } else {
584  if (fq > 1.70f) {
585  binq = 2;
586  } else {
587  binq = 3;
588  }
589  }
590  }
591 
592  // Calculate the Vavilov probability that the cluster charge is consistent with a merged cluster
593 
594  if (speed < 0) {
595  templ.vavilov2_pars(mpv, sigmaQ, kappa);
596 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
597  if ((sigmaQ <= 0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
598  throw cms::Exception("DataCorrupt") << "SiPixelTemplateSplit::Vavilov parameters mpv/sigmaQ/kappa = " << mpv
599  << "/" << sigmaQ << "/" << kappa << std::endl;
600  }
601 #else
602  assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
603 #endif
604  xvav = ((double)qtotal - mpv) / sigmaQ;
605  beta2 = 1.;
606  // VVIObj is a private port of CERNLIB VVIDIS
607  VVIObj vvidist(kappa, beta2, 1);
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

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, kappa, SiPixelTemplate2D::landau_par(), dqm-mbProfile::log, LOGDEBUG, LOGERROR, maxpix, 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().

◆ 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,
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.

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

References PixelTempSplit(), and HLTSiStripMonitoring_cff::speed.

◆ 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.

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

References PixelTempSplit(), and HLTSiStripMonitoring_cff::speed.

◆ 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.

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

References PixelTempSplit(), and HLTSiStripMonitoring_cff::speed.

SiPixelTemplate::sxmax
float sxmax()
average pixel signal for x-projection of cluster
Definition: SiPixelTemplate.h:417
BXM3
#define BXM3
Definition: SiPixelTemplateDefs.h:37
SiPixelTemplate::dxtwo
float dxtwo()
mean offset/correction for one double-pixel x-clusters
Definition: SiPixelTemplate.h:420
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:982
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
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
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
BXSIZE
#define BXSIZE
Definition: SiPixelTemplateDefs.h:33
min
T min(T a, T b)
Definition: MathUtil.h:58
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
muonClassificationByHits_cfi.pixel
pixel
Definition: muonClassificationByHits_cfi.py:9
cms::cuda::assert
assert(be >=bs)
ysize
const Int_t ysize
Definition: trackSplitPlot.h:43
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
SiPixelTemplate::cytemp
int cytemp()
Return central pixel of y template pixels above readout threshold.
Definition: SiPixelTemplate.cc:2649
SiPixelTemplate::ysize
float ysize()
pixel y-size (microns)
Definition: SiPixelTemplate.h:676
xsize
const Int_t xsize
Definition: trackSplitPlot.h:42
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:97
SiPixelTemplate::chi2xavgc2m
float chi2xavgc2m(int i)
1st pass chi2 min search: average x-chisq for merged clusters
Definition: SiPixelTemplate.h:624
TXSIZE
#define TXSIZE
Definition: SiPixelTemplateDefs.h:30
BXM1
#define BXM1
Definition: SiPixelTemplateDefs.h:35
SiPixelTemplate::s50
float s50()
1/2 of the pixel threshold signal in electrons
Definition: SiPixelTemplate.h:410
SiPixelTemplate::vavilov2_pars
void vavilov2_pars(double &mpv, double &sigma, double &kappa)
Definition: SiPixelTemplate.cc:4079
BYM2
#define BYM2
Definition: SiPixelTemplateDefs.h:28
BYM1
#define BYM1
Definition: SiPixelTemplateDefs.h:27
SiPixelTemplate::sytwo
float sytwo()
rms for one double-pixel y-clusters
Definition: SiPixelTemplate.h:416
BXM2
#define BXM2
Definition: SiPixelTemplateDefs.h:36
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
theVerboseLevel
static constexpr int theVerboseLevel
Definition: SiPixelTemplateSplit.cc:50
BHX
#define BHX
Definition: SiPixelTemplateDefs.h:34
SiPixelTemplate::dytwo
float dytwo()
mean offset/correction for one double-pixel y-clusters
Definition: SiPixelTemplate.h:415
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
dqmdumpme.k
k
Definition: dqmdumpme.py:60
SiPixelTemplate::chi2yavgone
float chi2yavgone()
//!< average y chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:657
BYM3
#define BYM3
Definition: SiPixelTemplateDefs.h:29
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
ENDL
#define ENDL
Definition: SiPixelTemplateSplit.cc:51
SiPixelTemplate::xtemp3d
void xtemp3d(int j, int k, std::vector< float > &xtemplate)
Definition: SiPixelTemplate.cc:2937
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
LOGDEBUG
#define LOGDEBUG(x)
Definition: SiPixelTemplateSplit.cc:49
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
SiPixelTemplate::ytemp3d_int
void ytemp3d_int(int nypix, int &nybins)
Definition: SiPixelTemplate.cc:2757
dumpMFGeometry_cfg.delta
delta
Definition: dumpMFGeometry_cfg.py:25
createfilelist.int
int
Definition: createfilelist.py:10
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
VVIObj
Definition: VVIObj.h:24
SiPixelTemplate::xsize
float xsize()
pixel x-size (microns)
Definition: SiPixelTemplate.h:675
HLTSiStripMonitoring_cff.speed
speed
Definition: HLTSiStripMonitoring_cff.py:22
LOGERROR
#define LOGERROR(x)
Definition: SiPixelTemplateSplit.cc:48
BHY
#define BHY
Definition: SiPixelTemplateDefs.h:26
SiPixelTemplate::qscale
float qscale()
charge scaling factor
Definition: SiPixelTemplate.h:409
SiPixelTemplateSplit::array_2d
boost::multi_array< float, 2 > array_2d
Definition: SiPixelTemplateSplit.h:47
SiPixelTemplate::chi2yavgc2m
float chi2yavgc2m(int i)
1st pass chi2 min search: average y-chisq for merged clusters
Definition: SiPixelTemplate.h:582
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:1751
Exception
Definition: hltDiff.cc:245
SiPixelTemplate::ytemp3d
void ytemp3d(int j, int k, std::vector< float > &ytemplate)
Definition: SiPixelTemplate.cc:2835
SiPixelTemplate::sxtwo
float sxtwo()
rms for one double-pixel x-clusters
Definition: SiPixelTemplate.h:421
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])
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
SiPixelTemplate::xtemp3d_int
void xtemp3d_int(int nxpix, int &nxbins)
Definition: SiPixelTemplate.cc:2858
kappa
static const G4double kappa
Definition: UrbanMscModel93.cc:35
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::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:3763
SiPixelTemplate::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
Definition: SiPixelTemplate.cc:1327
SiPixelTemplate::qmin
float qmin()
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
Definition: SiPixelTemplate.h:422
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16