Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-06-25 05:16:14

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 //  V10.30 - Use "good" end to center the y-projections [for high eta radiation damage],
0049 //         - use Gaussian fit paramters for y-bias and y-errors, remove chi2min_y warning
0050 //
0051 //
0052 //  Created by Morris Swartz on 10/27/06.
0053 //  VVIObjF object simplification by Tamas Vami
0054 //
0055 
0056 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0057 //#include <cmath.h>
0058 #else
0059 #include <math.h>
0060 #endif
0061 #include <algorithm>
0062 #include <vector>
0063 #include <utility>
0064 #include <iostream>
0065 // ROOT::Math has a c++ function that does the probability calc, but only in v5.12 and later
0066 #include "TMath.h"
0067 #include "Math/DistFunc.h"
0068 // Use current version of gsl instead of ROOT::Math
0069 //#include <gsl/gsl_cdf.h>
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 //! Reconstruct the best estimate of the hit position for pixel clusters.
0092 //! \param         id - (input) identifier of the template to use
0093 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
0094 //! \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)
0095 //!                    for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0
0096 //!                    for Phase 1 FPix IP-related tracks, see next comment
0097 //! \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)
0098 //!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0
0099 //!                    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)
0100 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
0101 //!           origin of local coords (0,0) at center of pixel cluster[0][0].  Set dead pixels to small non-zero values (0.1 e).
0102 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
0103 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
0104 //! \param      templ - (input) the template used in the reconstruction
0105 //! \param       yrec - (output) best estimate of y-coordinate of hit in microns
0106 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
0107 //! \param      proby - (output) probability describing goodness-of-fit for y-reco
0108 //! \param       xrec - (output) best estimate of x-coordinate of hit in microns
0109 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
0110 //! \param      probx - (output) probability describing goodness-of-fit for x-reco
0111 //! \param       qbin - (output) index (0-4) describing the charge of the cluster
0112 //!                       qbin = 0        Q/Q_avg > 1.5   [few % of all hits]
0113 //!                              1  1.5 > Q/Q_avg > 1.0   [~30% of all hits]
0114 //!                              2  1.0 > Q/Q_avg > 0.85  [~30% of all hits]
0115 //!                              3 0.85 > Q/Q_avg > min1  [~30% of all hits]
0116 //!                              4 min1 > Q/Q_avg > min2  [~0.1% of all hits]
0117 //!                              5 min2 > Q/Q_avg         [~0.1% of all hits]
0118 //! \param      speed - (input) switch (-2->5) trading speed vs robustness
0119 //!                     -2       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4),
0120 //!                              calculates Q probability w/ VVIObj (better but slower)
0121 //!                     -1       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4),
0122 //!                              calculates Q probability w/ TMath::VavilovI (poorer but faster)
0123 //!                      0       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4)
0124 //!                      1       faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density
0125 //!                      2       faster yet, searches same range as 1 but at 1/2 density
0126 //!                      3       fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
0127 //!                      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),
0128 //!                              calculates Q probability w/ VVIObj (better but slower)
0129 //!                      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),
0130 //!                              calculates Q probability w/ TMath::VavilovI (poorer but faster)
0131 //! \param    deadpix - (input)  bool to indicate that there are dead pixels to be included in the analysis
0132 //! \param    zeropix - (input)  vector of index pairs pointing to the dead pixels
0133 //! \param      probQ - (output) the Vavilov-distribution-based cluster charge probability
0134 //! \param      nypix - (output) the projected y-size of the cluster
0135 //! \param      nxpix - (output) the projected x-size of the cluster
0136 //! \param      goodEdgeAlgo - (input) bool to indicate whether to use centering based on the good end of the cluster (compatible only with templates derived using the same)
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   // Local variables
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   //int fypix2D, lypix2D, fxpix2D, lxpix2D;
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   // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
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   // First, interpolate the template needed to analyze this cluster
0190   // check to see of the track direction is in the physical range of the loaded template
0191 
0192   if (id >= 0) {  //if id < 0 bypass interpolation (used in calibration)
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   // Check to see if Q probability is selected
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   // Get pixel dimensions from the template (to allow multiple detectors in the future)
0223 
0224   xsize = templ.xsize();
0225   ysize = templ.ysize();
0226 
0227   // Allow Qbin Q/Q_avg fractions to vary to optimize error estimation
0228 
0229   for (i = 0; i < 3; ++i) {
0230     fbin[i] = templ.fbin(i);
0231   }
0232 
0233   // Define size of pseudopixel
0234 
0235   q50 = templ.s50();
0236   pseudopix = 0.2f * q50;
0237 
0238   // Get charge scaling factor
0239 
0240   qscale = templ.qscale();
0241 
0242   // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
0243 
0244   nclusx = cluster.mrow;
0245   nclusy = (int)cluster.mcol;
0246   auto const xdouble = cluster.xdouble;
0247   auto const ydouble = cluster.ydouble;
0248 
0249   // First, rescale all pixel charges and compute total charge
0250   qtotal = 0.f;
0251   for (i = 0; i < nclusx * nclusy; ++i) {
0252     cluster.matrix[i] *= qscale;
0253     qtotal += cluster.matrix[i];
0254   }
0255   // Next, sum the total charge and "decapitate" big pixels
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   // Do the cluster repair here
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       // Do preliminary cluster projection in y
0277       for (j = 0; j < nclusx; ++j) {
0278         ysum[i] += cluster(j, i);
0279       }
0280       if (ysum[i] > 0.f) {
0281         // identify ends of cluster to determine what the missing charge should be
0282         if (i < fypix) {
0283           fypix = i;
0284         }
0285         if (i > lypix) {
0286           lypix = i;
0287         }
0288       }
0289     }
0290 
0291     // Now loop over dead pixel list and "fix" everything
0292 
0293     //First see if the cluster ends are redefined and that we have only one dead pixel per column
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       // count the number of dead pixels in each column
0304       ++nyzero[i];
0305       // allow them to redefine the cluster ends
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     // 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
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       // Adjust the total cluster charge to reflect the charge of the "repaired" cluster
0343       qtotal += qpixel;
0344     }
0345     // End of cluster repair section
0346   }
0347 
0348   // Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container
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     // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
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   // Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container
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     // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
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   // next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx
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   //    dlengthy = (float)nypix - templ.clsleny();
0431 
0432   // Make sure cluster is continuous
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   // If cluster is longer than max template size, technique fails
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   // Remember these numbers for later
0464 
0465   //fypix2D = fypix;
0466   //lypix2D = lypix;
0467 
0468   // next, center the cluster on template center if necessary
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   }  //if(goodEdgeAlgo)
0478   else {
0479     midpix = (fypix + lypix) / 2;
0480     shifty = templ.cytemp() - midpix;
0481   }  //else
0482 
0483   // calculate new cluster boundaries
0484 
0485   int lytmp = lypix + shifty;
0486   int fytmp = fypix + shifty;
0487 
0488   // Check the boundaries
0489 
0490   if (fytmp <= 1) {
0491     return 8;
0492   }
0493   if (lytmp >= BYM2) {
0494     return 8;
0495   }
0496 
0497   //  If OK, shift everything
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   // add pseudopixels
0519 
0520   ysum[fypix - 1] = pseudopix;
0521   ysum[fypix - 2] = pseudopix;
0522   ysum[lypix + 1] = pseudopix;
0523   ysum[lypix + 2] = pseudopix;
0524 
0525   // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
0526 
0527   if (ydouble[0]) {
0528     originy = -0.5f;
0529   } else {
0530     originy = 0.f;
0531   }
0532 
0533   // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx
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   //    dlengthx = (float)nxpix - templ.clslenx();
0554 
0555   // Make sure cluster is continuous
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   // If cluster is longer than max template size, technique fails
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   // Remember these numbers for later
0587 
0588   //fxpix2D = fxpix;
0589   //lxpix2D = lxpix;
0590 
0591   // next, center the cluster on template center if necessary
0592 
0593   midpix = (fxpix + lxpix) / 2;
0594   shiftx = templ.cxtemp() - midpix;
0595 
0596   // calculate new cluster boundaries
0597 
0598   int lxtmp = lxpix + shiftx;
0599   int fxtmp = fxpix + shiftx;
0600 
0601   // Check the boundaries
0602 
0603   if (fxtmp <= 1) {
0604     return 9;
0605   }
0606   if (lxtmp >= BXM2) {
0607     return 9;
0608   }
0609 
0610   //  If OK, shift everything
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   // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
0637 
0638   if (xdouble[0]) {
0639     originx = -0.5f;
0640   } else {
0641     originx = 0.f;
0642   }
0643 
0644   // uncertainty and final corrections depend upon total charge bin
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   // Return the charge bin via the parameter list unless the charge is too small (then flag it)
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   // Next, copy the y- and x-templates to local arrays
0677 
0678   // First, decide on chi^2 min search parameters
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   // Now do the copies
0742 
0743   templ.ytemp(fybin, lybin, ytemp);
0744 
0745   templ.xtemp(fxbin, lxbin, xtemp);
0746 
0747   // Do the y-reconstruction first
0748 
0749   // Apply the first-pass template algorithm to all clusters
0750 
0751   // Modify the template if double pixels are present
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           // Sum the adjacent cells and put the average signal in both
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   // Define the maximum signal to allow before de-weighting a pixel
0772 
0773   sythr = 1.1 * (templ.symax());
0774 
0775   // Make sure that there will be at least two pixels that are not de-weighted
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   // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
0787 
0788   //    for(i=0; i<BYSIZE; ++i) { ysig2[i] = 0.;}
0789   templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
0790 
0791   // Find the template bin that minimizes the Chi^2
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   // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
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     // Do probability calculation for one-pixel clusters
0865 
0866     chi21max = fmax(chi21min, (double)templ.chi2yminone());
0867     chi2ymin -= chi21max;
0868     if (chi2ymin < 0.) {
0869       chi2ymin = 0.;
0870     }
0871     //     proby = gsl_cdf_chisq_Q(chi2ymin, mean1pix);
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     // For cluster > 1 pix, make the second, interpolating pass with the templates
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     // rat is the fraction of the "distance" from template a to template b
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     // Calculate the charges in the first and last pixels
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     //  Now calculate the mean bias correction and uncertainties
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     // uncertainty and final correction depend upon charge bin
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     // Do goodness of fit test in y
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     // gsl function that calculates the chi^2 tail prob for non-integral dof
0960     //     proby = gsl_cdf_chisq_Q(chi2y, meany);
0961     //     proby = ROOT::Math::chisquared_cdf_c(chi2y, meany);
0962     hchi2 = chi2y / 2.;
0963     hndof = meany / 2.;
0964     proby = 1. - TMath::Gamma(hndof, hchi2);
0965   }
0966 
0967   // Do the x-reconstruction next
0968 
0969   // Apply the first-pass template algorithm to all clusters
0970 
0971   // Modify the template if double pixels are present
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           // Sum the adjacent cells and put the average signal in both
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   // Define the maximum signal to allow before de-weighting a pixel
0992 
0993   sxthr = 1.1f * templ.sxmax();
0994 
0995   // Make sure that there will be at least two pixels that are not de-weighted
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   // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1006 
1007   //    for(i=0; i<BXSIZE; ++i) { xsig2[i] = 0.; }
1008   templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
1009 
1010   // Find the template bin that minimizes the Chi^2
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         //           std::cout << j << " ";
1032         for (i = fxpix - 2; i <= lxpix + 2; ++i) {
1033           //                std::cout << xtemp[j][i] << ", ";
1034           ssa += xsw[i] * xtemp[j][i];
1035           sa2 += xtemp[j][i] * xtemp[j][i] * xw2[i];
1036         }
1037         //           std::cout << std::endl;
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   // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
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     // Do probability calculation for one-pixel clusters
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     // Now make the second, interpolating pass with the templates
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     // rat is the fraction of the "distance" from template a to template b
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     // Calculate the charges in the first and last pixels
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     //  Now calculate the mean bias correction and uncertainties
1147 
1148     float qxfrac = (qfx - qlx) / (qfx + qlx);
1149     bias = templ.xflcorr(binq, qxfrac) + templ.xavg(binq);
1150 
1151     // uncertainty and final correction depend upon charge bin
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     // Do goodness of fit test in x
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     // gsl function that calculates the chi^2 tail prob for non-integral dof
1172     //     probx = gsl_cdf_chisq_Q(chi2x, meanx);
1173     //     probx = ROOT::Math::chisquared_cdf_c(chi2x, meanx, trx0);
1174     hchi2 = chi2x / 2.;
1175     hndof = meanx / 2.;
1176     probx = 1. - TMath::Gamma(hndof, hchi2);
1177   }
1178 
1179   //  Don't return exact zeros for the probability
1180 
1181   if (proby < probmin) {
1182     proby = probmin;
1183   }
1184   if (probx < probmin) {
1185     probx = probmin;
1186   }
1187 
1188   //  Decide whether to generate a cluster charge probability
1189 
1190   if (calc_probQ) {
1191     // Calculate the Vavilov probability that the cluster charge is OK
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       //  VVIObj is a private port of CERNLIB VVIDIS
1206       VVIObjF vvidist(kappa);
1207       prvav = vvidist.fcn(xvav);
1208     } else {
1209       //  Use faster but less accurate TMath Vavilov distribution function
1210       prvav = TMath::VavilovI(xvav, kappa, beta2);
1211     }
1212     //  Change to upper tail probability
1213     //      if(prvav > 0.5) prvav = 1. - prvav;
1214     //      probQ = (float)(2.*prvav);
1215     probQ = 1. - prvav;
1216     if (probQ < probQmin) {
1217       probQ = probQmin;
1218     }
1219   } else {
1220     probQ = -1;
1221   }
1222 
1223   return 0;
1224 }  // PixelTempReco1D
1225 
1226 // *************************************************************************************************************************************
1227 //  Overload parameter list for compatibility with older versions
1228 //! Reconstruct the best estimate of the hit position for pixel clusters.
1229 //! \param         id - (input) identifier of the template to use
1230 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1231 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1232 //! \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)
1233 //!                    for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0
1234 //!                    for Phase 1 FPix IP-related tracks, see next comment
1235 //! \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)
1236 //!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0
1237 //!                    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)
1238 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
1239 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
1240 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
1241 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
1242 //! \param      templ - (input) the template used in the reconstruction
1243 //! \param       yrec - (output) best estimate of y-coordinate of hit in microns
1244 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
1245 //! \param      proby - (output) probability describing goodness-of-fit for y-reco
1246 //! \param       xrec - (output) best estimate of x-coordinate of hit in microns
1247 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
1248 //! \param      probx - (output) probability describing goodness-of-fit for x-reco
1249 //! \param       qbin - (output) index (0-4) describing the charge of the cluster
1250 //!                     [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]
1251 //! \param      speed - (input) switch (-1-4) trading speed vs robustness
1252 //!                     -1       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability
1253 //!                      0       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4)
1254 //!                      1       faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density
1255 //!                      2       faster yet, searches same range as 1 but at 1/2 density
1256 //!                      3       fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
1257 //!                      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
1258 //! \param      probQ - (output) the Vavilov-distribution-based cluster charge probability
1259 //! \param      goodEdgeAlgo - (input) bool to indicate whether to use centering based on the good end of the cluster (compatible only with templates derived using the same)
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   // Local variables
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 }  // PixelTempReco1D
1308 
1309 // *************************************************************************************************************************************
1310 //  Overload parameter list for compatibility with older versions
1311 //! Reconstruct the best estimate of the hit position for pixel clusters.
1312 //! \param         id - (input) identifier of the template to use
1313 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1314 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1315 //! \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)
1316 //!                    for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0
1317 //!                    for Phase 1 FPix IP-related tracks, see next comment
1318 //! \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)
1319 //!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0
1320 //!                    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)
1321 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
1322 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
1323 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
1324 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
1325 //! \param      templ - (input) the template used in the reconstruction
1326 //! \param       yrec - (output) best estimate of y-coordinate of hit in microns
1327 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
1328 //! \param      proby - (output) probability describing goodness-of-fit for y-reco
1329 //! \param       xrec - (output) best estimate of x-coordinate of hit in microns
1330 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
1331 //! \param      probx - (output) probability describing goodness-of-fit for x-reco
1332 //! \param       qbin - (output) index (0-4) describing the charge of the cluster
1333 //!                     [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]
1334 //! \param      speed - (input) switch (-1-4) trading speed vs robustness
1335 //!                     -1       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability
1336 //!                      0       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4)
1337 //!                      1       faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density
1338 //!                      2       faster yet, searches same range as 1 but at 1/2 density
1339 //!                      3       fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
1340 //!                      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
1341 //! \param      probQ - (output) the Vavilov-distribution-based cluster charge probability
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   // Local variables
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 }  // PixelTempReco1D
1390 
1391 // *************************************************************************************************************************************
1392 //  Overload parameter list for compatibility with older versions
1393 //! Reconstruct the best estimate of the hit position for pixel clusters.
1394 //! \param         id - (input) identifier of the template to use
1395 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1396 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1397 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
1398 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
1399 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
1400 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
1401 //! \param      templ - (input) the template used in the reconstruction
1402 //! \param       yrec - (output) best estimate of y-coordinate of hit in microns
1403 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
1404 //! \param      proby - (output) probability describing goodness-of-fit for y-reco
1405 //! \param       xrec - (output) best estimate of x-coordinate of hit in microns
1406 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
1407 //! \param      probx - (output) probability describing goodness-of-fit for x-reco
1408 //! \param       qbin - (output) index (0-4) describing the charge of the cluster
1409 //!                     [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]
1410 //! \param      speed - (input) switch (-1-4) trading speed vs robustness
1411 //!                     -1       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability
1412 //!                      0       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4)
1413 //!                      1       faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density
1414 //!                      2       faster yet, searches same range as 1 but at 1/2 density
1415 //!                      3       fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
1416 //!                      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
1417 //! \param      probQ - (output) the Vavilov-distribution-based cluster charge probability
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   // Local variables
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 }  // PixelTempReco1D
1473 
1474 // *************************************************************************************************************************************
1475 //  Overload parameter list for compatibility with older versions
1476 //! Reconstruct the best estimate of the hit position for pixel clusters.
1477 //! \param         id - (input) identifier of the template to use
1478 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1479 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1480 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
1481 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
1482 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
1483 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
1484 //! \param      templ - (input) the template used in the reconstruction
1485 //! \param       yrec - (output) best estimate of y-coordinate of hit in microns
1486 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
1487 //! \param      proby - (output) probability describing goodness-of-fit for y-reco
1488 //! \param       xrec - (output) best estimate of x-coordinate of hit in microns
1489 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
1490 //! \param      probx - (output) probability describing goodness-of-fit for x-reco
1491 //! \param       qbin - (output) index (0-4) describing the charge of the cluster
1492 //!                     [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]
1493 //! \param      speed - (input) switch (0-3) trading speed vs robustness
1494 //!                      0       totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4)
1495 //!                      1       faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density
1496 //!                      2       faster yet, searches same range as 1 but at 1/2 density
1497 //!                      3       fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
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   // Local variables
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 }  // PixelTempReco1D