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 81 of file SiPixelTemplateReco2D.cc.

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

85 {
86  // Local variables
87 
88  float template2d[BXM2][BYM2], dpdx2d[2][BXM2][BYM2], fbin[3];
89 
90  // fraction of truncation signal to measure the cluster ends
91  const float fracpix = 0.45f;
92 
93  const int nilist = 9, njlist = 5;
94  const float ilist[nilist] = {0.f, -1.f, -0.75f, -0.5f, -0.25f, 0.25f, 0.5f, 0.75f, 1.f};
95  const float jlist[njlist] = {0.f, -0.5f, -0.25f, 0.25f, 0.5f};
96 
97 
98  // Extract some relevant info from the 2D template
99 
100  if(id > 0) {if(!templ2D.interpolate(id, cotalpha, cotbeta, locBz, locBx)) return 4;}
101  float xsize = templ2D.xsize();
102  float ysize = templ2D.ysize();
103 
104  // Allow Qbin Q/Q_avg fractions to vary to optimize error estimation
105 
106  for(int i=0; i<3; ++i) {fbin[i] = templ2D.fbin(i);}
107 
108  float q50 = templ2D.s50();
109  float pseudopix = 0.2f*q50;
110  float pseudopix2 = q50*q50;
111 
112 
113  // Get charge scaling factor
114 
115  float qscale = templ2D.qscale();
116 
117  // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
118 
119  int nclusx = cluster.mrow;
120  int nclusy = (int)cluster.mcol;
121  bool * xdouble = cluster.xdouble;
122  bool * ydouble = cluster.ydouble;
123 
124  // First, rescale all pixel charges and compute total charge
125  float qtotal = 0.f;
126  int imin=BYM2, imax=0, jmin=BXM2, jmax=0;
127  for(int k=0; k<nclusx*nclusy; ++k) {
128  cluster.matrix[k] *= qscale;
129  qtotal +=cluster.matrix[k];
130  int j = k/nclusy;
131  int i = k - j*nclusy;
132  if(cluster.matrix[k] > 0.f) {
133  if(j < jmin) {jmin = j;}
134  if(j > jmax) {jmax = j;}
135  if(i < imin) {imin = i;}
136  if(i > imax) {imax = i;}
137  }
138  }
139 
140 // Calculate the shifts needed to center the cluster in the buffer
141 
142  int shiftx = THXP1 - (jmin+jmax)/2;
143  int shifty = THYP1 - (imin+imax)/2;
144 // Always shift by at least one pixel
145  if(shiftx < 1) shiftx = 1;
146  if(shifty < 1) shifty = 1;
147 // Shift the min maxes too
148  jmin += shiftx;
149  jmax += shiftx;
150  imin += shifty;
151  imax += shifty;
152 
153  // uncertainty and final corrections depend upon total charge bin
154 
155  float fq0 = qtotal/templ2D.qavg();
156 
157  // Next, copy the cluster and double pixel flags into a shifted workspace to
158  // allow for significant zeros [pseudopixels] to be added
159 
160  float clusxy[BXM2][BYM2];
161  for(int j=0; j<BXM2; ++j) {for(int i=0; i<BYM2; ++i) {clusxy[j][i] = 0.f;}}
162 
163  const unsigned int NPIXMAX = 200;
164 
165  int indexxy[2][NPIXMAX];
166  float pixel[NPIXMAX];
167  float sigma2[NPIXMAX];
168  float minmax = templ2D.pixmax();
169  float ylow0 = 0.f, yhigh0 = 0.f, xlow0 = 0.f, xhigh0 = 0.f;
170  int npixel = 0;
171  float ysum[BYM2], ye[BYM2+1], ypos[BYM2], xpos[BXM2], xe[BXM2+1];
172  bool yd[BYM2], xd[BXM2];
173 
174 // Copy double pixel flags to shifted arrays
175 
176  for(int i=0; i<BYM2; ++i) {ysum[i] = 0.f; yd[i] = false;}
177  for(int j=0; j<BXM2; ++j) {xd[j] = false;}
178  for(int i=0; i<nclusy; ++i) {
179  if(ydouble[i]) {
180  int iy = i+shifty;
181  if(iy > -1 && iy < BYM2) yd[iy] = true;
182  }
183  }
184  for(int j=0; j<nclusx; ++j) {
185  if(xdouble[j]) {
186  int jx = j+shiftx;
187  if(jx > -1 && jx < BXM2) xd[jx] = true;
188  }
189  }
190 // Work out the positions of the rows+columns relative to the lower left edge of pixel[1,1]
191  xe[0] = -xsize;
192  ye[0] = -ysize;
193  for(int i=0; i<BYM2; ++i) {
194  float ypitch = ysize;
195  if(yd[i]) {ypitch += ysize;}
196  ye[i+1] = ye[i] + ypitch;
197  ypos[i] = ye[i] + ypitch/2.;
198  }
199  for(int j=0; j<BXM2; ++j) {
200  float xpitch = xsize;
201  if(xd[j]) {xpitch += xsize;}
202  xe[j+1] = xe[j] + xpitch;
203  xpos[j] = xe[j] + xpitch/2.;
204  }
205 // Shift the cluster center to the central pixel of the array, truncate big signals
206  for(int i=0; i<nclusy; ++i) {
207  int iy = i+shifty;
208  float maxpix = minmax;
209  if(ydouble[i]) {maxpix *=2.f;}
210  if(iy > -1 && iy < BYM2) {
211  for(int j=0; j<nclusx; ++j) {
212  int jx = j+shiftx;
213  if(jx > -1 && jx < BXM2) {
214  if(cluster(j,i) > maxpix) {
215  clusxy[jx][iy] = maxpix;
216  } else {
217  clusxy[jx][iy] = cluster(j,i);
218  }
219  if(clusxy[jx][iy] > 0.f) {
220  ysum[iy] += clusxy[jx][iy];
221  indexxy[0][npixel] = jx;
222  indexxy[1][npixel] = iy;
223  pixel[npixel] = clusxy[jx][iy];
224  ++npixel;
225  }
226  }
227  }
228  }
229  }
230 // Get the shifted coordinates of the cluster ends
231  xlow0 = xe[jmin];
232  xhigh0 = xe[jmax+1];
233  ylow0 = ye[imin];
234  yhigh0 = ye[imax+1];
235 
236  // Next, calculate the error^2 [need to know which is the middle y pixel of the cluster]
237 
238  int ypixoff = T2HYP1 - (imin+imax)/2;
239  for(int k=0; k<npixel; ++k) {
240  int ypixeff = ypixoff + indexxy[1][k];
241  templ2D.xysigma2(pixel[k], ypixeff, sigma2[k]);
242  }
243 
244  // Next, find the half-height edges of the y-projection and identify any
245  // missing columns to remove from the fit
246 
247  int imisslow = -1, imisshigh = -1, jmisslow = -1, jmisshigh = -1;
248  float ylow = -1.f, yhigh = -1.f;
249  float hmaxpix = fracpix*templ2D.sxymax();
250  for(int i=imin; i<=imax; ++i) {
251  if(ysum[i] > hmaxpix && ysum[i-1] < hmaxpix && ylow < 0.f) {
252  ylow = ypos[i-1] + (ypos[i]-ypos[i-1])*(hmaxpix-ysum[i-1])/(ysum[i]-ysum[i-1]);
253  }
254  if(ysum[i] < q50) {
255  if(imisslow < 0) imisslow = i;
256  imisshigh = i;
257  }
258  if(ysum[i] > hmaxpix && ysum[i+1] < hmaxpix) {
259  yhigh = ypos[i] + (ypos[i+1]-ypos[i])*(ysum[i]-hmaxpix)/(ysum[i]-ysum[i+1]);
260  }
261  }
262  if(ylow < 0.f || yhigh < 0.f) {
263  ylow = ylow0;
264  yhigh = yhigh0;
265  }
266 
267  float templeny = templ2D.clsleny();
268  deltay = templeny - (yhigh - ylow)/ysize;
269 
270 
271  // x0 and y0 are best guess seeds for the fit
272 
273  float x0 = 0.5f*(xlow0 + xhigh0) - templ2D.lorxdrift();
274  float y0 = 0.5f*(ylow + yhigh) - templ2D.lorydrift();
275 // float y1 = yhigh - halfy - templ2D.lorydrift();
276 // printf("y0 = %f, y1 = %f \n", y0, y1);
277  // float y0 = 0.5f*(ylow + yhigh);
278 
279  // If there are missing edge columns, set up missing column flags and number
280  // of minimization passes
281 
282  int npass = 1;
283 
284  switch(edgeflagy) {
285  case 0:
286  break;
287  case 1:
288  imisshigh = imin-1;
289  imisslow = -1;
290  break;
291  case 2:
292  imisshigh = -1;
293  imisslow = imax+1;
294  break;
295  case 3:
296  imisshigh = imin-1;
297  imisslow = -1;
298  npass = 2;
299  break;
300  default:
301 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
302  throw cms::Exception("DataCorrupt") << "PixelTemplateReco2D::illegal edgeflagy = " << edgeflagy << std::endl;
303 #else
304  std::cout << "PixelTemplate:2D:illegal edgeflagy = " << edgeflagy << std::endl;
305 #endif
306  }
307 
308  switch(edgeflagx) {
309  case 0:
310  break;
311  case 1:
312  jmisshigh = jmin-1;
313  jmisslow = -1;
314  break;
315  case 2:
316  jmisshigh = -1;
317  jmisslow = jmax+1;
318  break;
319  default:
320 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
321  throw cms::Exception("DataCorrupt") << "PixelTemplateReco2D::illegal edgeflagx = " << edgeflagx << std::endl;
322 #else
323  std::cout << "PixelTemplate:2D:illegal edgeflagx = " << edgeflagx << std::endl;
324 #endif
325  }
326 
327  // Define quantities to be saved for each pass
328 
329 
330  float chi2min[2], xerr2[2], yerr2[2];
331  float x2D0[2], y2D0[2], qtfrac0[2];
332  int ipass, tpixel;
333  // niter0[ipass] = niter; // Not used -- comment out for now.
334 
335  for(ipass = 0; ipass < npass; ++ipass) {
336 
337  if(ipass == 1) {
338 
339  // Now try again if both edges are possible
340 
341  imisshigh = -1;
342  imisslow = imax+1;
343 // erase pseudo pixels from previous pass
344  for(int k=npixel; k<tpixel; ++k) {
345  int j = indexxy[0][k];
346  int i = indexxy[1][k];
347  clusxy[j][i] = 0.f;
348  }
349  }
350 
351  // Next, add pseudo pixels around the periphery of the cluster
352 
353  tpixel = npixel;
354  for(int k=0; k<npixel; ++k) {
355  int j = indexxy[0][k];
356  int i = indexxy[1][k];
357  if((j-1) != jmisshigh) {
358  if(clusxy[j-1][i] < pseudopix) {
359  indexxy[0][tpixel] = j-1;
360  indexxy[1][tpixel] = i;
361  clusxy[j-1][i] = pseudopix;
362  pixel[tpixel] = pseudopix;
363  sigma2[tpixel] = pseudopix2;
364  ++tpixel;
365  }
366  }
367  if((j+1) != jmisslow) {
368  if(clusxy[j+1][i] < pseudopix) {
369  indexxy[0][tpixel] = j+1;
370  indexxy[1][tpixel] = i;
371  clusxy[j+1][i] = pseudopix;
372  pixel[tpixel] = pseudopix;
373  sigma2[tpixel] = pseudopix2;
374  ++tpixel;
375  }
376  }
377  // Don't add them if this is a dead column
378  if((i+1) != imisslow) {
379  if((j-1) != jmisshigh) {
380  if(clusxy[j-1][i+1] < pseudopix) {
381  indexxy[0][tpixel] = j-1;
382  indexxy[1][tpixel] = i+1;
383  clusxy[j-1][i+1] = pseudopix;
384  pixel[tpixel] = pseudopix;
385  sigma2[tpixel] = pseudopix2;
386  ++tpixel;
387  }
388  }
389  if(clusxy[j][i+1] < pseudopix) {
390  indexxy[0][tpixel] = j;
391  indexxy[1][tpixel] = i+1;
392  clusxy[j][i+1] = pseudopix;
393  pixel[tpixel] = pseudopix;
394  sigma2[tpixel] = pseudopix2;
395  ++tpixel;
396  }
397  if((j+1) != jmisslow) {
398  if(clusxy[j+1][i+1] < pseudopix) {
399  indexxy[0][tpixel] = j+1;
400  indexxy[1][tpixel] = i+1;
401  clusxy[j+1][i+1] = pseudopix;
402  pixel[tpixel] = pseudopix;
403  sigma2[tpixel] = pseudopix2;
404  ++tpixel;
405  }
406  }
407  }
408  // Don't add them if this is a dead column
409  if((i-1) != imisshigh) {
410  if((j-1) != jmisshigh) {
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  if(clusxy[j][i-1] < pseudopix) {
421  indexxy[0][tpixel] = j;
422  indexxy[1][tpixel] = i-1;
423  clusxy[j][i-1] = pseudopix;
424  pixel[tpixel] = pseudopix;
425  sigma2[tpixel] = pseudopix2;
426  ++tpixel;
427  }
428  if((j+1) != jmisslow) {
429  if(clusxy[j+1][i-1] < pseudopix) {
430  indexxy[0][tpixel] = j+1;
431  indexxy[1][tpixel] = i-1;
432  clusxy[j+1][i-1] = pseudopix;
433  pixel[tpixel] = pseudopix;
434  sigma2[tpixel] = pseudopix2;
435  ++tpixel;
436  }
437  }
438  }
439  }
440 
441 
442  // Calculate chi2 over a grid of several seeds and choose the smallest
443 
444  chi2min[ipass] = 1000000.f;
445  float chi2, qtemplate, qactive, qtfrac = 0.f, x2D = 0.f, y2D = 0.f;
446 
447  for(int is = 0; is<nilist; ++is) {
448  for(int js = 0; js<njlist; ++js) {
449  float xtry = x0 + jlist[js]*xsize;
450  float ytry = y0 + ilist[is]*ysize;
451  chi2 = 0.f;
452  qactive = 0.f;
453  for(int j=0; j<BXM2; ++j) {for(int i=0; i<BYM2; ++i) {template2d[j][i] = 0.f;}}
454  templ2D.xytemp(xtry, ytry, yd, xd, template2d, false, dpdx2d, qtemplate);
455  for(int k=0; k<tpixel; ++k) {
456  int jpix = indexxy[0][k];
457  int ipix = indexxy[1][k];
458  float qtpixel = template2d[jpix][ipix];
459  float pt = pixel[k]-qtpixel;
460  chi2 += pt*pt/sigma2[k];
461  if(k < npixel) {qactive += qtpixel;}
462  }
463  if(chi2 < chi2min[ipass]) {
464  chi2min[ipass] = chi2;
465  x2D = xtry;
466  y2D = ytry;
467  qtfrac = qactive/qtemplate;
468  }
469  }
470  }
471 
472  int niter = 0;
473  float xstep = 1.0f, ystep = 1.0f;
474  float minv11 = 1000.f, minv12 = 1000.f, minv22 = 1000.f;
475  chi2 = chi2min[ipass];
476  // int niter0[2];
477  while(chi2 <= chi2min[ipass] && niter < 15 && (niter < 2 || (fabs(xstep) > 0.2 && fabs(ystep) > 0.2))) {
478 
479  // Remember the present parameters
480  x2D0[ipass] = x2D;
481  y2D0[ipass] = y2D;
482  qtfrac0[ipass] = qtfrac;
483  xerr2[ipass] = minv11;
484  yerr2[ipass] = minv22;
485  chi2min[ipass] = chi2;
486 
487  // Calculate the initial template which also allows the error calculation for the struck pixels
488 
489  for(int j=0; j<BXM2; ++j) {for(int i=0; i<BYM2; ++i) {template2d[j][i] = 0.f;}}
490  templ2D.xytemp(x2D, y2D, yd, xd, template2d, true, dpdx2d, qtemplate);
491 
492  float sumptdt1 = 0., sumptdt2 = 0.;
493  float sumdtdt11 = 0., sumdtdt12 = 0., sumdtdt22 = 0.;
494  chi2 = 0.f;
495  qactive = 0.f;
496  // Loop over all pixels and calculate matrices
497 
498  for(int k=0; k<tpixel; ++k) {
499  int jpix = indexxy[0][k];
500  int ipix = indexxy[1][k];
501  float qtpixel = template2d[jpix][ipix];
502  float pt = pixel[k]-qtpixel;
503  float dtdx = dpdx2d[0][jpix][ipix];
504  float dtdy = dpdx2d[1][jpix][ipix];
505  chi2 += pt*pt/sigma2[k];
506  sumptdt1 += pt*dtdx/sigma2[k];
507  sumptdt2 += pt*dtdy/sigma2[k];
508  sumdtdt11 += dtdx*dtdx/sigma2[k];
509  sumdtdt12 += dtdx*dtdy/sigma2[k];
510  sumdtdt22 += dtdy*dtdy/sigma2[k];
511  if(k < npixel) {qactive += qtpixel;}
512  }
513 
514  // invert the parameter covariance matrix
515 
516  float D = sumdtdt11*sumdtdt22 - sumdtdt12*sumdtdt12;
517 
518  // If the matrix is non-singular invert and solve
519 
520  if(fabs(D) > 1.e-3) {
521 
522  minv11 = sumdtdt22/D;
523  minv12 = -sumdtdt12/D;
524  minv22 = sumdtdt11/D;
525 
526  qtfrac = qactive/qtemplate;
527 
528  //Calculate the step size
529 
530  xstep = minv11*sumptdt1 + minv12*sumptdt2;
531  ystep = minv12*sumptdt1 + minv22*sumptdt2;
532 
533  } else {
534 // Assume alternately that ystep = 0 and then xstep = 0
535  xstep = sumptdt1/sumdtdt11;
536  ystep = sumptdt2/sumdtdt22;
537  }
538  xstep *= 0.9f;
539  ystep *= 0.9f;
540  if(fabs(xstep) > 2.*xsize || fabs(ystep) > 2.*ysize) break;
541  x2D += xstep;
542  y2D += ystep;
543  ++niter;
544  }
545  }
546 
547  ipass = 0;
548 
549 //--- Somewhat experimental, keep this commented out:
550 // if(npass == 1) {
551 // // one pass, require that it have iterated
552 // if(niter0[0] == 0) {return 2;}
553 // } else {
554 // // two passes
555 // if(niter0[0] == 0 && niter0[1] == 0) {return 2;}
556 // if(niter0[0] > 0 && niter0[1] > 0) {
557 // // if both have iterated, take the smaller chi2
558 // if(chi2min[1] < chi2min[0]) {ipass = 1;}
559 // } else {
560 // // if one has iterated, take it
561 // if(niter0[1] > 0) {ipass = 1;}
562 // }
563 // }
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:329
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:958
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:1511
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
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:613
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:151
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:33
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:81
Generated for CMSSW Reference Manual by doxygen 1.8.11