Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-20 23:57:10

0001 //
0002 //  SiPixelTemplateReco.cc (Version 10.01)
0003 //
0004 //  Add goodness-of-fit to algorithm, include single pixel clusters in chi2 calculation
0005 //  Try "decapitation" of large single pixels
0006 //  Add correction for (Q_F-Q_L)/(Q_F+Q_L) bias
0007 //  Add cot(beta) reflection to reduce y-entries and more sophisticated x-interpolation
0008 //  Fix small double pixel bug with decapitation (2.41 5-Mar-2007).
0009 //  Fix pseudopixel bug causing possible memory overwrite (2.42 12-Mar-2007)
0010 //  Adjust template binning to span 3 (or 4) central pixels and implement improved (faster) chi2min search
0011 //  Replace internal containers with static arrays
0012 //  Add external threshold to calls to ysigma2 and xsigma2, use sorted signal heights to guarrantee min clust size = 2
0013 //  Use denser search over larger bin range for clusters with big pixels.
0014 //  Use single calls to template object to load template arrays (had been many)
0015 //  Add speed switch to trade-off speed and robustness
0016 //  Add qmin and re-define qbin to flag low-q clusters
0017 //  Add qscale to match charge scales
0018 //  Return error if no pixels in cluster
0019 //  Replace 4 cout's with LogError's
0020 //  Add LogDebug I/O to report various common errors
0021 //  Incorporate "cluster repair" to handle dead pixels
0022 //  Take truncation size from new pixmax information
0023 //  Change to allow template sizes to be changed at compile time
0024 //  Move interpolation range error to LogDebug
0025 //  Add qbin = 5 and change 1-pixel probability to use new template info
0026 //  Add floor for probabilities (no exact zeros)
0027 //  Replace asserts with exceptions in CMSSW
0028 //  Change calling sequence to handle cot(beta)<0 for FPix cosmics
0029 //
0030 //  V7.00 - Decouple BPix and FPix information into separate templates
0031 //  Pass all containers by alias to prevent excessive cpu-usage (V7.01)
0032 //  Slightly modify search bin range to avoid problem with single pixel clusters + large Lorentz drift (V7.02)
0033 //
0034 //  V8.00 - Add 2D probabilities, take pixel sizes from the template
0035 //  V8.05 - Shift 2-D cluster to center on the buffer
0036 //  V8.06 - Add locBz to the 2-D template (causes failover to the simple template when the cotbeta-locBz correlation is incorrect ... ie for non-IP tracks)
0037 //        - include minimum value for prob2D (1.e-30)
0038 //  V8.07 - Tune 2-d probability: consider only pixels above threshold and use threshold value for zero signal pixels (non-zero template)
0039 //  V8.10 - Remove 2-d probability for ineffectiveness and replace with simple cluster charge probability
0040 //  V8.11 - Change probQ to upper tail probability always (rather than two-sided tail probability)
0041 //  V8.20 - Use template cytemp/cxtemp methods to center the data cluster in the right place when the template becomes asymmetric after irradiation
0042 //  V8.25 - Incorporate VIs speed improvements
0043 //  V8.26 - Fix centering problem for small signals
0044 //  V9.00 - Set QProb = Q/Q_avg when calcultion is turned off, use fbin definitions of Qbin
0045 //  V10.00 - Use new template object to reco Phase 1 FPix hits
0046 //  V10.01 - Fix memory overwriting bug
0047 //  V10.10 - Change VVIObjF so it only reads kappa
0048 //
0049 //
0050 //  Created by Morris Swartz on 10/27/06.
0051 //  VVIObjF object simplification by Tamas Vami
0052 //
0053 
0054 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0055 //#include <cmath.h>
0056 #else
0057 #include <math.h>
0058 #endif
0059 #include <algorithm>
0060 #include <vector>
0061 #include <utility>
0062 #include <iostream>
0063 // ROOT::Math has a c++ function that does the probability calc, but only in v5.12 and later
0064 #include "TMath.h"
0065 #include "Math/DistFunc.h"
0066 // Use current version of gsl instead of ROOT::Math
0067 //#include <gsl/gsl_cdf.h>
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 //! Reconstruct the best estimate of the hit position for pixel clusters.
0090 //! \param         id - (input) identifier of the template to use
0091 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
0092 //! \param locBz - (input) the sign of this quantity is used to determine whether to flip cot(beta)<0 quantities from cot(beta)>0 (FPix only)
0093 //!                    for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0
0094 //!                    for Phase 1 FPix IP-related tracks, see next comment
0095 //! \param 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)
0096 //!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0
0097 //!                    for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0//! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
0098 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
0099 //!           origin of local coords (0,0) at center of pixel cluster[0][0].  Set dead pixels to small non-zero values (0.1 e).
0100 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
0101 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
0102 //! \param      templ - (input) the template used in the reconstruction
0103 //! \param       yrec - (output) best estimate of y-coordinate of hit in microns
0104 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
0105 //! \param      proby - (output) probability describing goodness-of-fit for y-reco
0106 //! \param       xrec - (output) best estimate of x-coordinate of hit in microns
0107 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
0108 //! \param      probx - (output) probability describing goodness-of-fit for x-reco
0109 //! \param       qbin - (output) index (0-4) describing the charge of the cluster
0110 //!                       qbin = 0        Q/Q_avg > 1.5   [few % of all hits]
0111 //!                              1  1.5 > Q/Q_avg > 1.0   [~30% of all hits]
0112 //!                              2  1.0 > Q/Q_avg > 0.85  [~30% of all hits]
0113 //!                              3 0.85 > Q/Q_avg > min1  [~30% of all hits]
0114 //!                              4 min1 > Q/Q_avg > min2  [~0.1% of all hits]
0115 //!                              5 min2 > Q/Q_avg         [~0.1% of all hits]
0116 //! \param      speed - (input) switch (-2->5) trading speed vs robustness
0117 //!                     -2       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4),
0118 //!                              calculates Q probability w/ VVIObj (better but slower)
0119 //!                     -1       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4),
0120 //!                              calculates Q probability w/ TMath::VavilovI (poorer but faster)
0121 //!                      0       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4)
0122 //!                      1       faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density
0123 //!                      2       faster yet, searches same range as 1 but at 1/2 density
0124 //!                      3       fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
0125 //!                      4       fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster),
0126 //!                              calculates Q probability w/ VVIObj (better but slower)
0127 //!                      5       fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster),
0128 //!                              calculates Q probability w/ TMath::VavilovI (poorer but faster)
0129 //! \param    deadpix - (input)  bool to indicate that there are dead pixels to be included in the analysis
0130 //! \param    zeropix - (input)  vector of index pairs pointing to the dead pixels
0131 //! \param      probQ - (output) the Vavilov-distribution-based cluster charge probability
0132 //! \param      nypix - (output) the projected y-size of the cluster
0133 //! \param      nxpix - (output) the projected x-size of the cluster
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   // Local variables
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   //int fypix2D, lypix2D, fxpix2D, lxpix2D;
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   // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
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   // First, interpolate the template needed to analyze this cluster
0186   // check to see of the track direction is in the physical range of the loaded template
0187 
0188   if (id >= 0) {  //if id < 0 bypass interpolation (used in calibration)
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   // Check to see if Q probability is selected
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   // Get pixel dimensions from the template (to allow multiple detectors in the future)
0219 
0220   xsize = templ.xsize();
0221   ysize = templ.ysize();
0222 
0223   // Allow Qbin Q/Q_avg fractions to vary to optimize error estimation
0224 
0225   for (i = 0; i < 3; ++i) {
0226     fbin[i] = templ.fbin(i);
0227   }
0228 
0229   // Define size of pseudopixel
0230 
0231   q50 = templ.s50();
0232   pseudopix = 0.2f * q50;
0233 
0234   // Get charge scaling factor
0235 
0236   qscale = templ.qscale();
0237 
0238   // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
0239 
0240   nclusx = cluster.mrow;
0241   nclusy = (int)cluster.mcol;
0242   auto const xdouble = cluster.xdouble;
0243   auto const ydouble = cluster.ydouble;
0244 
0245   // First, rescale all pixel charges and compute total charge
0246   qtotal = 0.f;
0247   for (i = 0; i < nclusx * nclusy; ++i) {
0248     cluster.matrix[i] *= qscale;
0249     qtotal += cluster.matrix[i];
0250   }
0251   // Next, sum the total charge and "decapitate" big pixels
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   // Do the cluster repair here
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       // Do preliminary cluster projection in y
0273       for (j = 0; j < nclusx; ++j) {
0274         ysum[i] += cluster(j, i);
0275       }
0276       if (ysum[i] > 0.f) {
0277         // identify ends of cluster to determine what the missing charge should be
0278         if (i < fypix) {
0279           fypix = i;
0280         }
0281         if (i > lypix) {
0282           lypix = i;
0283         }
0284       }
0285     }
0286 
0287     // Now loop over dead pixel list and "fix" everything
0288 
0289     //First see if the cluster ends are redefined and that we have only one dead pixel per column
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       // count the number of dead pixels in each column
0300       ++nyzero[i];
0301       // allow them to redefine the cluster ends
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     // Now adjust the charge in the dead pixels to sum to 0.5*truncation value in the end columns and the truncation value in the interior columns
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       // Adjust the total cluster charge to reflect the charge of the "repaired" cluster
0339       qtotal += qpixel;
0340     }
0341     // End of cluster repair section
0342   }
0343 
0344   // Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container
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     // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
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   // Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container
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     // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
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   // next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx
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   //    dlengthy = (float)nypix - templ.clsleny();
0427 
0428   // Make sure cluster is continuous
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   // If cluster is longer than max template size, technique fails
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   // Remember these numbers for later
0460 
0461   //fypix2D = fypix;
0462   //lypix2D = lypix;
0463 
0464   // next, center the cluster on template center if necessary
0465 
0466   midpix = (fypix + lypix) / 2;
0467   shifty = templ.cytemp() - midpix;
0468 
0469   // calculate new cluster boundaries
0470 
0471   int lytmp = lypix + shifty;
0472   int fytmp = fypix + shifty;
0473 
0474   // Check the boundaries
0475 
0476   if (fytmp <= 1) {
0477     return 8;
0478   }
0479   if (lytmp >= BYM2) {
0480     return 8;
0481   }
0482 
0483   //  If OK, shift everything
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   // add pseudopixels
0505 
0506   ysum[fypix - 1] = pseudopix;
0507   ysum[fypix - 2] = pseudopix;
0508   ysum[lypix + 1] = pseudopix;
0509   ysum[lypix + 2] = pseudopix;
0510 
0511   // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
0512 
0513   if (ydouble[0]) {
0514     originy = -0.5f;
0515   } else {
0516     originy = 0.f;
0517   }
0518 
0519   // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx
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   //    dlengthx = (float)nxpix - templ.clslenx();
0540 
0541   // Make sure cluster is continuous
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   // If cluster is longer than max template size, technique fails
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   // Remember these numbers for later
0573 
0574   //fxpix2D = fxpix;
0575   //lxpix2D = lxpix;
0576 
0577   // next, center the cluster on template center if necessary
0578 
0579   midpix = (fxpix + lxpix) / 2;
0580   shiftx = templ.cxtemp() - midpix;
0581 
0582   // calculate new cluster boundaries
0583 
0584   int lxtmp = lxpix + shiftx;
0585   int fxtmp = fxpix + shiftx;
0586 
0587   // Check the boundaries
0588 
0589   if (fxtmp <= 1) {
0590     return 9;
0591   }
0592   if (lxtmp >= BXM2) {
0593     return 9;
0594   }
0595 
0596   //  If OK, shift everything
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   // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
0623 
0624   if (xdouble[0]) {
0625     originx = -0.5f;
0626   } else {
0627     originx = 0.f;
0628   }
0629 
0630   // uncertainty and final corrections depend upon total charge bin
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   // Return the charge bin via the parameter list unless the charge is too small (then flag it)
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   // Next, copy the y- and x-templates to local arrays
0663 
0664   // First, decide on chi^2 min search parameters
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   // Now do the copies
0728 
0729   templ.ytemp(fybin, lybin, ytemp);
0730 
0731   templ.xtemp(fxbin, lxbin, xtemp);
0732 
0733   // Do the y-reconstruction first
0734 
0735   // Apply the first-pass template algorithm to all clusters
0736 
0737   // Modify the template if double pixels are present
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           // Sum the adjacent cells and put the average signal in both
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   // Define the maximum signal to allow before de-weighting a pixel
0758 
0759   sythr = 1.1 * (templ.symax());
0760 
0761   // Make sure that there will be at least two pixels that are not de-weighted
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   // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
0773 
0774   //    for(i=0; i<BYSIZE; ++i) { ysig2[i] = 0.;}
0775   templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
0776 
0777   // Find the template bin that minimizes the Chi^2
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   // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
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     // Do probability calculation for one-pixel clusters
0851 
0852     chi21max = fmax(chi21min, (double)templ.chi2yminone());
0853     chi2ymin -= chi21max;
0854     if (chi2ymin < 0.) {
0855       chi2ymin = 0.;
0856     }
0857     //     proby = gsl_cdf_chisq_Q(chi2ymin, mean1pix);
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     // For cluster > 1 pix, make the second, interpolating pass with the templates
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     // rat is the fraction of the "distance" from template a to template b
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     // Calculate the charges in the first and last pixels
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     //  Now calculate the mean bias correction and uncertainties
0914 
0915     float qyfrac = (qfy - qly) / (qfy + qly);
0916     bias = templ.yflcorr(binq, qyfrac) + templ.yavg(binq);
0917 
0918     // uncertainty and final correction depend upon charge bin
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     // Do goodness of fit test in y
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     // gsl function that calculates the chi^2 tail prob for non-integral dof
0939     //     proby = gsl_cdf_chisq_Q(chi2y, meany);
0940     //     proby = ROOT::Math::chisquared_cdf_c(chi2y, meany);
0941     hchi2 = chi2y / 2.;
0942     hndof = meany / 2.;
0943     proby = 1. - TMath::Gamma(hndof, hchi2);
0944   }
0945 
0946   // Do the x-reconstruction next
0947 
0948   // Apply the first-pass template algorithm to all clusters
0949 
0950   // Modify the template if double pixels are present
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           // Sum the adjacent cells and put the average signal in both
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   // Define the maximum signal to allow before de-weighting a pixel
0971 
0972   sxthr = 1.1f * templ.sxmax();
0973 
0974   // Make sure that there will be at least two pixels that are not de-weighted
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   // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
0985 
0986   //    for(i=0; i<BXSIZE; ++i) { xsig2[i] = 0.; }
0987   templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
0988 
0989   // Find the template bin that minimizes the Chi^2
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         //           std::cout << j << " ";
1011         for (i = fxpix - 2; i <= lxpix + 2; ++i) {
1012           //                std::cout << xtemp[j][i] << ", ";
1013           ssa += xsw[i] * xtemp[j][i];
1014           sa2 += xtemp[j][i] * xtemp[j][i] * xw2[i];
1015         }
1016         //           std::cout << std::endl;
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   // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
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     // Do probability calculation for one-pixel clusters
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     // Now make the second, interpolating pass with the templates
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     // rat is the fraction of the "distance" from template a to template b
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     // Calculate the charges in the first and last pixels
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     //  Now calculate the mean bias correction and uncertainties
1126 
1127     float qxfrac = (qfx - qlx) / (qfx + qlx);
1128     bias = templ.xflcorr(binq, qxfrac) + templ.xavg(binq);
1129 
1130     // uncertainty and final correction depend upon charge bin
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     // Do goodness of fit test in x
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     // gsl function that calculates the chi^2 tail prob for non-integral dof
1151     //     probx = gsl_cdf_chisq_Q(chi2x, meanx);
1152     //     probx = ROOT::Math::chisquared_cdf_c(chi2x, meanx, trx0);
1153     hchi2 = chi2x / 2.;
1154     hndof = meanx / 2.;
1155     probx = 1. - TMath::Gamma(hndof, hchi2);
1156   }
1157 
1158   //  Don't return exact zeros for the probability
1159 
1160   if (proby < probmin) {
1161     proby = probmin;
1162   }
1163   if (probx < probmin) {
1164     probx = probmin;
1165   }
1166 
1167   //  Decide whether to generate a cluster charge probability
1168 
1169   if (calc_probQ) {
1170     // Calculate the Vavilov probability that the cluster charge is OK
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       //  VVIObj is a private port of CERNLIB VVIDIS
1185       VVIObjF vvidist(kappa);
1186       prvav = vvidist.fcn(xvav);
1187     } else {
1188       //  Use faster but less accurate TMath Vavilov distribution function
1189       prvav = TMath::VavilovI(xvav, kappa, beta2);
1190     }
1191     //  Change to upper tail probability
1192     //      if(prvav > 0.5) prvav = 1. - prvav;
1193     //      probQ = (float)(2.*prvav);
1194     probQ = 1. - prvav;
1195     if (probQ < probQmin) {
1196       probQ = probQmin;
1197     }
1198   } else {
1199     probQ = -1;
1200   }
1201 
1202   return 0;
1203 }  // PixelTempReco1D
1204 
1205 // *************************************************************************************************************************************
1206 //  Overload parameter list for compatibility with older versions
1207 //! Reconstruct the best estimate of the hit position for pixel clusters.
1208 //! \param         id - (input) identifier of the template to use
1209 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1210 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1211 //! \param locBz - (input) the sign of this quantity is used to determine whether to flip cot(beta)<0 quantities from cot(beta)>0 (FPix only)
1212 //!                    for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0
1213 //!                    for Phase 1 FPix IP-related tracks, see next comment
1214 //! \param 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)
1215 //!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0
1216 //!                    for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0//! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1217 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
1218 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
1219 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
1220 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
1221 //! \param      templ - (input) the template used in the reconstruction
1222 //! \param       yrec - (output) best estimate of y-coordinate of hit in microns
1223 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
1224 //! \param      proby - (output) probability describing goodness-of-fit for y-reco
1225 //! \param       xrec - (output) best estimate of x-coordinate of hit in microns
1226 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
1227 //! \param      probx - (output) probability describing goodness-of-fit for x-reco
1228 //! \param       qbin - (output) index (0-4) describing the charge of the cluster
1229 //!                     [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin]
1230 //! \param      speed - (input) switch (-1-4) trading speed vs robustness
1231 //!                     -1       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability
1232 //!                      0       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4)
1233 //!                      1       faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density
1234 //!                      2       faster yet, searches same range as 1 but at 1/2 density
1235 //!                      3       fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
1236 //!                      4       fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster), calculates Q probability
1237 //! \param      probQ - (output) the Vavilov-distribution-based cluster charge probability
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   // Local variables
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 }  // PixelTempReco1D
1284 
1285 // *************************************************************************************************************************************
1286 //  Overload parameter list for compatibility with older versions
1287 //! Reconstruct the best estimate of the hit position for pixel clusters.
1288 //! \param         id - (input) identifier of the template to use
1289 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1290 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1291 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
1292 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
1293 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
1294 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
1295 //! \param      templ - (input) the template used in the reconstruction
1296 //! \param       yrec - (output) best estimate of y-coordinate of hit in microns
1297 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
1298 //! \param      proby - (output) probability describing goodness-of-fit for y-reco
1299 //! \param       xrec - (output) best estimate of x-coordinate of hit in microns
1300 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
1301 //! \param      probx - (output) probability describing goodness-of-fit for x-reco
1302 //! \param       qbin - (output) index (0-4) describing the charge of the cluster
1303 //!                     [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin]
1304 //! \param      speed - (input) switch (-1-4) trading speed vs robustness
1305 //!                     -1       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability
1306 //!                      0       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4)
1307 //!                      1       faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density
1308 //!                      2       faster yet, searches same range as 1 but at 1/2 density
1309 //!                      3       fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
1310 //!                      4       fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster), calculates Q probability
1311 //! \param      probQ - (output) the Vavilov-distribution-based cluster charge probability
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   // Local variables
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 }  // PixelTempReco1D
1365 
1366 // *************************************************************************************************************************************
1367 //  Overload parameter list for compatibility with older versions
1368 //! Reconstruct the best estimate of the hit position for pixel clusters.
1369 //! \param         id - (input) identifier of the template to use
1370 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1371 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1372 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
1373 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
1374 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
1375 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
1376 //! \param      templ - (input) the template used in the reconstruction
1377 //! \param       yrec - (output) best estimate of y-coordinate of hit in microns
1378 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
1379 //! \param      proby - (output) probability describing goodness-of-fit for y-reco
1380 //! \param       xrec - (output) best estimate of x-coordinate of hit in microns
1381 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
1382 //! \param      probx - (output) probability describing goodness-of-fit for x-reco
1383 //! \param       qbin - (output) index (0-4) describing the charge of the cluster
1384 //!                     [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin]
1385 //! \param      speed - (input) switch (0-3) trading speed vs robustness
1386 //!                      0       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4)
1387 //!                      1       faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density
1388 //!                      2       faster yet, searches same range as 1 but at 1/2 density
1389 //!                      3       fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
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   // Local variables
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 }  // PixelTempReco1D