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