CMS 3D CMS Logo
Classes | Functions
SiPixelTemplateReco2D Namespace Reference

Classes

struct  ClusMatrix
 

Functions

int PixelTempReco2D (int id, float cotalpha, float cotbeta, float locBz, float locBx, int edgeflagy, int edgeflagx, ClusMatrix &cluster, SiPixelTemplate2D &templ, float &yrec, float &sigmay, float &xrec, float &sigmax, float &probxy, float &probQ, int &qbin, float &deltay, int &npixel)
 
int PixelTempReco2D (int id, float cotalpha, float cotbeta, float locBz, float locBx, int edgeflagy, int edgeflagx, ClusMatrix &cluster, SiPixelTemplate2D &templ, float &yrec, float &sigmay, float &xrec, float &sigmax, float &probxy, float &probQ, int &qbin, float &deltay)
 

Function Documentation

int SiPixelTemplateReco2D::PixelTempReco2D ( int  id,
float  cotalpha,
float  cotbeta,
float  locBz,
float  locBx,
int  edgeflagy,
int  edgeflagx,
ClusMatrix &  cluster,
SiPixelTemplate2D &  templ2D,
float &  yrec,
float &  sigmay,
float &  xrec,
float &  sigmax,
float &  probxy,
float &  probQ,
int &  qbin,
float &  deltay,
int &  npixels 
)

Reconstruct the best estimate of the hit position for pixel clusters.

