File indexing completed on 2025-06-25 05:16:14
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0057
0058 #else
0059 #include <math.h>
0060 #endif
0061 #include <algorithm>
0062 #include <vector>
0063 #include <utility>
0064 #include <iostream>
0065
0066 #include "TMath.h"
0067 #include "Math/DistFunc.h"
0068
0069
0070
0071 static const int theVerboseLevel = 2;
0072 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0073 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateReco.h"
0074 #include "RecoLocalTracker/SiPixelRecHits/interface/VVIObjF.h"
0075 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0076 #define LOGERROR(x) edm::LogError(x)
0077 #define LOGDEBUG(x) LogDebug(x)
0078 #define ENDL " "
0079 #include "FWCore/Utilities/interface/Exception.h"
0080 #else
0081 #include "SiPixelTemplateReco.h"
0082 #include "VVIObjF.h"
0083 #define LOGERROR(x) std::cout << x << ": "
0084 #define LOGDEBUG(x) std::cout << x << ": "
0085 #define ENDL std::endl
0086 #endif
0087
0088 using namespace SiPixelTemplateReco;
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138 int SiPixelTemplateReco::PixelTempReco1D(int id,
0139 float cotalpha,
0140 float cotbeta,
0141 float locBz,
0142 float locBx,
0143 ClusMatrix& cluster,
0144 SiPixelTemplate& templ,
0145 float& yrec,
0146 float& sigmay,
0147 float& proby,
0148 float& xrec,
0149 float& sigmax,
0150 float& probx,
0151 int& qbin,
0152 int speed,
0153 bool deadpix,
0154 std::vector<std::pair<int, int> >& zeropix,
0155 float& probQ,
0156 int& nypix,
0157 int& nxpix,
0158 bool goodEdgeAlgo)
0159
0160 {
0161
0162 int i, j, k, minbin, binl, binh, binq, midpix, fypix, lypix, logypx;
0163 int fxpix, lxpix, logxpx, shifty, shiftx, nyzero[TYSIZE];
0164 int nclusx, nclusy;
0165 int deltaj, jmin, jmax, fxbin, lxbin, fybin, lybin, djy, djx;
0166
0167 float sythr, sxthr, rnorm, delta, sigma, sigavg, pseudopix, qscale, q50;
0168 float ss2, ssa, sa2, ssba, saba, sba2, rat, fq, qtotal, qpixel, fbin[3];
0169 float originx, originy, qfy, qly, qfx, qlx, bias, maxpix, minmax;
0170 double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
0171 double hchi2, hndof, prvav, mpv, sigmaQ, kappa, xvav, beta2;
0172 float ytemp[41][BYSIZE], xtemp[41][BXSIZE], ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
0173 float chi2ybin[41], chi2xbin[41], ysig2[BYSIZE], xsig2[BXSIZE];
0174 float yw2[BYSIZE], xw2[BXSIZE], ysw[BYSIZE], xsw[BXSIZE];
0175 bool yd[BYSIZE], xd[BXSIZE], anyyd, anyxd, calc_probQ, use_VVIObj;
0176 float ysize, xsize;
0177 const float probmin = {1.110223e-16};
0178 const float probQmin = {1.e-5};
0179
0180
0181
0182 const double mean1pix = {0.100};
0183 #ifdef SI_PIXEL_TEMPLATE_STANDALONE
0184 const double chi21min = {0.};
0185 #else
0186 const double chi21min = {0.160};
0187 #endif
0188
0189
0190
0191
0192 if (id >= 0) {
0193 if (!templ.interpolate(id, cotalpha, cotbeta, locBz, locBx, goodEdgeAlgo)) {
0194 if (theVerboseLevel > 2) {
0195 LOGDEBUG("SiPixelTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha
0196 << ", cot(beta) = " << cotbeta << ", local B_z = " << locBz
0197 << ", local B_x = " << locBx << ", template ID = " << id
0198 << ", no reconstruction performed" << ENDL;
0199 }
0200 return 20;
0201 }
0202 }
0203
0204
0205
0206 calc_probQ = false;
0207 use_VVIObj = false;
0208 if (speed < 0) {
0209 calc_probQ = true;
0210 if (speed < -1)
0211 use_VVIObj = true;
0212 speed = 0;
0213 }
0214
0215 if (speed > 3) {
0216 calc_probQ = true;
0217 if (speed < 5)
0218 use_VVIObj = true;
0219 speed = 3;
0220 }
0221
0222
0223
0224 xsize = templ.xsize();
0225 ysize = templ.ysize();
0226
0227
0228
0229 for (i = 0; i < 3; ++i) {
0230 fbin[i] = templ.fbin(i);
0231 }
0232
0233
0234
0235 q50 = templ.s50();
0236 pseudopix = 0.2f * q50;
0237
0238
0239
0240 qscale = templ.qscale();
0241
0242
0243
0244 nclusx = cluster.mrow;
0245 nclusy = (int)cluster.mcol;
0246 auto const xdouble = cluster.xdouble;
0247 auto const ydouble = cluster.ydouble;
0248
0249
0250 qtotal = 0.f;
0251 for (i = 0; i < nclusx * nclusy; ++i) {
0252 cluster.matrix[i] *= qscale;
0253 qtotal += cluster.matrix[i];
0254 }
0255
0256 minmax = templ.pixmax();
0257 for (j = 0; j < nclusx; ++j)
0258 for (i = 0; i < nclusy; ++i) {
0259 maxpix = minmax;
0260 if (ydouble[i]) {
0261 maxpix *= 2.f;
0262 }
0263 if (cluster(j, i) > maxpix) {
0264 cluster(j, i) = maxpix;
0265 }
0266 }
0267
0268
0269
0270 if (deadpix) {
0271 fypix = BYM3;
0272 lypix = -1;
0273 memset(nyzero, 0, TYSIZE * sizeof(int));
0274 memset(ysum, 0, BYSIZE * sizeof(float));
0275 for (i = 0; i < nclusy; ++i) {
0276
0277 for (j = 0; j < nclusx; ++j) {
0278 ysum[i] += cluster(j, i);
0279 }
0280 if (ysum[i] > 0.f) {
0281
0282 if (i < fypix) {
0283 fypix = i;
0284 }
0285 if (i > lypix) {
0286 lypix = i;
0287 }
0288 }
0289 }
0290
0291
0292
0293
0294
0295 std::vector<std::pair<int, int> >::const_iterator zeroIter = zeropix.begin(), zeroEnd = zeropix.end();
0296 for (; zeroIter != zeroEnd; ++zeroIter) {
0297 i = zeroIter->second;
0298 if (i < 0 || i > TYSIZE - 1) {
0299 LOGERROR("SiPixelTemplateReco") << "dead pixel column y-index " << i << ", no reconstruction performed" << ENDL;
0300 return 11;
0301 }
0302
0303
0304 ++nyzero[i];
0305
0306 if (i < fypix) {
0307 fypix = i;
0308 }
0309 if (i > lypix) {
0310 lypix = i;
0311 }
0312 }
0313
0314 nypix = lypix - fypix + 1;
0315
0316
0317
0318 for (zeroIter = zeropix.begin(); zeroIter != zeroEnd; ++zeroIter) {
0319 i = zeroIter->second;
0320 j = zeroIter->first;
0321 if (j < 0 || j > TXSIZE - 1) {
0322 LOGERROR("SiPixelTemplateReco") << "dead pixel column x-index " << j << ", no reconstruction performed" << ENDL;
0323 return 12;
0324 }
0325 if ((i == fypix || i == lypix) && nypix > 1) {
0326 maxpix = templ.symax() / 2.;
0327 } else {
0328 maxpix = templ.symax();
0329 }
0330 if (ydouble[i]) {
0331 maxpix *= 2.;
0332 }
0333 if (nyzero[i] > 0 && nyzero[i] < 3) {
0334 qpixel = (maxpix - ysum[i]) / (float)nyzero[i];
0335 } else {
0336 qpixel = 1.;
0337 }
0338 if (qpixel < 1.) {
0339 qpixel = 1.;
0340 }
0341 cluster(j, i) = qpixel;
0342
0343 qtotal += qpixel;
0344 }
0345
0346 }
0347
0348
0349
0350 for (i = 0; i < BYSIZE; ++i) {
0351 ysum[i] = 0.f;
0352 yd[i] = false;
0353 }
0354 k = 0;
0355 anyyd = false;
0356 for (i = 0; i < nclusy; ++i) {
0357 for (j = 0; j < nclusx; ++j) {
0358 ysum[k] += cluster(j, i);
0359 }
0360
0361
0362
0363 if (ydouble[i]) {
0364 ysum[k] /= 2.f;
0365 ysum[k + 1] = ysum[k];
0366 yd[k] = true;
0367 yd[k + 1] = false;
0368 k = k + 2;
0369 anyyd = true;
0370 } else {
0371 yd[k] = false;
0372 ++k;
0373 }
0374 if (k > BYM1) {
0375 break;
0376 }
0377 }
0378
0379
0380
0381 for (i = 0; i < BXSIZE; ++i) {
0382 xsum[i] = 0.f;
0383 xd[i] = false;
0384 }
0385 k = 0;
0386 anyxd = false;
0387 for (j = 0; j < nclusx; ++j) {
0388 for (i = 0; i < nclusy; ++i) {
0389 xsum[k] += cluster(j, i);
0390 }
0391
0392
0393
0394 if (xdouble[j]) {
0395 xsum[k] /= 2.;
0396 xsum[k + 1] = xsum[k];
0397 xd[k] = true;
0398 xd[k + 1] = false;
0399 k = k + 2;
0400 anyxd = true;
0401 } else {
0402 xd[k] = false;
0403 ++k;
0404 }
0405 if (k > BXM1) {
0406 break;
0407 }
0408 }
0409
0410
0411
0412 fypix = -1;
0413 nypix = 0;
0414 lypix = 0;
0415 logypx = 0;
0416 for (i = 0; i < BYSIZE; ++i) {
0417 if (ysum[i] > 0.f) {
0418 if (fypix == -1) {
0419 fypix = i;
0420 }
0421 if (!yd[i]) {
0422 ysort[logypx] = ysum[i];
0423 ++logypx;
0424 }
0425 ++nypix;
0426 lypix = i;
0427 }
0428 }
0429
0430
0431
0432
0433
0434 if ((lypix - fypix + 1) != nypix || nypix == 0) {
0435 LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster doesn't agree with number of pixels above threshold"
0436 << ENDL;
0437 if (theVerboseLevel > 2) {
0438 LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
0439 for (i = 0; i < BYSIZE - 1; ++i) {
0440 LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";
0441 }
0442 LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE - 1] << ENDL;
0443 }
0444
0445 return 1;
0446 }
0447
0448
0449
0450 if (nypix > TYSIZE) {
0451 LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster is larger than maximum template size" << ENDL;
0452 if (theVerboseLevel > 2) {
0453 LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
0454 for (i = 0; i < BYSIZE - 1; ++i) {
0455 LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";
0456 }
0457 LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE - 1] << ENDL;
0458 }
0459
0460 return 6;
0461 }
0462
0463
0464
0465
0466
0467
0468
0469
0470 if (goodEdgeAlgo) {
0471 if (cotbeta >= 0) {
0472 midpix = (int)(fypix + 0.5 * cotbeta * templ.zsize() / ysize);
0473 } else {
0474 midpix = (int)(lypix + 0.5 * cotbeta * templ.zsize() / ysize);
0475 }
0476 shifty = BHY - midpix;
0477 }
0478 else {
0479 midpix = (fypix + lypix) / 2;
0480 shifty = templ.cytemp() - midpix;
0481 }
0482
0483
0484
0485 int lytmp = lypix + shifty;
0486 int fytmp = fypix + shifty;
0487
0488
0489
0490 if (fytmp <= 1) {
0491 return 8;
0492 }
0493 if (lytmp >= BYM2) {
0494 return 8;
0495 }
0496
0497
0498
0499 if (shifty > 0) {
0500 for (i = lypix; i >= fypix; --i) {
0501 ysum[i + shifty] = ysum[i];
0502 ysum[i] = 0.;
0503 yd[i + shifty] = yd[i];
0504 yd[i] = false;
0505 }
0506 } else if (shifty < 0) {
0507 for (i = fypix; i <= lypix; ++i) {
0508 ysum[i + shifty] = ysum[i];
0509 ysum[i] = 0.;
0510 yd[i + shifty] = yd[i];
0511 yd[i] = false;
0512 }
0513 }
0514
0515 lypix = lytmp;
0516 fypix = fytmp;
0517
0518
0519
0520 ysum[fypix - 1] = pseudopix;
0521 ysum[fypix - 2] = pseudopix;
0522 ysum[lypix + 1] = pseudopix;
0523 ysum[lypix + 2] = pseudopix;
0524
0525
0526
0527 if (ydouble[0]) {
0528 originy = -0.5f;
0529 } else {
0530 originy = 0.f;
0531 }
0532
0533
0534
0535 fxpix = -1;
0536 nxpix = 0;
0537 lxpix = 0;
0538 logxpx = 0;
0539 for (i = 0; i < BXSIZE; ++i) {
0540 if (xsum[i] > 0.) {
0541 if (fxpix == -1) {
0542 fxpix = i;
0543 }
0544 if (!xd[i]) {
0545 xsort[logxpx] = xsum[i];
0546 ++logxpx;
0547 }
0548 ++nxpix;
0549 lxpix = i;
0550 }
0551 }
0552
0553
0554
0555
0556
0557 if ((lxpix - fxpix + 1) != nxpix) {
0558 LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold"
0559 << ENDL;
0560 if (theVerboseLevel > 2) {
0561 LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
0562 for (i = 0; i < BXSIZE - 1; ++i) {
0563 LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";
0564 }
0565 LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE - 1] << ENDL;
0566 }
0567
0568 return 2;
0569 }
0570
0571
0572
0573 if (nxpix > TXSIZE) {
0574 LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
0575 if (theVerboseLevel > 2) {
0576 LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
0577 for (i = 0; i < BXSIZE - 1; ++i) {
0578 LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";
0579 }
0580 LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE - 1] << ENDL;
0581 }
0582
0583 return 7;
0584 }
0585
0586
0587
0588
0589
0590
0591
0592
0593 midpix = (fxpix + lxpix) / 2;
0594 shiftx = templ.cxtemp() - midpix;
0595
0596
0597
0598 int lxtmp = lxpix + shiftx;
0599 int fxtmp = fxpix + shiftx;
0600
0601
0602
0603 if (fxtmp <= 1) {
0604 return 9;
0605 }
0606 if (lxtmp >= BXM2) {
0607 return 9;
0608 }
0609
0610
0611
0612 if (shiftx > 0) {
0613 for (i = lxpix; i >= fxpix; --i) {
0614 xsum[i + shiftx] = xsum[i];
0615 xsum[i] = 0.;
0616 xd[i + shiftx] = xd[i];
0617 xd[i] = false;
0618 }
0619 } else if (shiftx < 0) {
0620 for (i = fxpix; i <= lxpix; ++i) {
0621 xsum[i + shiftx] = xsum[i];
0622 xsum[i] = 0.;
0623 xd[i + shiftx] = xd[i];
0624 xd[i] = false;
0625 }
0626 }
0627
0628 lxpix = lxtmp;
0629 fxpix = fxtmp;
0630
0631 xsum[fxpix - 1] = pseudopix;
0632 xsum[fxpix - 2] = pseudopix;
0633 xsum[lxpix + 1] = pseudopix;
0634 xsum[lxpix + 2] = pseudopix;
0635
0636
0637
0638 if (xdouble[0]) {
0639 originx = -0.5f;
0640 } else {
0641 originx = 0.f;
0642 }
0643
0644
0645
0646 fq = qtotal / templ.qavg();
0647 if (fq > fbin[0]) {
0648 binq = 0;
0649 } else {
0650 if (fq > fbin[1]) {
0651 binq = 1;
0652 } else {
0653 if (fq > fbin[2]) {
0654 binq = 2;
0655 } else {
0656 binq = 3;
0657 }
0658 }
0659 }
0660
0661
0662
0663 qbin = binq;
0664 if (!deadpix && qtotal < 0.95f * templ.qmin()) {
0665 qbin = 5;
0666 } else {
0667 if (!deadpix && qtotal < 0.95f * templ.qmin(1)) {
0668 qbin = 4;
0669 }
0670 }
0671 if (theVerboseLevel > 9) {
0672 LOGDEBUG("SiPixelTemplateReco") << "ID = " << id << " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta
0673 << " nclusx = " << nclusx << " nclusy = " << nclusy << ENDL;
0674 }
0675
0676
0677
0678
0679
0680 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0681 if (speed < 0 || speed > 3) {
0682 throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::PixelTempReco1D called with illegal speed = " << speed
0683 << std::endl;
0684 }
0685 #else
0686 assert(speed >= 0 && speed < 4);
0687 #endif
0688 fybin = 3;
0689 lybin = 37;
0690 fxbin = 3;
0691 lxbin = 37;
0692 djy = 1;
0693 djx = 1;
0694 if (speed > 0) {
0695 fybin = 8;
0696 lybin = 32;
0697 if (yd[fypix]) {
0698 fybin = 4;
0699 lybin = 36;
0700 }
0701 if (lypix > fypix) {
0702 if (yd[lypix - 1]) {
0703 fybin = 4;
0704 lybin = 36;
0705 }
0706 }
0707 fxbin = 8;
0708 lxbin = 32;
0709 if (xd[fxpix]) {
0710 fxbin = 4;
0711 lxbin = 36;
0712 }
0713 if (lxpix > fxpix) {
0714 if (xd[lxpix - 1]) {
0715 fxbin = 4;
0716 lxbin = 36;
0717 }
0718 }
0719 }
0720
0721 if (speed > 1) {
0722 djy = 2;
0723 djx = 2;
0724 if (speed > 2) {
0725 if (!anyyd) {
0726 djy = 4;
0727 }
0728 if (!anyxd) {
0729 djx = 4;
0730 }
0731 }
0732 }
0733
0734 if (theVerboseLevel > 9) {
0735 LOGDEBUG("SiPixelTemplateReco") << "fypix " << fypix << " lypix = " << lypix << " fybin = " << fybin
0736 << " lybin = " << lybin << " djy = " << djy << " logypx = " << logypx << ENDL;
0737 LOGDEBUG("SiPixelTemplateReco") << "fxpix " << fxpix << " lxpix = " << lxpix << " fxbin = " << fxbin
0738 << " lxbin = " << lxbin << " djx = " << djx << " logxpx = " << logxpx << ENDL;
0739 }
0740
0741
0742
0743 templ.ytemp(fybin, lybin, ytemp);
0744
0745 templ.xtemp(fxbin, lxbin, xtemp);
0746
0747
0748
0749
0750
0751
0752
0753 if (nypix > logypx) {
0754 i = fypix;
0755 while (i < lypix) {
0756 if (yd[i] && !yd[i + 1]) {
0757 for (j = fybin; j <= lybin; ++j) {
0758
0759
0760 sigavg = (ytemp[j][i] + ytemp[j][i + 1]) / 2.f;
0761 ytemp[j][i] = sigavg;
0762 ytemp[j][i + 1] = sigavg;
0763 }
0764 i += 2;
0765 } else {
0766 ++i;
0767 }
0768 }
0769 }
0770
0771
0772
0773 sythr = 1.1 * (templ.symax());
0774
0775
0776
0777 std::sort(&ysort[0], &ysort[logypx]);
0778 if (logypx == 1) {
0779 sythr = 1.01f * ysort[0];
0780 } else {
0781 if (ysort[1] > sythr) {
0782 sythr = 1.01f * ysort[1];
0783 }
0784 }
0785
0786
0787
0788
0789 templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
0790
0791
0792
0793 chi2ymin = 1.e15;
0794 for (i = fybin; i <= lybin; ++i) {
0795 chi2ybin[i] = -1.e15f;
0796 }
0797 ss2 = 0.f;
0798 for (i = fypix - 2; i <= lypix + 2; ++i) {
0799 yw2[i] = 1.f / ysig2[i];
0800 ysw[i] = ysum[i] * yw2[i];
0801 ss2 += ysum[i] * ysw[i];
0802 }
0803
0804 minbin = -1;
0805 deltaj = djy;
0806 jmin = fybin;
0807 jmax = lybin;
0808 while (deltaj > 0) {
0809 for (j = jmin; j <= jmax; j += deltaj) {
0810 if (chi2ybin[j] < -100.f) {
0811 ssa = 0.f;
0812 sa2 = 0.f;
0813 for (i = fypix - 2; i <= lypix + 2; ++i) {
0814 ssa += ysw[i] * ytemp[j][i];
0815 sa2 += ytemp[j][i] * ytemp[j][i] * yw2[i];
0816 }
0817 rat = ssa / ss2;
0818 if (rat <= 0.f) {
0819 LOGERROR("SiPixelTemplateReco") << "illegal chi2ymin normalization (1) = " << rat << ENDL;
0820 rat = 1.;
0821 }
0822 chi2ybin[j] = ss2 - 2.f * ssa / rat + sa2 / (rat * rat);
0823 }
0824 if (chi2ybin[j] < chi2ymin) {
0825 chi2ymin = chi2ybin[j];
0826 minbin = j;
0827 }
0828 }
0829 deltaj /= 2;
0830 if (minbin > fybin) {
0831 jmin = minbin - deltaj;
0832 } else {
0833 jmin = fybin;
0834 }
0835 if (minbin < lybin) {
0836 jmax = minbin + deltaj;
0837 } else {
0838 jmax = lybin;
0839 }
0840 }
0841
0842 if (theVerboseLevel > 9) {
0843 LOGDEBUG("SiPixelTemplateReco") << "minbin " << minbin << " chi2ymin = " << chi2ymin << ENDL;
0844 }
0845
0846
0847
0848 if (logypx == 1) {
0849 if (nypix == 1) {
0850 delta = templ.dyone();
0851 sigma = templ.syone();
0852 } else {
0853 delta = templ.dytwo();
0854 sigma = templ.sytwo();
0855 }
0856
0857 yrec = 0.5f * (fypix + lypix - 2 * shifty + 2.f * originy) * ysize - delta;
0858 if (sigma <= 0.f) {
0859 sigmay = 43.3f;
0860 } else {
0861 sigmay = sigma;
0862 }
0863
0864
0865
0866 chi21max = fmax(chi21min, (double)templ.chi2yminone());
0867 chi2ymin -= chi21max;
0868 if (chi2ymin < 0.) {
0869 chi2ymin = 0.;
0870 }
0871
0872 meany = fmax(mean1pix, (double)templ.chi2yavgone());
0873 hchi2 = chi2ymin / 2.;
0874 hndof = meany / 2.;
0875 proby = 1. - TMath::Gamma(hndof, hchi2);
0876
0877 } else {
0878
0879
0880 binl = minbin - 1;
0881 binh = binl + 2;
0882 if (binl < fybin) {
0883 binl = fybin;
0884 }
0885 if (binh > lybin) {
0886 binh = lybin;
0887 }
0888 ssa = 0.;
0889 sa2 = 0.;
0890 ssba = 0.;
0891 saba = 0.;
0892 sba2 = 0.;
0893 for (i = fypix - 2; i <= lypix + 2; ++i) {
0894 ssa += ysw[i] * ytemp[binl][i];
0895 sa2 += ytemp[binl][i] * ytemp[binl][i] * yw2[i];
0896 ssba += ysw[i] * (ytemp[binh][i] - ytemp[binl][i]);
0897 saba += ytemp[binl][i] * (ytemp[binh][i] - ytemp[binl][i]) * yw2[i];
0898 sba2 += (ytemp[binh][i] - ytemp[binl][i]) * (ytemp[binh][i] - ytemp[binl][i]) * yw2[i];
0899 }
0900
0901
0902
0903 rat = (ssba * ssa - ss2 * saba) / (ss2 * sba2 - ssba * ssba);
0904 if (rat < 0.f) {
0905 rat = 0.f;
0906 }
0907 if (rat > 1.f) {
0908 rat = 1.0f;
0909 }
0910 rnorm = (ssa + rat * ssba) / ss2;
0911
0912
0913
0914 qfy = ysum[fypix];
0915 if (yd[fypix]) {
0916 qfy += ysum[fypix + 1];
0917 }
0918 if (logypx > 1) {
0919 qly = ysum[lypix];
0920 if (yd[lypix - 1]) {
0921 qly += ysum[lypix - 1];
0922 }
0923 } else {
0924 qly = qfy;
0925 }
0926
0927
0928
0929 float qyfrac = (qfy - qly) / (qfy + qly);
0930 if (goodEdgeAlgo) {
0931 bias = templ.yflcorr(binq, qyfrac) + templ.ygx0(binq);
0932 } else {
0933 bias = templ.yflcorr(binq, qyfrac) + templ.yavg(binq);
0934 }
0935
0936
0937 yrec = (0.125f * binl + BHY - 2.5f + rat * (binh - binl) * 0.125f - (float)shifty + originy) * ysize - bias;
0938 if (goodEdgeAlgo) {
0939 sigmay = templ.ygsig(binq);
0940 } else {
0941 sigmay = templ.yrms(binq);
0942 }
0943
0944
0945
0946 if (rnorm <= 0.) {
0947 LOGERROR("SiPixelTemplateReco") << "illegal chi2y normalization (2) = " << rnorm << ENDL;
0948 rnorm = 1.;
0949 }
0950 chi2y = ss2 - 2. / rnorm * ssa - 2. / rnorm * rat * ssba +
0951 (sa2 + 2. * rat * saba + rat * rat * sba2) / (rnorm * rnorm) - templ.chi2ymin(binq);
0952 if (chi2y < 0.0) {
0953 chi2y = 0.0;
0954 }
0955 meany = templ.chi2yavg(binq);
0956 if (meany < 0.01) {
0957 meany = 0.01;
0958 }
0959
0960
0961
0962 hchi2 = chi2y / 2.;
0963 hndof = meany / 2.;
0964 proby = 1. - TMath::Gamma(hndof, hchi2);
0965 }
0966
0967
0968
0969
0970
0971
0972
0973 if (nxpix > logxpx) {
0974 i = fxpix;
0975 while (i < lxpix) {
0976 if (xd[i] && !xd[i + 1]) {
0977 for (j = fxbin; j <= lxbin; ++j) {
0978
0979
0980 sigavg = (xtemp[j][i] + xtemp[j][i + 1]) / 2.f;
0981 xtemp[j][i] = sigavg;
0982 xtemp[j][i + 1] = sigavg;
0983 }
0984 i += 2;
0985 } else {
0986 ++i;
0987 }
0988 }
0989 }
0990
0991
0992
0993 sxthr = 1.1f * templ.sxmax();
0994
0995
0996 std::sort(&xsort[0], &xsort[logxpx]);
0997 if (logxpx == 1) {
0998 sxthr = 1.01f * xsort[0];
0999 } else {
1000 if (xsort[1] > sxthr) {
1001 sxthr = 1.01f * xsort[1];
1002 }
1003 }
1004
1005
1006
1007
1008 templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
1009
1010
1011
1012 chi2xmin = 1.e15;
1013 for (i = fxbin; i <= lxbin; ++i) {
1014 chi2xbin[i] = -1.e15f;
1015 }
1016 ss2 = 0.f;
1017 for (i = fxpix - 2; i <= lxpix + 2; ++i) {
1018 xw2[i] = 1.f / xsig2[i];
1019 xsw[i] = xsum[i] * xw2[i];
1020 ss2 += xsum[i] * xsw[i];
1021 }
1022 minbin = -1;
1023 deltaj = djx;
1024 jmin = fxbin;
1025 jmax = lxbin;
1026 while (deltaj > 0) {
1027 for (j = jmin; j <= jmax; j += deltaj) {
1028 if (chi2xbin[j] < -100.f) {
1029 ssa = 0.f;
1030 sa2 = 0.f;
1031
1032 for (i = fxpix - 2; i <= lxpix + 2; ++i) {
1033
1034 ssa += xsw[i] * xtemp[j][i];
1035 sa2 += xtemp[j][i] * xtemp[j][i] * xw2[i];
1036 }
1037
1038 rat = ssa / ss2;
1039 if (rat <= 0.f) {
1040 LOGERROR("SiPixelTemplateReco") << "illegal chi2xmin normalization (1) = " << rat << ENDL;
1041 rat = 1.;
1042 }
1043 chi2xbin[j] = ss2 - 2.f * ssa / rat + sa2 / (rat * rat);
1044 }
1045 if (chi2xbin[j] < chi2xmin) {
1046 chi2xmin = chi2xbin[j];
1047 minbin = j;
1048 }
1049 }
1050 deltaj /= 2;
1051 if (minbin > fxbin) {
1052 jmin = minbin - deltaj;
1053 } else {
1054 jmin = fxbin;
1055 }
1056 if (minbin < lxbin) {
1057 jmax = minbin + deltaj;
1058 } else {
1059 jmax = lxbin;
1060 }
1061 }
1062
1063 if (theVerboseLevel > 9) {
1064 LOGDEBUG("SiPixelTemplateReco") << "minbin " << minbin << " chi2xmin = " << chi2xmin << ENDL;
1065 }
1066
1067
1068
1069 if (logxpx == 1) {
1070 if (nxpix == 1) {
1071 delta = templ.dxone();
1072 sigma = templ.sxone();
1073 } else {
1074 delta = templ.dxtwo();
1075 sigma = templ.sxtwo();
1076 }
1077 xrec = 0.5 * (fxpix + lxpix - 2 * shiftx + 2. * originx) * xsize - delta;
1078 if (sigma <= 0.) {
1079 sigmax = 28.9;
1080 } else {
1081 sigmax = sigma;
1082 }
1083
1084
1085
1086 chi21max = fmax(chi21min, (double)templ.chi2xminone());
1087 chi2xmin -= chi21max;
1088 if (chi2xmin < 0.) {
1089 chi2xmin = 0.;
1090 }
1091 meanx = fmax(mean1pix, (double)templ.chi2xavgone());
1092 hchi2 = chi2xmin / 2.;
1093 hndof = meanx / 2.;
1094 probx = 1. - TMath::Gamma(hndof, hchi2);
1095
1096 } else {
1097
1098
1099 binl = minbin - 1;
1100 binh = binl + 2;
1101 if (binl < fxbin) {
1102 binl = fxbin;
1103 }
1104 if (binh > lxbin) {
1105 binh = lxbin;
1106 }
1107 ssa = 0.;
1108 sa2 = 0.;
1109 ssba = 0.;
1110 saba = 0.;
1111 sba2 = 0.;
1112 for (i = fxpix - 2; i <= lxpix + 2; ++i) {
1113 ssa += xsw[i] * xtemp[binl][i];
1114 sa2 += xtemp[binl][i] * xtemp[binl][i] * xw2[i];
1115 ssba += xsw[i] * (xtemp[binh][i] - xtemp[binl][i]);
1116 saba += xtemp[binl][i] * (xtemp[binh][i] - xtemp[binl][i]) * xw2[i];
1117 sba2 += (xtemp[binh][i] - xtemp[binl][i]) * (xtemp[binh][i] - xtemp[binl][i]) * xw2[i];
1118 }
1119
1120
1121
1122 rat = (ssba * ssa - ss2 * saba) / (ss2 * sba2 - ssba * ssba);
1123 if (rat < 0.f) {
1124 rat = 0.f;
1125 }
1126 if (rat > 1.f) {
1127 rat = 1.0f;
1128 }
1129 rnorm = (ssa + rat * ssba) / ss2;
1130
1131
1132
1133 qfx = xsum[fxpix];
1134 if (xd[fxpix]) {
1135 qfx += xsum[fxpix + 1];
1136 }
1137 if (logxpx > 1) {
1138 qlx = xsum[lxpix];
1139 if (xd[lxpix - 1]) {
1140 qlx += xsum[lxpix - 1];
1141 }
1142 } else {
1143 qlx = qfx;
1144 }
1145
1146
1147
1148 float qxfrac = (qfx - qlx) / (qfx + qlx);
1149 bias = templ.xflcorr(binq, qxfrac) + templ.xavg(binq);
1150
1151
1152
1153 xrec = (0.125f * binl + BHX - 2.5f + rat * (binh - binl) * 0.125f - (float)shiftx + originx) * xsize - bias;
1154 sigmax = templ.xrms(binq);
1155
1156
1157
1158 if (rnorm <= 0.) {
1159 LOGERROR("SiPixelTemplateReco") << "illegal chi2x normalization (2) = " << rnorm << ENDL;
1160 rnorm = 1.;
1161 }
1162 chi2x = ss2 - 2. / rnorm * ssa - 2. / rnorm * rat * ssba +
1163 (sa2 + 2. * rat * saba + rat * rat * sba2) / (rnorm * rnorm) - templ.chi2xmin(binq);
1164 if (chi2x < 0.0) {
1165 chi2x = 0.0;
1166 }
1167 meanx = templ.chi2xavg(binq);
1168 if (meanx < 0.01) {
1169 meanx = 0.01;
1170 }
1171
1172
1173
1174 hchi2 = chi2x / 2.;
1175 hndof = meanx / 2.;
1176 probx = 1. - TMath::Gamma(hndof, hchi2);
1177 }
1178
1179
1180
1181 if (proby < probmin) {
1182 proby = probmin;
1183 }
1184 if (probx < probmin) {
1185 probx = probmin;
1186 }
1187
1188
1189
1190 if (calc_probQ) {
1191
1192
1193 templ.vavilov_pars(mpv, sigmaQ, kappa);
1194 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1195 if ((sigmaQ <= 0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
1196 throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::Vavilov parameters mpv/sigmaQ/kappa = " << mpv << "/"
1197 << sigmaQ << "/" << kappa << std::endl;
1198 }
1199 #else
1200 assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
1201 #endif
1202 xvav = ((double)qtotal - mpv) / sigmaQ;
1203 beta2 = 1.;
1204 if (use_VVIObj) {
1205
1206 VVIObjF vvidist(kappa);
1207 prvav = vvidist.fcn(xvav);
1208 } else {
1209
1210 prvav = TMath::VavilovI(xvav, kappa, beta2);
1211 }
1212
1213
1214
1215 probQ = 1. - prvav;
1216 if (probQ < probQmin) {
1217 probQ = probQmin;
1218 }
1219 } else {
1220 probQ = -1;
1221 }
1222
1223 return 0;
1224 }
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261 int SiPixelTemplateReco::PixelTempReco1D(int id,
1262 float cotalpha,
1263 float cotbeta,
1264 float locBz,
1265 float locBx,
1266 ClusMatrix& cluster,
1267 SiPixelTemplate& templ,
1268 float& yrec,
1269 float& sigmay,
1270 float& proby,
1271 float& xrec,
1272 float& sigmax,
1273 float& probx,
1274 int& qbin,
1275 int speed,
1276 float& probQ,
1277 bool goodEdgeAlgo)
1278
1279 {
1280
1281 const bool deadpix = false;
1282 std::vector<std::pair<int, int> > zeropix;
1283 int nypix, nxpix;
1284
1285 return SiPixelTemplateReco::PixelTempReco1D(id,
1286 cotalpha,
1287 cotbeta,
1288 locBz,
1289 locBx,
1290 cluster,
1291 templ,
1292 yrec,
1293 sigmay,
1294 proby,
1295 xrec,
1296 sigmax,
1297 probx,
1298 qbin,
1299 speed,
1300 deadpix,
1301 zeropix,
1302 probQ,
1303 nypix,
1304 nxpix,
1305 goodEdgeAlgo);
1306
1307 }
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343 int SiPixelTemplateReco::PixelTempReco1D(int id,
1344 float cotalpha,
1345 float cotbeta,
1346 float locBz,
1347 float locBx,
1348 ClusMatrix& cluster,
1349 SiPixelTemplate& templ,
1350 float& yrec,
1351 float& sigmay,
1352 float& proby,
1353 float& xrec,
1354 float& sigmax,
1355 float& probx,
1356 int& qbin,
1357 int speed,
1358 float& probQ)
1359
1360 {
1361
1362 const bool deadpix = false;
1363 std::vector<std::pair<int, int> > zeropix;
1364 int nypix, nxpix;
1365 const bool goodEdgeAlgo = false;
1366
1367 return SiPixelTemplateReco::PixelTempReco1D(id,
1368 cotalpha,
1369 cotbeta,
1370 locBz,
1371 locBx,
1372 cluster,
1373 templ,
1374 yrec,
1375 sigmay,
1376 proby,
1377 xrec,
1378 sigmax,
1379 probx,
1380 qbin,
1381 speed,
1382 deadpix,
1383 zeropix,
1384 probQ,
1385 nypix,
1386 nxpix,
1387 goodEdgeAlgo);
1388
1389 }
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419 int SiPixelTemplateReco::PixelTempReco1D(int id,
1420 float cotalpha,
1421 float cotbeta,
1422 ClusMatrix& cluster,
1423 SiPixelTemplate& templ,
1424 float& yrec,
1425 float& sigmay,
1426 float& proby,
1427 float& xrec,
1428 float& sigmax,
1429 float& probx,
1430 int& qbin,
1431 int speed,
1432 float& probQ)
1433
1434 {
1435
1436 const bool deadpix = false;
1437 const bool goodEdgeAlgo = false;
1438 std::vector<std::pair<int, int> > zeropix;
1439 int nypix, nxpix;
1440 float locBx, locBz;
1441 locBx = 1.f;
1442 if (cotbeta < 0.f) {
1443 locBx = -1.f;
1444 }
1445 locBz = locBx;
1446 if (cotalpha < 0.f) {
1447 locBz = -locBx;
1448 }
1449
1450 return SiPixelTemplateReco::PixelTempReco1D(id,
1451 cotalpha,
1452 cotbeta,
1453 locBz,
1454 locBx,
1455 cluster,
1456 templ,
1457 yrec,
1458 sigmay,
1459 proby,
1460 xrec,
1461 sigmax,
1462 probx,
1463 qbin,
1464 speed,
1465 deadpix,
1466 zeropix,
1467 probQ,
1468 nypix,
1469 nxpix,
1470 goodEdgeAlgo);
1471
1472 }
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499 int SiPixelTemplateReco::PixelTempReco1D(int id,
1500 float cotalpha,
1501 float cotbeta,
1502 ClusMatrix& cluster,
1503 SiPixelTemplate& templ,
1504 float& yrec,
1505 float& sigmay,
1506 float& proby,
1507 float& xrec,
1508 float& sigmax,
1509 float& probx,
1510 int& qbin,
1511 int speed)
1512
1513 {
1514
1515 const bool deadpix = false;
1516 const bool goodEdgeAlgo = false;
1517 std::vector<std::pair<int, int> > zeropix;
1518 int nypix, nxpix;
1519 float locBx, locBz;
1520 locBx = 1.f;
1521 if (cotbeta < 0.f) {
1522 locBx = -1.f;
1523 }
1524 locBz = locBx;
1525 if (cotalpha < 0.f) {
1526 locBz = -locBx;
1527 }
1528 float probQ;
1529 if (speed < 0)
1530 speed = 0;
1531 if (speed > 3)
1532 speed = 3;
1533
1534 return SiPixelTemplateReco::PixelTempReco1D(id,
1535 cotalpha,
1536 cotbeta,
1537 locBz,
1538 locBx,
1539 cluster,
1540 templ,
1541 yrec,
1542 sigmay,
1543 proby,
1544 xrec,
1545 sigmax,
1546 probx,
1547 qbin,
1548 speed,
1549 deadpix,
1550 zeropix,
1551 probQ,
1552 nypix,
1553 nxpix,
1554 goodEdgeAlgo);
1555
1556 }