CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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

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

Definition at line 50 of file SiPixelTemplateSplit.h.

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

Definition at line 51 of file SiPixelTemplateSplit.h.

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

Definition at line 52 of file SiPixelTemplateSplit.h.

Function Documentation

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.

References 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(), delta, SiPixelTemplate::dxone(), SiPixelTemplate::dxtwo(), SiPixelTemplate::dyone(), SiPixelTemplate::dytwo(), alignCSCRings::e, ENDL, edm::hlt::Exception, f, VVIObj::fcn(), Gamma, i, SiPixelTemplate::interpolate(), j, roll_playback::k, kappa, SiPixelTemplate2D::landau_par(), create_public_lumi_plots::log, LOGDEBUG, LOGERROR, bookConverter::max, maxpix, min(), SiPixelTemplate::pixmax(), SiPixelTemplate::qavg(), SiPixelTemplate::qmin(), SiPixelTemplate::qscale(), SiPixelTemplate::s50(), SiPixelTemplate::simpletemplate2D(), python.multivaluedict::sort(), 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 PixelCPETemplateReco::localPosition(), PixelTempSplit(), and TrackClusterSplitter::splitCluster().

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

References PixelTempSplit().

1096 {
1097  // Local variables
1098  const bool deadpix = false;
1099  std::vector<std::pair<int, int> > zeropix;
1100 
1101  return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ,
1102  yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);
1103 
1104 } // 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:97
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 1132 of file SiPixelTemplateSplit.cc.

References PixelTempSplit().

1137 {
1138  // Local variables
1139  const bool deadpix = false;
1140  std::vector<std::pair<int, int> > zeropix;
1141  const int speed = 1;
1142 
1143  return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ,
1144  yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);
1145 
1146 } // 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:97
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 1171 of file SiPixelTemplateSplit.cc.

References PixelTempSplit().

1176 {
1177  // Local variables
1178  const bool deadpix = false;
1179  const bool resolve = true;
1180  float dchisq;
1181  std::vector<std::pair<int, int> > zeropix;
1182  const int speed = 1;
1183 
1184  return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ,
1185  yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);
1186 
1187 } // 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:97
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 1213 of file SiPixelTemplateSplit.cc.

References PixelTempSplit().

1218 {
1219  // Local variables
1220  const bool deadpix = false;
1221  const bool resolve = true;
1222  float dchisq, prob2Q;
1223  std::vector<std::pair<int, int> > zeropix;
1224  const int speed = 1;
1225 
1226  return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ,
1227  yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);
1228 
1229 } // 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:97
Generated for CMSSW Reference Manual by doxygen 1.8.5