Parameters
id- (input) identifier of the template to use
cotalpha- (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
cotbeta- (input) the cotangent of the beta track angle (see CMS IN 2004/014)
locBz- (input) the sign of this quantity is used to determine whether to flip cot(beta)<0 quantities from cot(beta)>0 (FPix only) for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0 for Phase 1 FPix IP-related tracks, see next comment
locBx- (input) the sign of this quantity is used to determine whether to flip cot(alpha/beta)<0 quantities from cot(alpha/beta)>0 (FPix only) for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0 for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0
edgeflagy- (input) flag to indicate the present of edges in y: 0-none (or interior gap), 1-edge at small y, 2-edge at large y, 3-edge at either end
edgeflagx- (input) flag to indicate the present of edges in x: 0-none, 1-edge at small x, 2-edge at large x
cluster- (input) pixel array struct with double pixel flags and bounds origin of local coords (0,0) at center of pixel cluster[0][0].
templ2D- (input) the 2D template used in the reconstruction
yrec- (output) best estimate of y-coordinate of hit in microns
sigmay- (output) best estimate of uncertainty on yrec in microns
xrec- (output) best estimate of x-coordinate of hit in microns
sigmax- (output) best estimate of uncertainty on xrec in microns
probxy- (output) probability describing goodness-of-fit
probQ- (output) probability describing upper cluster charge tail
qbin- (output) index (0-4) describing the charge of the cluster qbin = 0 Q/Q_avg > 1.5 [few % of all hits] 1 1.5 > Q/Q_avg > 1.0 [~30% of all hits] 2 1.0 > Q/Q_avg > 0.85 [~30% of all hits] 3 0.85 > Q/Q_avg > min1 [~30% of all hits]
deltay- (output) template y-length - cluster length [when > 0, possibly missing end]

Definition at line 83 of file SiPixelTemplateReco2D.cc.

References funct::abs(), BXM2, BYM2, vertices_cff::chi2, SiPixelTemplate2D::chi2ppix(), SiPixelTemplate2D::chi2scale(), SiPixelTemplate2D::clsleny(), gather_cfg::cout, MillePedeFileConverter_cfg::e, Exception, f, SiPixelTemplate2D::fbin(), VVIObjF::fcn(), objects.autophobj::float, Gamma, mps_fire::i, createfilelist::int, SiPixelTemplate2D::interpolate(), jlist, gen::k, kappa, SiPixelTemplate2D::kappavav(), SiPixelTemplate2D::lorxdrift(), SiPixelTemplate2D::lorydrift(), SiPixelTemplateReco2D::ClusMatrix::matrix, maxpix, SiPixelTemplateReco2D::ClusMatrix::mcol, SiPixelTemplate2D::mpvvav(), SiPixelTemplateReco2D::ClusMatrix::mrow, SiPixelTemplate2D::offsetx(), SiPixelTemplate2D::offsety(), digitizers_cfi::pixel, SiPixelTemplate2D::pixmax(), EnergyCorrector::pt, SiPixelTemplate2D::qavg(), SiPixelTemplate2D::qscale(), SiPixelTemplate2D::s50(), SiPixelTemplate2D::scalex(), SiPixelTemplate2D::scaley(), SiPixelTemplate2D::sigmavav(), mathSSE::sqrt(), SiPixelTemplate2D::sxymax(), T2HYP1, THXP1, THYP1, SiPixelTemplateReco2D::ClusMatrix::xdouble, xsize, SiPixelTemplate2D::xsize(), SiPixelTemplate2D::xysigma2(), SiPixelTemplate2D::xytemp(), SiPixelTemplateReco2D::ClusMatrix::ydouble, ysize, and SiPixelTemplate2D::ysize().

Referenced by PixelCPEClusterRepair::callTempReco2D(), and PixelTempReco2D().

87 {
88  // Local variables
89 
90  float template2d[BXM2][BYM2], dpdx2d[2][BXM2][BYM2], fbin[3];
91 
92  // fraction of truncation signal to measure the cluster ends
93  const float fracpix = 0.45f;
94 
95  const int nilist = 9, njlist = 5;
96  const float ilist[nilist] = {0.f, -1.f, -0.75f, -0.5f, -0.25f, 0.25f, 0.5f, 0.75f, 1.f};
97  const float jlist[njlist] = {0.f, -0.5f, -0.25f, 0.25f, 0.50f};
98 
99 
100  // Extract some relevant info from the 2D template
101 
102  if(id > 0) {if(!templ2D.interpolate(id, cotalpha, cotbeta, locBz, locBx)) return 4;}
103  float xsize = templ2D.xsize();
104  float ysize = templ2D.ysize();
105 
106  // Allow Qbin Q/Q_avg fractions to vary to optimize error estimation
107 
108  for(int i=0; i<3; ++i) {fbin[i] = templ2D.fbin(i);}
109 
110  float q50 = templ2D.s50();
111  float pseudopix = 0.2f*q50;
112  float pseudopix2 = q50*q50;
113 
114 
115  // Get charge scaling factor
116 
117  float qscale = templ2D.qscale();
118 
119  // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
120 
121  int nclusx = cluster.mrow;
122  int nclusy = (int)cluster.mcol;
123  bool * xdouble = cluster.xdouble;
124  bool * ydouble = cluster.ydouble;
125 
126  // First, rescale all pixel charges and compute total charge
127  float qtotal = 0.f;
128  int imin=BYM2, imax=0, jmin=BXM2, jmax=0;
129  for(int k=0; k<nclusx*nclusy; ++k) {
130  cluster.matrix[k] *= qscale;
131  qtotal +=cluster.matrix[k];
132  int j = k/nclusy;
133  int i = k - j*nclusy;
134  if(cluster.matrix[k] > 0.f) {
135  if(j < jmin) {jmin = j;}
136  if(j > jmax) {jmax = j;}
137  if(i < imin) {imin = i;}
138  if(i > imax) {imax = i;}
139  }
140  }
141 
142 // Calculate the shifts needed to center the cluster in the buffer
143 
144  int shiftx = THXP1 - (jmin+jmax)/2;
145  int shifty = THYP1 - (imin+imax)/2;
146 // Always shift by at least one pixel
147  if(shiftx < 1) shiftx = 1;
148  if(shifty < 1) shifty = 1;
149 // Shift the min maxes too
150  jmin += shiftx;
151  jmax += shiftx;
152  imin += shifty;
153  imax += shifty;
154 
155  // uncertainty and final corrections depend upon total charge bin
156 
157  float fq0 = qtotal/templ2D.qavg();
158 
159  // Next, copy the cluster and double pixel flags into a shifted workspace to
160  // allow for significant zeros [pseudopixels] to be added
161 
162  float clusxy[BXM2][BYM2];
163  for(int j=0; j<BXM2; ++j) {for(int i=0; i<BYM2; ++i) {clusxy[j][i] = 0.f;}}
164 
165  const unsigned int NPIXMAX = 200;
166 
167  int indexxy[2][NPIXMAX];
168  float pixel[NPIXMAX];
169  float sigma2[NPIXMAX];
170  float minmax = templ2D.pixmax();
171  float ylow0 = 0.f, yhigh0 = 0.f, xlow0 = 0.f, xhigh0 = 0.f;
172  int npixel = 0;
173  float ysum[BYM2], ye[BYM2+1], ypos[BYM2], xpos[BXM2], xe[BXM2+1];
174  bool yd[BYM2], xd[BXM2];
175 
176 // Copy double pixel flags to shifted arrays
177 
178  for(int i=0; i<BYM2; ++i) {ysum[i] = 0.f; yd[i] = false;}
179  for(int j=0; j<BXM2; ++j) {xd[j] = false;}
180  for(int i=0; i<nclusy; ++i) {
181  if(ydouble[i]) {
182  int iy = i+shifty;
183  if(iy > -1 && iy < BYM2) yd[iy] = true;
184  }
185  }
186  for(int j=0; j<nclusx; ++j) {
187  if(xdouble[j]) {
188  int jx = j+shiftx;
189  if(jx > -1 && jx < BXM2) xd[jx] = true;
190  }
191  }
192 // Work out the positions of the rows+columns relative to the lower left edge of pixel[1,1]
193  xe[0] = -xsize;
194  ye[0] = -ysize;
195  for(int i=0; i<BYM2; ++i) {
196  float ypitch = ysize;
197  if(yd[i]) {ypitch += ysize;}
198  ye[i+1] = ye[i] + ypitch;
199  ypos[i] = ye[i] + ypitch/2.;
200  }
201  for(int j=0; j<BXM2; ++j) {
202  float xpitch = xsize;
203  if(xd[j]) {xpitch += xsize;}
204  xe[j+1] = xe[j] + xpitch;
205  xpos[j] = xe[j] + xpitch/2.;
206  }
207 // Shift the cluster center to the central pixel of the array, truncate big signals
208  for(int i=0; i<nclusy; ++i) {
209  int iy = i+shifty;
210  float maxpix = minmax;
211  if(ydouble[i]) {maxpix *=2.f;}
212  if(iy > -1 && iy < BYM2) {
213  for(int j=0; j<nclusx; ++j) {
214  int jx = j+shiftx;
215  if(jx > -1 && jx < BXM2) {
216  if(cluster(j,i) > maxpix) {
217  clusxy[jx][iy] = maxpix;
218  } else {
219  clusxy[jx][iy] = cluster(j,i);
220  }
221  if(clusxy[jx][iy] > 0.f) {
222  ysum[iy] += clusxy[jx][iy];
223  indexxy[0][npixel] = jx;
224  indexxy[1][npixel] = iy;
225  pixel[npixel] = clusxy[jx][iy];
226  ++npixel;
227  }
228  }
229  }
230  }
231  }
232 
233 // Make sure that we find at least one pixel
234  if(npixel < 1) {
235 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
236  throw cms::Exception("DataCorrupt") << "PixelTemplateReco2D::number of pixels above threshold = " << npixel << std::endl;
237 #else
238  std::cout << "PixelTemplateReco2D::number of pixels above threshold = " << npixel << std::endl;
239 #endif
240  return 1;
241  }
242 
243 // Get the shifted coordinates of the cluster ends
244  xlow0 = xe[jmin];
245  xhigh0 = xe[jmax+1];
246  ylow0 = ye[imin];
247  yhigh0 = ye[imax+1];
248 
249  // Next, calculate the error^2 [need to know which is the middle y pixel of the cluster]
250 
251  int ypixoff = T2HYP1 - (imin+imax)/2;
252  for(int k=0; k<npixel; ++k) {
253  int ypixeff = ypixoff + indexxy[1][k];
254  templ2D.xysigma2(pixel[k], ypixeff, sigma2[k]);
255  }
256 
257  // Next, find the half-height edges of the y-projection and identify any
258  // missing columns to remove from the fit
259 
260  int imisslow = -1, imisshigh = -1, jmisslow = -1, jmisshigh = -1;
261  float ylow = -1.f, yhigh = -1.f;
262  float hmaxpix = fracpix*templ2D.sxymax();
263  for(int i=imin; i<=imax; ++i) {
264  if(ysum[i] > hmaxpix && ysum[i-1] < hmaxpix && ylow < 0.f) {
265  ylow = ypos[i-1] + (ypos[i]-ypos[i-1])*(hmaxpix-ysum[i-1])/(ysum[i]-ysum[i-1]);
266  }
267  if(ysum[i] < q50) {
268  if(imisslow < 0) imisslow = i;
269  imisshigh = i;
270  }
271  if(ysum[i] > hmaxpix && ysum[i+1] < hmaxpix) {
272  yhigh = ypos[i] + (ypos[i+1]-ypos[i])*(ysum[i]-hmaxpix)/(ysum[i]-ysum[i+1]);
273  }
274  }
275  if(ylow < 0.f || yhigh < 0.f) {
276  ylow = ylow0;
277  yhigh = yhigh0;
278  }
279 
280  float templeny = templ2D.clsleny();
281  deltay = templeny - (yhigh - ylow)/ysize;
282 
283 
284  // x0 and y0 are best guess seeds for the fit
285 
286  float x0 = 0.5f*(xlow0 + xhigh0) - templ2D.lorxdrift();
287  float y0 = 0.5f*(ylow + yhigh) - templ2D.lorydrift();
288 // float y1 = yhigh - halfy - templ2D.lorydrift();
289 // printf("y0 = %f, y1 = %f \n", y0, y1);
290  // float y0 = 0.5f*(ylow + yhigh);
291 
292  // If there are missing edge columns, set up missing column flags and number
293  // of minimization passes
294 
295  int npass = 1;
296 
297  switch(edgeflagy) {
298  case 0:
299  break;
300  case 1:
301  imisshigh = imin-1;
302  imisslow = -1;
303  break;
304  case 2:
305  imisshigh = -1;
306  imisslow = imax+1;
307  break;
308  case 3:
309  imisshigh = imin-1;
310  imisslow = -1;
311  npass = 2;
312  break;
313  default:
314 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
315  throw cms::Exception("DataCorrupt") << "PixelTemplateReco2D::illegal edgeflagy = " << edgeflagy << std::endl;
316 #else
317  std::cout << "PixelTemplate:2D:illegal edgeflagy = " << edgeflagy << std::endl;
318 #endif
319  }
320 
321  switch(edgeflagx) {
322  case 0:
323  break;
324  case 1:
325  jmisshigh = jmin-1;
326  jmisslow = -1;
327  break;
328  case 2:
329  jmisshigh = -1;
330  jmisslow = jmax+1;
331  break;
332  default:
333 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
334  throw cms::Exception("DataCorrupt") << "PixelTemplateReco2D::illegal edgeflagx = " << edgeflagx << std::endl;
335 #else
336  std::cout << "PixelTemplate:2D:illegal edgeflagx = " << edgeflagx << std::endl;
337 #endif
338  }
339 
340  // Define quantities to be saved for each pass
341 
342 
343  float chi2min[2], xerr2[2], yerr2[2];
344  float x2D0[2], y2D0[2], qtfrac0[2];
345  int ipass, tpixel;
346  //int niter0[2];
347 
348  for(ipass = 0; ipass < npass; ++ipass) {
349 
350  if(ipass == 1) {
351 
352  // Now try again if both edges are possible
353 
354  imisshigh = -1;
355  imisslow = imax+1;
356 // erase pseudo pixels from previous pass
357  for(int k=npixel; k<tpixel; ++k) {
358  int j = indexxy[0][k];
359  int i = indexxy[1][k];
360  clusxy[j][i] = 0.f;
361  }
362  }
363 
364  // Next, add pseudo pixels around the periphery of the cluster
365 
366  tpixel = npixel;
367  for(int k=0; k<npixel; ++k) {
368  int j = indexxy[0][k];
369  int i = indexxy[1][k];
370  if((j-1) != jmisshigh) {
371  if(clusxy[j-1][i] < pseudopix) {
372  indexxy[0][tpixel] = j-1;
373  indexxy[1][tpixel] = i;
374  clusxy[j-1][i] = pseudopix;
375  pixel[tpixel] = pseudopix;
376  sigma2[tpixel] = pseudopix2;
377  ++tpixel;
378  }
379  }
380  if((j+1) != jmisslow) {
381  if(clusxy[j+1][i] < pseudopix) {
382  indexxy[0][tpixel] = j+1;
383  indexxy[1][tpixel] = i;
384  clusxy[j+1][i] = pseudopix;
385  pixel[tpixel] = pseudopix;
386  sigma2[tpixel] = pseudopix2;
387  ++tpixel;
388  }
389  }
390  // Don't add them if this is a dead column
391  if((i+1) != imisslow) {
392  if((j-1) != jmisshigh) {
393  if(clusxy[j-1][i+1] < pseudopix) {
394  indexxy[0][tpixel] = j-1;
395  indexxy[1][tpixel] = i+1;
396  clusxy[j-1][i+1] = pseudopix;
397  pixel[tpixel] = pseudopix;
398  sigma2[tpixel] = pseudopix2;
399  ++tpixel;
400  }
401  }
402  if(clusxy[j][i+1] < pseudopix) {
403  indexxy[0][tpixel] = j;
404  indexxy[1][tpixel] = i+1;
405  clusxy[j][i+1] = pseudopix;
406  pixel[tpixel] = pseudopix;
407  sigma2[tpixel] = pseudopix2;
408  ++tpixel;
409  }
410  if((j+1) != jmisslow) {
411  if(clusxy[j+1][i+1] < pseudopix) {
412  indexxy[0][tpixel] = j+1;
413  indexxy[1][tpixel] = i+1;
414  clusxy[j+1][i+1] = pseudopix;
415  pixel[tpixel] = pseudopix;
416  sigma2[tpixel] = pseudopix2;
417  ++tpixel;
418  }
419  }
420  }
421  // Don't add them if this is a dead column
422  if((i-1) != imisshigh) {
423  if((j-1) != jmisshigh) {
424  if(clusxy[j-1][i-1] < pseudopix) {
425  indexxy[0][tpixel] = j-1;
426  indexxy[1][tpixel] = i-1;
427  clusxy[j-1][i-1] = pseudopix;
428  pixel[tpixel] = pseudopix;
429  sigma2[tpixel] = pseudopix2;
430  ++tpixel;
431  }
432  }
433  if(clusxy[j][i-1] < pseudopix) {
434  indexxy[0][tpixel] = j;
435  indexxy[1][tpixel] = i-1;
436  clusxy[j][i-1] = pseudopix;
437  pixel[tpixel] = pseudopix;
438  sigma2[tpixel] = pseudopix2;
439  ++tpixel;
440  }
441  if((j+1) != jmisslow) {
442  if(clusxy[j+1][i-1] < pseudopix) {
443  indexxy[0][tpixel] = j+1;
444  indexxy[1][tpixel] = i-1;
445  clusxy[j+1][i-1] = pseudopix;
446  pixel[tpixel] = pseudopix;
447  sigma2[tpixel] = pseudopix2;
448  ++tpixel;
449  }
450  }
451  }
452  }
453 
454 
455  // Calculate chi2 over a grid of several seeds and choose the smallest
456 
457  chi2min[ipass] = 1000000.f;
458  float chi2, qtemplate, qactive, qtfrac = 0.f, x2D = 0.f, y2D = 0.f;
459 // Scale the y search grid for long clusters [longer than 7 pixels]
460  float ygridscale = 0.271*cotbeta;
461  if(ygridscale < 1.f) ygridscale = 1.f;
462  for(int is = 0; is<nilist; ++is) {
463  for(int js = 0; js<njlist; ++js) {
464  float xtry = x0 + jlist[js]*xsize;
465  float ytry = y0 + ilist[is]*ygridscale*ysize;
466  chi2 = 0.f;
467  qactive = 0.f;
468  for(int j=0; j<BXM2; ++j) {for(int i=0; i<BYM2; ++i) {template2d[j][i] = 0.f;}}
469  templ2D.xytemp(xtry, ytry, yd, xd, template2d, false, dpdx2d, qtemplate);
470  for(int k=0; k<tpixel; ++k) {
471  int jpix = indexxy[0][k];
472  int ipix = indexxy[1][k];
473  float qtpixel = template2d[jpix][ipix];
474  float pt = pixel[k]-qtpixel;
475  chi2 += pt*pt/sigma2[k];
476  if(k < npixel) {qactive += qtpixel;}
477  }
478  if(chi2 < chi2min[ipass]) {
479  chi2min[ipass] = chi2;
480  x2D = xtry;
481  y2D = ytry;
482  qtfrac = qactive/qtemplate;
483  }
484  }
485  }
486 
487  int niter = 0;
488  float xstep = 1.0f, ystep = 1.0f;
489  float minv11 = 1000.f, minv12 = 1000.f, minv22 = 1000.f;
490  chi2 = chi2min[ipass];
491  while(chi2 <= chi2min[ipass] && niter < 15 && (niter < 2 || (std::abs(xstep) > 0.2 || std::abs(ystep) > 0.2))) {
492 
493  // Remember the present parameters
494  x2D0[ipass] = x2D;
495  y2D0[ipass] = y2D;
496  qtfrac0[ipass] = qtfrac;
497  xerr2[ipass] = minv11;
498  yerr2[ipass] = minv22;
499  chi2min[ipass] = chi2;
500  //niter0[ipass] = niter;
501 
502  // Calculate the initial template which also allows the error calculation for the struck pixels
503 
504  for(int j=0; j<BXM2; ++j) {for(int i=0; i<BYM2; ++i) {template2d[j][i] = 0.f;}}
505  templ2D.xytemp(x2D, y2D, yd, xd, template2d, true, dpdx2d, qtemplate);
506 
507  float sumptdt1 = 0., sumptdt2 = 0.;
508  float sumdtdt11 = 0., sumdtdt12 = 0., sumdtdt22 = 0.;
509  chi2 = 0.f;
510  qactive = 0.f;
511  // Loop over all pixels and calculate matrices
512 
513  for(int k=0; k<tpixel; ++k) {
514  int jpix = indexxy[0][k];
515  int ipix = indexxy[1][k];
516  float qtpixel = template2d[jpix][ipix];
517  float pt = pixel[k]-qtpixel;
518  float dtdx = dpdx2d[0][jpix][ipix];
519  float dtdy = dpdx2d[1][jpix][ipix];
520  chi2 += pt*pt/sigma2[k];
521  sumptdt1 += pt*dtdx/sigma2[k];
522  sumptdt2 += pt*dtdy/sigma2[k];
523  sumdtdt11 += dtdx*dtdx/sigma2[k];
524  sumdtdt12 += dtdx*dtdy/sigma2[k];
525  sumdtdt22 += dtdy*dtdy/sigma2[k];
526  if(k < npixel) {qactive += qtpixel;}
527  }
528 
529  // invert the parameter covariance matrix
530 
531  float D = sumdtdt11*sumdtdt22 - sumdtdt12*sumdtdt12;
532 
533  // If the matrix is non-singular invert and solve
534 
535  if(std::abs(D) > 1.e-3) {
536 
537  minv11 = sumdtdt22/D;
538  minv12 = -sumdtdt12/D;
539  minv22 = sumdtdt11/D;
540 
541  qtfrac = qactive/qtemplate;
542 
543  //Calculate the step size
544 
545  xstep = minv11*sumptdt1 + minv12*sumptdt2;
546  ystep = minv12*sumptdt1 + minv22*sumptdt2;
547 
548  } else {
549 
550 // Assume alternately that ystep = 0 and then xstep = 0
551  if(sumdtdt11 > 0.0001f) {xstep = sumptdt1/sumdtdt11;} else {xstep = 0.f;}
552  if(sumdtdt22 > 0.0001f) {ystep = sumptdt2/sumdtdt22;} else {ystep = 0.f;}
553  }
554  xstep *= 0.9f;
555  ystep *= 0.9f;
556  if(std::abs(xstep) > 2.*xsize || std::abs(ystep) > 2.*ysize) break;
557  x2D += xstep;
558  y2D += ystep;
559  ++niter;
560  }
561  }
562 
563  ipass = 0;
564 
565  if(npass > 1) {
566  // two passes, take smaller chisqared
567  if(chi2min[1] < chi2min[0]) {ipass = 1;}
568  }
569 
570 
571 
572 
573  // Correct the charge ratio for missing pixels
574 
575  float fq;
576  if(qtfrac0[ipass] < 0.10f || qtfrac0[ipass] > 1.f) {qtfrac0[ipass] = 1.f;}
577  fq = fq0/qtfrac0[ipass];
578 
579  // printf("qtfrac0 = %f \n", qtfrac0);
580 
581  if(fq > fbin[0]) {
582  qbin=0;
583  } else {
584  if(fq > fbin[1]) {
585  qbin=1;
586  } else {
587  if(fq > fbin[2]) {
588  qbin=2;
589  } else {
590  qbin=3;
591  }
592  }
593  }
594 
595  // Get charge related quantities
596 
597  float scalex = templ2D.scalex(qbin);
598  float scaley = templ2D.scaley(qbin);
599  float offsetx = templ2D.offsetx(qbin);
600  float offsety = templ2D.offsety(qbin);
601 
602  // This 2D code has the origin (0,0) at the lower left edge of the input cluster
603  // That is now pixel [shiftx,shifty] and the template reco convention is the middle
604  // of that pixel, so we need to correct
605 
606  xrec = x2D0[ipass] - xpos[shiftx] - offsetx;
607  yrec = y2D0[ipass] - ypos[shifty] - offsety;
608  if(xerr2[ipass] > 0.f) {
609  sigmax = scalex*sqrt(xerr2[ipass]);
610  if(sigmax < 3.f) sigmax = 3.f;
611  } else {
612  sigmax = 10000.f;
613  }
614  if(yerr2[ipass] > 0.f) {
615  sigmay = scaley*sqrt(yerr2[ipass]);
616  if(sigmay < 3.f) sigmay = 3.f;
617  } else {
618  sigmay = 10000.f;
619  }
620  if(id > 0) {
621  // Do chi2 probability calculation
622  double meanxy = (double)(npixel*templ2D.chi2ppix());
623  double chi2scale = (double)templ2D.chi2scale();
624  if(meanxy < 0.01) {meanxy = 0.01;}
625  double hndof = meanxy/2.f;
626  double hchi2 = chi2scale*chi2min[ipass]/2.f;
627  probxy = (float)(1. - TMath::Gamma(hndof, hchi2));
628  // Do the charge probability
629  float mpv = templ2D.mpvvav();
630  float sigmaQ = templ2D.sigmavav();
631  float kappa = templ2D.kappavav();
632  float xvav = (qtotal/qtfrac0[ipass]-mpv)/sigmaQ;
633  float beta2 = 1.f;
634  // VVIObj is a private port of CERNLIB VVIDIS
635  VVIObjF vvidist(kappa, beta2, 1);
636  float prvav = vvidist.fcn(xvav);
637  probQ = 1.f - prvav;
638  } else {
639  probxy = chi2min[ipass];
640  npixels = npixel;
641  probQ = 0.f;
642  }
643 
644  return 0;
645 } // PixelTempReco2D
T2HYP1
#define T2HYP1
Definition: SiPixelTemplateDefs.h:41
SiPixelTemplate2D::xsize
float xsize()
pixel x-size (microns)
Definition: SiPixelTemplate2D.h:208
mps_fire.i
i
Definition: mps_fire.py:338
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:965
SiPixelTemplate2D::chi2ppix
float chi2ppix()
average chi^2 per struck pixel
Definition: SiPixelTemplate2D.h:193
SiPixelTemplate2D::xysigma2
void xysigma2(float qpixel, int index, float &xysig2)
Definition: SiPixelTemplate2D.cc:1518
jlist
int jlist[30]
Definition: TauolaWrapper.h:79
digitizers_cfi.pixel
pixel
Definition: digitizers_cfi.py:16
maxpix
static const int maxpix
Definition: SiPixelLorentzAngle.h:49
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
ysize
const Int_t ysize
Definition: trackSplitPlot.h:43
SiPixelTemplate2D::s50
float s50()
1/2 of the pixel threshold signal in adc units
Definition: SiPixelTemplate2D.h:184
SiPixelTemplate2D::offsetx
float offsetx(int i)
x-offset in 4 charge bins
Definition: SiPixelTemplate2D.h:188
BXM2
#define BXM2
Definition: SiPixelTemplateDefs.h:36
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
BYM2
#define BYM2
Definition: SiPixelTemplateDefs.h:28
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
SiPixelTemplate2D::scalex
float scalex(int i)
x-error scale factor in 4 charge bins
Definition: SiPixelTemplate2D.h:186
SiPixelTemplate2D::chi2scale
float chi2scale()
scale factor for chi^2 distribution
Definition: SiPixelTemplate2D.h:194
THYP1
#define THYP1
Definition: SiPixelTemplateDefs.h:23
SiPixelTemplate2D::lorxdrift
float lorxdrift()
signed lorentz x-width (microns)
Definition: SiPixelTemplate2D.h:201
SiPixelTemplate2D::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
Definition: SiPixelTemplate2D.cc:614
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
SiPixelTemplate2D::mpvvav
float mpvvav()
most probable Q in Vavilov distribution
Definition: SiPixelTemplate2D.h:197
funct::D
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:152
SiPixelTemplate2D::kappavav
float kappavav()
kappa parameter in Vavilov distribution
Definition: SiPixelTemplate2D.h:199
SiPixelTemplate2D::pixmax
float pixmax()
maximum pixel charge
Definition: SiPixelTemplate2D.h:182
SiPixelTemplate2D::lorydrift
float lorydrift()
signed lorentz y-width (microns)
Definition: SiPixelTemplate2D.h:200
SiPixelTemplate2D::sxymax
float sxymax()
max pixel signal for pixel error calculation
Definition: SiPixelTemplate2D.h:185
SiPixelTemplate2D::ysize
float ysize()
pixel y-size (microns)
Definition: SiPixelTemplate2D.h:209
SiPixelTemplate2D::qscale
float qscale()
charge scaling factor
Definition: SiPixelTemplate2D.h:183
SiPixelTemplate2D::qavg
float qavg()
average cluster charge for this set of track angles
Definition: SiPixelTemplate2D.h:181
SiPixelTemplate2D::fbin
float fbin(int i)
Return lower bound of Qbin definition.
Definition: SiPixelTemplate2D.h:190
Gamma
dbl * Gamma
Definition: mlp_gen.cc:38
kappa
static const G4double kappa
Definition: UrbanMscModel93.cc:35
gather_cfg.cout
cout
Definition: gather_cfg.py:144
SiPixelTemplate2D::offsety
float offsety(int i)
y-offset in 4 charge bins
Definition: SiPixelTemplate2D.h:189
vertices_cff.chi2
chi2
Definition: vertices_cff.py:34
xsize
const Int_t xsize
Definition: trackSplitPlot.h:42
SiPixelTemplate2D::clsleny
float clsleny()
cluster y-size
Definition: SiPixelTemplate2D.h:202
THXP1
#define THXP1
Definition: SiPixelTemplateDefs.h:32
SiPixelTemplate2D::sigmavav
float sigmavav()
scale factor in Vavilov distribution
Definition: SiPixelTemplate2D.h:198
SiPixelTemplate2D::scaley
float scaley(int i)
y-error scale factor in 4 charge bins
Definition: SiPixelTemplate2D.h:187
int SiPixelTemplateReco2D::PixelTempReco2D ( int  id,
float  cotalpha,
float  cotbeta,
float  locBz,
float  locBx,
int  edgeflagy,
int  edgeflagx,
ClusMatrix &  cluster,
SiPixelTemplate2D &  templ,
float &  yrec,
float &  sigmay,
float &  xrec,
float &  sigmax,
float &  probxy,
float &  probQ,
int &  qbin,
float &  deltay 
)

Definition at line 647 of file SiPixelTemplateReco2D.cc.

References PixelTempReco2D().

651 {
652  // Local variables
653  int npixels;
654  return SiPixelTemplateReco2D::PixelTempReco2D(id, cotalpha, cotbeta, locBz, locBx, edgeflagy, edgeflagx, cluster,
655  templ2D, yrec, sigmay, xrec, sigmax, probxy, probQ, qbin, deltay, npixels);
656 } // PixelTempReco2D
SiPixelTemplateReco2D::PixelTempReco2D
int PixelTempReco2D(int id, float cotalpha, float cotbeta, float locBz, float locBx, int edgeflagy, int edgeflagx, ClusMatrix &cluster, SiPixelTemplate2D &templ, float &yrec, float &sigmay, float &xrec, float &sigmax, float &probxy, float &probQ, int &qbin, float &deltay, int &npixel)
Definition: SiPixelTemplateReco2D.cc:83
Generated for CMSSW Reference Manual by doxygen 1.8.11