24 #ifdef SI_PIXEL_TEMPLATE_STANDALONE 33 #include "Math/DistFunc.h" 35 #ifndef SI_PIXEL_TEMPLATE_STANDALONE 39 #define LOGERROR(x) edm::LogError(x) 40 #define LOGDEBUG(x) LogDebug(x) 46 #define LOGERROR(x) std::cout << x << ": " 47 #define LOGDEBUG(x) std::cout << x << ": " 48 #define ENDL std::endl 85 float& xrec,
float& sigmax,
float& probxy,
float& probQ,
int& qbin,
float& deltay,
int& npixels)
93 const float fracpix = 0.45f;
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};
102 if(
id > 0) {
if(!templ2D.
interpolate(
id, cotalpha, cotbeta, locBz, locBx))
return 4;}
108 for(
int i=0;
i<3; ++
i) {fbin[
i] = templ2D.
fbin(
i);}
110 float q50 = templ2D.
s50();
111 float pseudopix = 0.2f*q50;
112 float pseudopix2 = q50*q50;
117 float qscale = templ2D.
qscale();
121 int nclusx = cluster.
mrow;
122 int nclusy = (
int)cluster.
mcol;
123 bool * xdouble = cluster.
xdouble;
124 bool * ydouble = cluster.
ydouble;
128 int imin=
BYM2, imax=0, jmin=
BXM2, jmax=0;
129 for(
int k=0;
k<nclusx*nclusy; ++
k) {
133 int i =
k - j*nclusy;
135 if(j < jmin) {jmin = j;}
136 if(j > jmax) {jmax = j;}
137 if(i < imin) {imin =
i;}
138 if(i > imax) {imax =
i;}
144 int shiftx =
THXP1 - (jmin+jmax)/2;
145 int shifty =
THYP1 - (imin+imax)/2;
147 if(shiftx < 1) shiftx = 1;
148 if(shifty < 1) shifty = 1;
157 float fq0 = qtotal/templ2D.
qavg();
163 for(
int j=0; j<
BXM2; ++j) {
for(
int i=0;
i<
BYM2; ++
i) {clusxy[j][
i] = 0.f;}}
165 const unsigned int NPIXMAX = 200;
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;
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) {
183 if(iy > -1 && iy < BYM2) yd[iy] =
true;
186 for(
int j=0; j<nclusx; ++j) {
189 if(jx > -1 && jx < BXM2) xd[jx] =
true;
196 float ypitch =
ysize;
197 if(yd[
i]) {ypitch +=
ysize;}
198 ye[i+1] = ye[
i] + ypitch;
199 ypos[
i] = ye[
i] + ypitch/2.;
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.;
208 for(
int i=0;
i<nclusy; ++
i) {
211 if(ydouble[
i]) {maxpix *=2.f;}
212 if(iy > -1 && iy < BYM2) {
213 for(
int j=0; j<nclusx; ++j) {
215 if(jx > -1 && jx < BXM2) {
216 if(cluster(j,i) > maxpix) {
219 clusxy[jx][iy] = cluster(j,i);
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];
235 #ifndef SI_PIXEL_TEMPLATE_STANDALONE 236 throw cms::Exception(
"DataCorrupt") <<
"PixelTemplateReco2D::number of pixels above threshold = " << npixel << std::endl;
238 std::cout <<
"PixelTemplateReco2D::number of pixels above threshold = " << npixel << std::endl;
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]);
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]);
268 if(imisslow < 0) imisslow =
i;
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]);
275 if(ylow < 0.
f || yhigh < 0.
f) {
280 float templeny = templ2D.
clsleny();
281 deltay = templeny - (yhigh - ylow)/ysize;
286 float x0 = 0.5f*(xlow0 + xhigh0) - templ2D.
lorxdrift();
287 float y0 = 0.5f*(ylow + yhigh) - templ2D.
lorydrift();
314 #ifndef SI_PIXEL_TEMPLATE_STANDALONE 315 throw cms::Exception(
"DataCorrupt") <<
"PixelTemplateReco2D::illegal edgeflagy = " << edgeflagy << std::endl;
317 std::cout <<
"PixelTemplate:2D:illegal edgeflagy = " << edgeflagy << std::endl;
333 #ifndef SI_PIXEL_TEMPLATE_STANDALONE 334 throw cms::Exception(
"DataCorrupt") <<
"PixelTemplateReco2D::illegal edgeflagx = " << edgeflagx << std::endl;
336 std::cout <<
"PixelTemplate:2D:illegal edgeflagx = " << edgeflagx << std::endl;
343 float chi2min[2], xerr2[2], yerr2[2];
344 float x2D0[2], y2D0[2], qtfrac0[2];
348 for(ipass = 0; ipass < npass; ++ipass) {
357 for(
int k=npixel;
k<tpixel; ++
k) {
358 int j = indexxy[0][
k];
359 int i = indexxy[1][
k];
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;
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;
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;
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;
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;
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;
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;
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;
457 chi2min[ipass] = 1000000.f;
458 float chi2, qtemplate, qactive, qtfrac = 0.f, x2D = 0.f, y2D = 0.f;
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;
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;}
478 if(chi2 < chi2min[ipass]) {
479 chi2min[ipass] =
chi2;
482 qtfrac = qactive/qtemplate;
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))) {
496 qtfrac0[ipass] = qtfrac;
497 xerr2[ipass] = minv11;
498 yerr2[ipass] = minv22;
499 chi2min[ipass] =
chi2;
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);
507 float sumptdt1 = 0., sumptdt2 = 0.;
508 float sumdtdt11 = 0., sumdtdt12 = 0., sumdtdt22 = 0.;
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;}
531 float D = sumdtdt11*sumdtdt22 - sumdtdt12*sumdtdt12;
537 minv11 = sumdtdt22/
D;
538 minv12 = -sumdtdt12/
D;
539 minv22 = sumdtdt11/
D;
541 qtfrac = qactive/qtemplate;
545 xstep = minv11*sumptdt1 + minv12*sumptdt2;
546 ystep = minv12*sumptdt1 + minv22*sumptdt2;
551 if(sumdtdt11 > 0.0001
f) {xstep = sumptdt1/sumdtdt11;}
else {xstep = 0.f;}
552 if(sumdtdt22 > 0.0001
f) {ystep = sumptdt2/sumdtdt22;}
else {ystep = 0.f;}
567 if(chi2min[1] < chi2min[0]) {ipass = 1;}
576 if(qtfrac0[ipass] < 0.10
f || qtfrac0[ipass] > 1.
f) {qtfrac0[ipass] = 1.f;}
577 fq = fq0/qtfrac0[ipass];
597 float scalex = templ2D.
scalex(qbin);
598 float scaley = templ2D.
scaley(qbin);
599 float offsetx = templ2D.
offsetx(qbin);
600 float offsety = templ2D.
offsety(qbin);
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;
614 if(yerr2[ipass] > 0.
f) {
615 sigmay = scaley*
sqrt(yerr2[ipass]);
616 if(sigmay < 3.
f) sigmay = 3.f;
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;
629 float mpv = templ2D.
mpvvav();
632 float xvav = (qtotal/qtfrac0[ipass]-mpv)/sigmaQ;
635 VVIObjF vvidist(kappa, beta2, 1);
636 float prvav = vvidist.
fcn(xvav);
639 probxy = chi2min[ipass];
649 float& probxy,
float& probQ,
int& qbin,
float& deltay)
655 templ2D, yrec, sigmay, xrec, sigmax, probxy, probQ, qbin, deltay, npixels);
float xsize()
pixel x-size (microns)
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)
float chi2ppix()
average chi^2 per struck pixel
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)
void xysigma2(float qpixel, int index, float &xysig2)
float s50()
1/2 of the pixel threshold signal in adc units
float offsetx(int i)
x-offset in 4 charge bins
Abs< T >::type abs(const T &t)
float scalex(int i)
x-error scale factor in 4 charge bins
float chi2scale()
scale factor for chi^2 distribution
float lorxdrift()
signed lorentz x-width (microns)
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
float mpvvav()
most probable Q in Vavilov distribution
DecomposeProduct< arg, typename Div::arg > D
float kappavav()
kappa parameter in Vavilov distribution
float pixmax()
maximum pixel charge
float lorydrift()
signed lorentz y-width (microns)
float sxymax()
max pixel signal for pixel error calculation
float ysize()
pixel y-size (microns)
float qscale()
charge scaling factor
float qavg()
average cluster charge for this set of track angles
float fbin(int i)
Return lower bound of Qbin definition.
static const G4double kappa
float offsety(int i)
y-offset in 4 charge bins
float clsleny()
cluster y-size
float sigmavav()
scale factor in Vavilov distribution
float scaley(int i)
y-error scale factor in 4 charge bins