Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 //
0002 //  SiPixelTemplateSplit.cc (Version 2.30)
0003 //
0004 //  Procedure to fit two templates (same angle hypotheses) to a single cluster
0005 //  Return two x- and two y-coordinates for the cluster
0006 //
0007 //  Created by Morris Swartz on 04/10/08.
0008 //
0009 //  Incorporate "cluster repair" to handle dead pixels
0010 //  Take truncation size from new pixmax information
0011 //  Change to allow template sizes to be changed at compile time
0012 //  Move interpolation range error to LogDebug
0013 //  Add q2bin = 5 and change 1-pixel probability to use new template info
0014 //  Add floor for probabilities (no exact zeros)
0015 //  Add ambiguity resolution with crude 2-D templates (v2.00)
0016 //  Pass all containers by alias to prevent excessive cpu-usage (v2.01)
0017 //  Add ambiguity resolution with fancy 2-D templates (v2.10)
0018 //  Make small change to indices for ambiguity resolution (v2.11)
0019 //  Tune x and y errors separately (v2.12)
0020 //  Use template cytemp/cxtemp methods to center the data cluster in the right place when the templates become asymmetric after irradiation (v2.20)
0021 //  Add charge probability to the splitter [tests consistency with a two-hit merged cluster hypothesis]  (v2.20)
0022 //  Improve likelihood normalization slightly (v2.21)
0023 //  Replace hardwired pixel size derived errors with ones from templated pixel sizes (v2.22)
0024 //  Add shape and charge probabilities for the merged cluster hypothesis (v2.23)
0025 //  Incorporate VI-like speed improvements (v2.25)
0026 //  Improve speed by eliminating the triple index boost::multiarray objects and add speed switch to optimize the algorithm (v2.30)
0027 //  Change VVIObjF so it only reads kappa
0028 //
0029 
0030 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0031 //#include <cmath.h>
0032 #else
0033 #include <math.h>
0034 #endif
0035 #include <algorithm>
0036 #include <vector>
0037 #include <list>
0038 #include <iostream>
0039 // ROOT::Math has a c++ function that does the probability calc, but only in v5.12 and later
0040 #include "Math/DistFunc.h"
0041 #include "TMath.h"
0042 // Use current version of gsl instead of ROOT::Math
0043 //#include <gsl/gsl_cdf.h>
0044 
0045 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0046 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateSplit.h"
0047 #include "RecoLocalTracker/SiPixelRecHits/interface/VVIObj.h"
0048 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0049 #define LOGERROR(x) edm::LogError(x)
0050 #define LOGDEBUG(x) LogDebug(x)
0051 static constexpr int theVerboseLevel = 2;
0052 #define ENDL " "
0053 #else
0054 #include "SiPixelTemplateSplit.h"
0055 #include "VVIObj.h"
0056 //#include "SiPixelTemplate2D.h"
0057 //#include "SimpleTemplate2D.h"
0058 //static int theVerboseLevel = {2};
0059 #define LOGERROR(x) std::cout << x << ": "
0060 #define LOGDEBUG(x) std::cout << x << ": "
0061 #define ENDL std::endl
0062 #endif
0063 
0064 using namespace SiPixelTemplateSplit;
0065 
0066 // *************************************************************************************************************************************
0067 //! Reconstruct the best estimate of the hit positions for pixel clusters.
0068 //! \param         id - (input) identifier of the template to use
0069 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
0070 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
0071 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
0072 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
0073 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
0074 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
0075 //! \param      templ - (input) the template used in the reconstruction
0076 //! \param      yrec1 - (output) best estimate of first y-coordinate of hit in microns
0077 //! \param      yrec2 - (output) best estimate of second y-coordinate of hit in microns
0078 //! \param     sigmay - (output) best estimate of uncertainty on yrec1 and yrec2 in microns
0079 //! \param     prob2y - (output) probability describing goodness-of-fit to a merged cluster hypothesis for y-reco
0080 //! \param      xrec1 - (output) best estimate of first x-coordinate of hit in microns
0081 //! \param      xrec2 - (output) best estimate of second x-coordinate of hit in microns
0082 //! \param     sigmax - (output) best estimate of uncertainty on xrec1 and xrec2 in microns
0083 //! \param     prob2x - (output) probability describing goodness-of-fit to a merged cluster hypothesis for x-reco
0084 //! \param       q2bin - (output) index (0-4) describing the charge of the cluster assuming a merged 2-hit cluster hypothesis
0085 //!                     [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]
0086 //! \param     prob2Q - (output) probability that the cluster charge is compatible with a 2-hit merging
0087 //! \param    resolve - (input) if true, use 2-D chisquare to resolve the 2-fold x-y association ambiguity (costs a factor of two in speed)
0088 //! \param      speed - (input) switch (-1->2) trading speed vs robustness
0089 //!                     -1       totally bombproof, searches the entire ranges of template bins,
0090 //!                              calculates Q probability w/ VVIObj
0091 //!                      0       totally bombproof, searches the entire template bin range at full density (no Qprob)
0092 //!                      1       faster, searches same range as 0 but at 1/2 density
0093 //!                      2       fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
0094 //! \param    dchisq - (output) the delta chisquare estimator used to break the degeneracy (0 means no discrimination, larger than 0.1 is good)
0095 //! \param    deadpix - (input) bool to indicate that there are dead pixels to be included in the analysis
0096 //! \param    zeropix - (input) vector of index pairs pointing to the dead pixels
0097 // *************************************************************************************************************************************
0098 int SiPixelTemplateSplit::PixelTempSplit(int id,
0099                                          float cotalpha,
0100                                          float cotbeta,
0101                                          array_2d& clust,
0102                                          std::vector<bool>& ydouble,
0103                                          std::vector<bool>& xdouble,
0104                                          SiPixelTemplate& templ,
0105                                          float& yrec1,
0106                                          float& yrec2,
0107                                          float& sigmay,
0108                                          float& prob2y,
0109                                          float& xrec1,
0110                                          float& xrec2,
0111                                          float& sigmax,
0112                                          float& prob2x,
0113                                          int& q2bin,
0114                                          float& prob2Q,
0115                                          bool resolve,
0116                                          int speed,
0117                                          float& dchisq,
0118                                          bool deadpix,
0119                                          std::vector<std::pair<int, int> >& zeropix,
0120                                          SiPixelTemplate2D& templ2D)
0121 
0122 {
0123   // Local variables
0124   int i, j, k, binq, midpix, fypix, nypix, lypix, logypx, lparm;
0125   int fxpix, nxpix, lxpix, logxpx, shifty, shiftx;
0126   int nclusx, nclusy;
0127   int nybin, ycbin, nxbin, xcbin, minbinj, minbink;
0128   int deltaj, jmin, jmax, kmin, kmax, km, fxbin, lxbin, fybin, lybin, djy, djx;
0129   float sythr, sxthr, delta, sigma, sigavg, pseudopix, xsize, ysize, qscale, lanpar[2][5];
0130   float ss2, ssa, sa2, rat, fq, qtotal, qpixel, qavg;
0131   float x1p, x2p, y1p, y2p, deltachi2;
0132   float originx, originy, bias, maxpix, minmax;
0133   double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
0134   double hchi2, hndof, sigmal1, sigmal2, mpv1, mpv2, arg1, arg2, q05, like, loglike1, loglike2;
0135   double prvav, mpv, sigmaQ, kappa, xvav;
0136   float ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
0137   float ysig2[BYSIZE], xsig2[BXSIZE];
0138   float yw2[BYSIZE], xw2[BXSIZE], ysw[BYSIZE], xsw[BXSIZE];
0139   float cluster2[BXM2][BYM2], temp2d1[BXM2][BYM2], temp2d2[BXM2][BYM2];
0140   bool yd[BYSIZE], xd[BXSIZE], anyyd, anyxd, any2dfail;
0141   const float sqrt2x = {3.0000}, sqrt2y = {1.7000};
0142   const float sqrt12 = {3.4641};
0143   const float probmin = {1.110223e-16};
0144   const float prob2Qmin = {1.e-5};
0145   std::pair<int, int> pixel;
0146 
0147   //    bool SiPixelTemplateSplit::SimpleTemplate2D(float cotalpha, float cotbeta, float xhit, float yhit, float thick, float lorxwidth, float lorywidth,
0148   //                          float qavg, std::vector<bool> ydouble, std::vector<bool> xdouble, float template2d[BXM2][BYM2]);
0149 
0150   // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
0151 
0152   const double mean1pix = {0.100}, chi21min = {0.160};
0153 
0154   // First, interpolate the template needed to analyze this cluster
0155   // check to see of the track direction is in the physical range of the loaded template
0156 
0157   if (!templ.interpolate(id, cotalpha, cotbeta)) {
0158     LOGDEBUG("SiPixelTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha
0159                                     << ", cot(beta) = " << cotbeta
0160                                     << " is not within the acceptance of template ID = " << id
0161                                     << ", no reconstruction performed" << ENDL;
0162     return 20;
0163   }
0164 
0165   // Get pixel dimensions from the template (to allow multiple detectors in the future)
0166 
0167   xsize = templ.xsize();
0168   ysize = templ.ysize();
0169 
0170   // Define size of pseudopixel
0171 
0172   pseudopix = templ.s50();
0173   //    q05 = 0.28*pseudopix;
0174   q05 = 0.05f * pseudopix;
0175 
0176   // Get charge scaling factor
0177 
0178   qscale = templ.qscale();
0179 
0180   // Make a local copy of the cluster container so that we can muck with it
0181 
0182   array_2d cluster = clust;
0183 
0184   // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
0185 
0186   if (cluster.num_dimensions() != 2) {
0187     LOGERROR("SiPixelTemplateReco") << "input cluster container (BOOST Multiarray) has wrong number of dimensions"
0188                                     << ENDL;
0189     return 3;
0190   }
0191   nclusx = (int)cluster.shape()[0];
0192   nclusy = (int)cluster.shape()[1];
0193   if (nclusx != (int)xdouble.size()) {
0194     LOGERROR("SiPixelTemplateReco") << "input cluster container x-size is not equal to double pixel flag container size"
0195                                     << ENDL;
0196     return 4;
0197   }
0198   if (nclusy != (int)ydouble.size()) {
0199     LOGERROR("SiPixelTemplateReco") << "input cluster container y-size is not equal to double pixel flag container size"
0200                                     << ENDL;
0201     return 5;
0202   }
0203 
0204   // enforce maximum size
0205 
0206   if (nclusx > TXSIZE) {
0207     nclusx = TXSIZE;
0208   }
0209   if (nclusy > TYSIZE) {
0210     nclusy = TYSIZE;
0211   }
0212 
0213   // First, rescale all pixel charges
0214 
0215   for (i = 0; i < nclusy; ++i) {
0216     for (j = 0; j < nclusx; ++j) {
0217       if (cluster[j][i] > 0) {
0218         cluster[j][i] *= qscale;
0219       }
0220     }
0221   }
0222 
0223   // Next, sum the total charge and "decapitate" big pixels
0224 
0225   qtotal = 0.f;
0226   minmax = 2.0f * templ.pixmax();
0227   for (i = 0; i < nclusy; ++i) {
0228     maxpix = minmax;
0229     if (ydouble[i]) {
0230       maxpix *= 2.f;
0231     }
0232     for (j = 0; j < nclusx; ++j) {
0233       qtotal += cluster[j][i];
0234       if (cluster[j][i] > maxpix) {
0235         cluster[j][i] = maxpix;
0236       }
0237     }
0238   }
0239 
0240   // Do the cluster repair here
0241 
0242   if (deadpix) {
0243     fypix = BYM3;
0244     lypix = -1;
0245     for (i = 0; i < nclusy; ++i) {
0246       ysum[i] = 0.f;
0247       // Do preliminary cluster projection in y
0248       for (j = 0; j < nclusx; ++j) {
0249         ysum[i] += cluster[j][i];
0250       }
0251       if (ysum[i] > 0.f) {
0252         // identify ends of cluster to determine what the missing charge should be
0253         if (i < fypix) {
0254           fypix = i;
0255         }
0256         if (i > lypix) {
0257           lypix = i;
0258         }
0259       }
0260     }
0261 
0262     // Now loop over dead pixel list and "fix" everything
0263 
0264     //First see if the cluster ends are redefined and that we have only one dead pixel per column
0265     int nyzero[TYSIZE]{};
0266     std::vector<std::pair<int, int> >::const_iterator zeroIter = zeropix.begin(), zeroEnd = zeropix.end();
0267     for (; zeroIter != zeroEnd; ++zeroIter) {
0268       i = zeroIter->second;
0269       if (i < 0 || i > TYSIZE - 1) {
0270         LOGERROR("SiPixelTemplateReco") << "dead pixel column y-index " << i << ", no reconstruction performed" << ENDL;
0271         return 11;
0272       }
0273 
0274       // count the number of dead pixels in each column
0275       ++nyzero[i];
0276       // allow them to redefine the cluster ends
0277       if (i < fypix) {
0278         fypix = i;
0279       }
0280       if (i > lypix) {
0281         lypix = i;
0282       }
0283     }
0284 
0285     nypix = lypix - fypix + 1;
0286 
0287     // 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
0288 
0289     for (zeroIter = zeropix.begin(); zeroIter != zeroEnd; ++zeroIter) {
0290       i = zeroIter->second;
0291       j = zeroIter->first;
0292       if (j < 0 || j > TXSIZE - 1) {
0293         LOGERROR("SiPixelTemplateReco") << "dead pixel column x-index " << j << ", no reconstruction performed" << ENDL;
0294         return 12;
0295       }
0296       if ((i == fypix || i == lypix) && nypix > 1) {
0297         maxpix = templ.symax() / 2.;
0298       } else {
0299         maxpix = templ.symax();
0300       }
0301       if (ydouble[i]) {
0302         maxpix *= 2.;
0303       }
0304       if (nyzero[i] > 0 && nyzero[i] < 3) {
0305         qpixel = (maxpix - ysum[i]) / (float)nyzero[i];
0306       } else {
0307         qpixel = 1.;
0308       }
0309       if (qpixel < 1.f) {
0310         qpixel = 1.f;
0311       }
0312       cluster[j][i] = qpixel;
0313       // Adjust the total cluster charge to reflect the charge of the "repaired" cluster
0314       qtotal += qpixel;
0315     }
0316     // End of cluster repair section
0317   }
0318 
0319   // Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container
0320 
0321   for (i = 0; i < BYSIZE; ++i) {
0322     ysum[i] = 0.f;
0323     yd[i] = false;
0324   }
0325   k = 0;
0326   anyyd = false;
0327   for (i = 0; i < nclusy; ++i) {
0328     for (j = 0; j < nclusx; ++j) {
0329       ysum[k] += cluster[j][i];
0330     }
0331 
0332     // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
0333 
0334     if (ydouble[i]) {
0335       ysum[k] /= 2.f;
0336       ysum[k + 1] = ysum[k];
0337       yd[k] = true;
0338       yd[k + 1] = false;
0339       k = k + 2;
0340       anyyd = true;
0341     } else {
0342       yd[k] = false;
0343       ++k;
0344     }
0345     if (k > BYM1) {
0346       break;
0347     }
0348   }
0349 
0350   // Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container
0351 
0352   for (i = 0; i < BXSIZE; ++i) {
0353     xsum[i] = 0.f;
0354     xd[i] = false;
0355   }
0356   k = 0;
0357   anyxd = false;
0358   for (j = 0; j < nclusx; ++j) {
0359     for (i = 0; i < nclusy; ++i) {
0360       xsum[k] += cluster[j][i];
0361     }
0362 
0363     // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
0364 
0365     if (xdouble[j]) {
0366       xsum[k] /= 2.f;
0367       xsum[k + 1] = xsum[k];
0368       xd[k] = true;
0369       xd[k + 1] = false;
0370       k = k + 2;
0371       anyxd = true;
0372     } else {
0373       xd[k] = false;
0374       ++k;
0375     }
0376     if (k > BXM1) {
0377       break;
0378     }
0379   }
0380 
0381   // next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx
0382 
0383   fypix = -1;
0384   nypix = 0;
0385   lypix = 0;
0386   logypx = 0;
0387   for (i = 0; i < BYSIZE; ++i) {
0388     if (ysum[i] > 0.) {
0389       if (fypix == -1) {
0390         fypix = i;
0391       }
0392       if (!yd[i]) {
0393         ysort[logypx] = ysum[i];
0394         ++logypx;
0395       }
0396       ++nypix;
0397       lypix = i;
0398     }
0399   }
0400 
0401   // Make sure cluster is continuous
0402 
0403   if ((lypix - fypix + 1) != nypix || nypix == 0) {
0404     LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster doesn't agree with number of pixels above threshold"
0405                                     << ENDL;
0406     if (theVerboseLevel > 2) {
0407       LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
0408       for (i = 0; i < BYSIZE - 1; ++i) {
0409         LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";
0410       }
0411       LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE - 1] << ENDL;
0412     }
0413 
0414     return 1;
0415   }
0416 
0417   // If cluster is longer than max template size, technique fails
0418 
0419   if (nypix > TYSIZE) {
0420     LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster is larger than maximum template size" << ENDL;
0421     if (theVerboseLevel > 2) {
0422       LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
0423       for (i = 0; i < BYSIZE - 1; ++i) {
0424         LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";
0425       }
0426       LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE - 1] << ENDL;
0427     }
0428 
0429     return 6;
0430   }
0431 
0432   // next, center the cluster on template center if necessary
0433 
0434   midpix = (fypix + lypix) / 2;
0435   //    shifty = BHY - midpix;
0436   shifty = templ.cytemp() - midpix;
0437   if (shifty > 0) {
0438     for (i = lypix; i >= fypix; --i) {
0439       ysum[i + shifty] = ysum[i];
0440       ysum[i] = 0.;
0441       yd[i + shifty] = yd[i];
0442       yd[i] = false;
0443     }
0444   } else if (shifty < 0) {
0445     for (i = fypix; i <= lypix; ++i) {
0446       ysum[i + shifty] = ysum[i];
0447       ysum[i] = 0.;
0448       yd[i + shifty] = yd[i];
0449       yd[i] = false;
0450     }
0451   }
0452   lypix += shifty;
0453   fypix += shifty;
0454 
0455   // If the cluster boundaries are OK, add pesudopixels, otherwise quit
0456 
0457   if (fypix > 1 && fypix < BYM2) {
0458     ysum[fypix - 1] = pseudopix;
0459     ysum[fypix - 2] = 0.2f * pseudopix;
0460   } else {
0461     return 8;
0462   }
0463   if (lypix > 1 && lypix < BYM2) {
0464     ysum[lypix + 1] = pseudopix;
0465     ysum[lypix + 2] = 0.2f * pseudopix;
0466   } else {
0467     return 8;
0468   }
0469 
0470   // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
0471 
0472   if (ydouble[0]) {
0473     originy = -0.5f;
0474   } else {
0475     originy = 0.f;
0476   }
0477 
0478   // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx
0479 
0480   fxpix = -1;
0481   nxpix = 0;
0482   lxpix = 0;
0483   logxpx = 0;
0484   for (i = 0; i < BXSIZE; ++i) {
0485     if (xsum[i] > 0.) {
0486       if (fxpix == -1) {
0487         fxpix = i;
0488       }
0489       if (!xd[i]) {
0490         xsort[logxpx] = xsum[i];
0491         ++logxpx;
0492       }
0493       ++nxpix;
0494       lxpix = i;
0495     }
0496   }
0497 
0498   // Make sure cluster is continuous
0499 
0500   if ((lxpix - fxpix + 1) != nxpix) {
0501     LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold"
0502                                     << ENDL;
0503     if (theVerboseLevel > 2) {
0504       LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
0505       for (i = 0; i < BXSIZE - 1; ++i) {
0506         LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";
0507       }
0508       LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE - 1] << ENDL;
0509     }
0510 
0511     return 2;
0512   }
0513 
0514   // If cluster is longer than max template size, technique fails
0515 
0516   if (nxpix > TXSIZE) {
0517     LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
0518     if (theVerboseLevel > 2) {
0519       LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
0520       for (i = 0; i < BXSIZE - 1; ++i) {
0521         LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";
0522       }
0523       LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE - 1] << ENDL;
0524     }
0525 
0526     return 7;
0527   }
0528 
0529   // next, center the cluster on template center if necessary
0530 
0531   midpix = (fxpix + lxpix) / 2;
0532   //    shiftx = BHX - midpix;
0533   shiftx = templ.cxtemp() - midpix;
0534   if (shiftx > 0) {
0535     for (i = lxpix; i >= fxpix; --i) {
0536       xsum[i + shiftx] = xsum[i];
0537       xsum[i] = 0.f;
0538       xd[i + shiftx] = xd[i];
0539       xd[i] = false;
0540     }
0541   } else if (shiftx < 0) {
0542     for (i = fxpix; i <= lxpix; ++i) {
0543       xsum[i + shiftx] = xsum[i];
0544       xsum[i] = 0.f;
0545       xd[i + shiftx] = xd[i];
0546       xd[i] = false;
0547     }
0548   }
0549   lxpix += shiftx;
0550   fxpix += shiftx;
0551 
0552   // If the cluster boundaries are OK, add pesudopixels, otherwise quit
0553 
0554   if (fxpix > 1 && fxpix < BXM2) {
0555     xsum[fxpix - 1] = pseudopix;
0556     xsum[fxpix - 2] = 0.2f * pseudopix;
0557   } else {
0558     return 9;
0559   }
0560   if (lxpix > 1 && lxpix < BXM2) {
0561     xsum[lxpix + 1] = pseudopix;
0562     xsum[lxpix + 2] = 0.2f * pseudopix;
0563   } else {
0564     return 9;
0565   }
0566 
0567   // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
0568 
0569   if (xdouble[0]) {
0570     originx = -0.5f;
0571   } else {
0572     originx = 0.f;
0573   }
0574 
0575   // uncertainty and final corrections depend upon total charge bin
0576 
0577   qavg = templ.qavg();
0578   fq = qtotal / qavg;
0579   if (fq > 3.0f) {
0580     binq = 0;
0581   } else {
0582     if (fq > 2.0f) {
0583       binq = 1;
0584     } else {
0585       if (fq > 1.70f) {
0586         binq = 2;
0587       } else {
0588         binq = 3;
0589       }
0590     }
0591   }
0592 
0593   // Calculate the Vavilov probability that the cluster charge is consistent with a merged cluster
0594 
0595   if (speed < 0) {
0596     templ.vavilov2_pars(mpv, sigmaQ, kappa);
0597 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0598     if ((sigmaQ <= 0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
0599       throw cms::Exception("DataCorrupt") << "SiPixelTemplateSplit::Vavilov parameters mpv/sigmaQ/kappa = " << mpv
0600                                           << "/" << sigmaQ << "/" << kappa << std::endl;
0601     }
0602 #else
0603     assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
0604 #endif
0605     xvav = ((double)qtotal - mpv) / sigmaQ;
0606     //  VVIObj is a private port of CERNLIB VVIDIS
0607     VVIObj vvidist(kappa);
0608     prvav = vvidist.fcn(xvav);
0609     prob2Q = 1. - prvav;
0610     if (prob2Q < prob2Qmin) {
0611       prob2Q = prob2Qmin;
0612     }
0613   } else {
0614     prob2Q = -1.f;
0615   }
0616 
0617   // Return the charge bin via the parameter list unless the charge is too small (then flag it)
0618 
0619   q2bin = binq;
0620   if (!deadpix && qtotal < 1.9f * templ.qmin()) {
0621     q2bin = 5;
0622   } else {
0623     if (!deadpix && qtotal < 1.9f * templ.qmin(1)) {
0624       q2bin = 4;
0625     }
0626   }
0627 
0628   if (theVerboseLevel > 9) {
0629     LOGDEBUG("SiPixelTemplateSplit") << "ID = " << id << " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta
0630                                      << " nclusx = " << nclusx << " nclusy = " << nclusy << ENDL;
0631   }
0632 
0633   // Next, generate the 3d y- and x-templates
0634 
0635   templ.ytemp3d_int(nypix, nybin);
0636 
0637   ycbin = nybin / 2;
0638 
0639   templ.xtemp3d_int(nxpix, nxbin);
0640 
0641   // retrieve the number of x-bins
0642 
0643   xcbin = nxbin / 2;
0644 
0645   // First, decide on chi^2 min search parameters
0646 
0647 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0648   if (speed < -1 || speed > 2) {
0649     throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::PixelTempReco2D called with illegal speed = " << speed
0650                                         << std::endl;
0651   }
0652 #else
0653   assert(speed >= -1 && speed < 3);
0654 #endif
0655   fybin = 0;
0656   lybin = nybin - 1;
0657   fxbin = 0;
0658   lxbin = nxbin - 1;
0659   djy = 1;
0660   djx = 1;
0661   if (speed > 0) {
0662     djy = 2;
0663     djx = 2;
0664     if (speed > 1) {
0665       if (!anyyd) {
0666         djy = 4;
0667       }
0668       if (!anyxd) {
0669         djx = 4;
0670       }
0671     }
0672   }
0673 
0674   if (theVerboseLevel > 9) {
0675     LOGDEBUG("SiPixelTemplateReco") << "fypix " << fypix << " lypix = " << lypix << " fybin = " << fybin
0676                                     << " lybin = " << lybin << " djy = " << djy << " logypx = " << logypx << ENDL;
0677     LOGDEBUG("SiPixelTemplateReco") << "fxpix " << fxpix << " lxpix = " << lxpix << " fxbin = " << fxbin
0678                                     << " lxbin = " << lxbin << " djx = " << djx << " logxpx = " << logxpx << ENDL;
0679   }
0680 
0681   // Do the y-reconstruction first
0682 
0683   // Define the maximum signal to allow before de-weighting a pixel
0684 
0685   sythr = 2.1f * (templ.symax());
0686 
0687   // Make sure that there will be at least two pixels that are not de-weighted
0688 
0689   std::sort(&ysort[0], &ysort[logypx]);
0690   if (logypx == 1) {
0691     sythr = 1.01f * ysort[0];
0692   } else {
0693     if (ysort[1] > sythr) {
0694       sythr = 1.01f * ysort[1];
0695     }
0696   }
0697 
0698   // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
0699 
0700   //    for(i=0; i<BYSIZE; ++i) { ysig2[i] = 0.;}
0701   templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
0702 
0703   // Find the template bin that minimizes the Chi^2
0704 
0705   chi2ymin = 1.e15;
0706   ss2 = 0.f;
0707   for (i = fypix - 2; i <= lypix + 2; ++i) {
0708     yw2[i] = 1.f / ysig2[i];
0709     ysw[i] = ysum[i] * yw2[i];
0710     ss2 += ysum[i] * ysw[i];
0711   }
0712   minbinj = -1;
0713   minbink = -1;
0714   deltaj = djy;
0715   jmin = fybin;
0716   jmax = lybin;
0717   kmin = fybin;
0718   kmax = lybin;
0719   std::vector<float> ytemp(BYSIZE);
0720   while (deltaj > 0) {
0721     for (j = jmin; j < jmax; j += deltaj) {
0722       km = std::min(kmax, j);
0723       for (k = kmin; k <= km; k += deltaj) {
0724         // Get the template for this set of indices
0725 
0726         templ.ytemp3d(j, k, ytemp);
0727 
0728         // Modify the template if double pixels are present
0729 
0730         if (nypix > logypx) {
0731           i = fypix;
0732           while (i < lypix) {
0733             if (yd[i] && !yd[i + 1]) {
0734               // Sum the adjacent cells and put the average signal in both
0735 
0736               sigavg = (ytemp[i] + ytemp[i + 1]) / 2.f;
0737               ytemp[i] = sigavg;
0738               ytemp[i + 1] = sigavg;
0739               i += 2;
0740             } else {
0741               ++i;
0742             }
0743           }
0744         }
0745         ssa = 0.f;
0746         sa2 = 0.f;
0747         for (i = fypix - 2; i <= lypix + 2; ++i) {
0748           ssa += ysw[i] * ytemp[i];
0749           sa2 += ytemp[i] * ytemp[i] * yw2[i];
0750         }
0751         rat = ssa / ss2;
0752         if (rat <= 0.) {
0753           LOGERROR("SiPixelTemplateSplit") << "illegal chi2ymin normalization = " << rat << ENDL;
0754           rat = 1.;
0755         }
0756         chi2y = ss2 - 2.f * ssa / rat + sa2 / (rat * rat);
0757         if (chi2y < chi2ymin) {
0758           chi2ymin = chi2y;
0759           minbinj = j;
0760           minbink = k;
0761         }
0762       }
0763     }
0764     deltaj /= 2;
0765     if (minbinj > fybin) {
0766       jmin = minbinj - deltaj;
0767     } else {
0768       jmin = fybin;
0769     }
0770     if (minbinj < lybin) {
0771       jmax = minbinj + deltaj;
0772     } else {
0773       jmax = lybin;
0774     }
0775     if (minbink > fybin) {
0776       kmin = minbink - deltaj;
0777     } else {
0778       kmin = fybin;
0779     }
0780     if (minbink < lybin) {
0781       kmax = minbink + deltaj;
0782     } else {
0783       kmax = lybin;
0784     }
0785   }
0786 
0787   if (theVerboseLevel > 9) {
0788     LOGDEBUG("SiPixelTemplateReco") << "minbins " << minbinj << "," << minbink << " chi2ymin = " << chi2ymin << ENDL;
0789   }
0790 
0791   // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
0792 
0793   if (logypx == 1) {
0794     if (nypix == 1) {
0795       delta = templ.dyone();
0796       sigma = templ.syone();
0797     } else {
0798       delta = templ.dytwo();
0799       sigma = templ.sytwo();
0800     }
0801 
0802     yrec1 = 0.5f * (fypix + lypix - 2 * shifty + 2.f * originy) * ysize - delta;
0803     yrec2 = yrec1;
0804 
0805     if (sigma <= 0.) {
0806       sigmay = ysize / sqrt12;
0807     } else {
0808       sigmay = sigma;
0809     }
0810 
0811     // Do probability calculation for one-pixel clusters
0812 
0813     chi21max = fmax(chi21min, (double)templ.chi2yminone());
0814     chi2ymin -= chi21max;
0815     if (chi2ymin < 0.) {
0816       chi2ymin = 0.;
0817     }
0818     //     prob2y = gsl_cdf_chisq_Q(chi2ymin, mean1pix);
0819     meany = fmax(mean1pix, (double)templ.chi2yavgone());
0820     hchi2 = chi2ymin / 2.;
0821     hndof = meany / 2.;
0822     prob2y = 1. - TMath::Gamma(hndof, hchi2);
0823 
0824   } else {
0825     // For cluster > 1 pix, use chi^2 minimm to recontruct the two y-positions
0826 
0827     // at small eta, the templates won't actually work on two pixel y-clusters so just return the pixel centers
0828 
0829     if (logypx == 2 && fabsf(cotbeta) < 0.25f) {
0830       switch (nypix) {
0831         case 2:
0832           //  Both pixels are small
0833           yrec1 = (fypix - shifty + originy) * ysize;
0834           yrec2 = (lypix - shifty + originy) * ysize;
0835           sigmay = ysize / sqrt12;
0836           break;
0837         case 3:
0838           //  One big pixel and one small pixel
0839           if (yd[fypix]) {
0840             yrec1 = (fypix + 0.5f - shifty + originy) * ysize;
0841             yrec2 = (lypix - shifty + originy) * ysize;
0842             sigmay = ysize / sqrt12;
0843           } else {
0844             yrec1 = (fypix - shifty + originy) * ysize;
0845             yrec2 = (lypix - 0.5f - shifty + originy) * ysize;
0846             sigmay = 1.5f * ysize / sqrt12;
0847           }
0848           break;
0849         case 4:
0850           //  Two big pixels
0851           yrec1 = (fypix + 0.5f - shifty + originy) * ysize;
0852           yrec2 = (lypix - 0.5f - shifty + originy) * ysize;
0853           sigmay = 2.f * ysize / sqrt12;
0854           break;
0855         default:
0856           //  Something is screwy ...
0857           LOGERROR("SiPixelTemplateReco")
0858               << "weird problem: logical y-pixels = " << logypx << ", total ysize in normal pixels = " << nypix << ENDL;
0859           return 10;
0860       }
0861     } else {
0862       // uncertainty and final correction depend upon charge bin
0863 
0864       bias = templ.yavgc2m(binq);
0865       yrec1 = (0.125f * (minbink - ycbin) + BHY - (float)shifty + originy) * ysize - bias;
0866       yrec2 = (0.125f * (minbinj - ycbin) + BHY - (float)shifty + originy) * ysize - bias;
0867       sigmay = sqrt2y * templ.yrmsc2m(binq);
0868     }
0869 
0870     // Do goodness of fit test in y
0871 
0872     if (chi2ymin < 0.0) {
0873       chi2ymin = 0.0;
0874     }
0875     meany = templ.chi2yavgc2m(binq);
0876     if (meany < 0.01) {
0877       meany = 0.01;
0878     }
0879     // gsl function that calculates the chi^2 tail prob for non-integral dof
0880     //     prob2y = gsl_cdf_chisq_Q(chi2y, meany);
0881     //     prob2y = ROOT::Math::chisquared_cdf_c(chi2y, meany);
0882     hchi2 = chi2ymin / 2.;
0883     hndof = meany / 2.;
0884     prob2y = 1. - TMath::Gamma(hndof, hchi2);
0885   }
0886 
0887   // Do the x-reconstruction next
0888 
0889   // Define the maximum signal to allow before de-weighting a pixel
0890 
0891   sxthr = 2.1f * templ.sxmax();
0892 
0893   // Make sure that there will be at least two pixels that are not de-weighted
0894   std::sort(&xsort[0], &xsort[logxpx]);
0895   if (logxpx == 1) {
0896     sxthr = 1.01f * xsort[0];
0897   } else {
0898     if (xsort[1] > sxthr) {
0899       sxthr = 1.01f * xsort[1];
0900     }
0901   }
0902 
0903   // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
0904 
0905   //    for(i=0; i<BYSIZE; ++i) { xsig2[i] = 0.; }
0906   templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
0907 
0908   // Find the template bin that minimizes the Chi^2
0909 
0910   chi2xmin = 1.e15;
0911   ss2 = 0.f;
0912   for (i = fxpix - 2; i <= lxpix + 2; ++i) {
0913     xw2[i] = 1.f / xsig2[i];
0914     xsw[i] = xsum[i] * xw2[i];
0915     ss2 += xsum[i] * xsw[i];
0916   }
0917   minbinj = -1;
0918   minbink = -1;
0919   deltaj = djx;
0920   jmin = fxbin;
0921   jmax = lxbin;
0922   kmin = fxbin;
0923   kmax = lxbin;
0924   std::vector<float> xtemp(BXSIZE);
0925   while (deltaj > 0) {
0926     for (j = jmin; j < jmax; j += deltaj) {
0927       km = std::min(kmax, j);
0928       for (k = kmin; k <= km; k += deltaj) {
0929         // Get the template for this set of indices
0930 
0931         templ.xtemp3d(j, k, xtemp);
0932 
0933         // Modify the template if double pixels are present
0934 
0935         if (nxpix > logxpx) {
0936           i = fxpix;
0937           while (i < lxpix) {
0938             if (xd[i] && !xd[i + 1]) {
0939               // Sum the adjacent cells and put the average signal in both
0940 
0941               sigavg = (xtemp[i] + xtemp[i + 1]) / 2.f;
0942               xtemp[i] = sigavg;
0943               xtemp[i + 1] = sigavg;
0944               i += 2;
0945             } else {
0946               ++i;
0947             }
0948           }
0949         }
0950         ssa = 0.f;
0951         sa2 = 0.f;
0952         for (i = fxpix - 2; i <= lxpix + 2; ++i) {
0953           ssa += xsw[i] * xtemp[i];
0954           sa2 += xtemp[i] * xtemp[i] * xw2[i];
0955         }
0956         rat = ssa / ss2;
0957         if (rat <= 0.f) {
0958           LOGERROR("SiPixelTemplateSplit") << "illegal chi2xmin normalization = " << rat << ENDL;
0959           rat = 1.;
0960         }
0961         chi2x = ss2 - 2.f * ssa / rat + sa2 / (rat * rat);
0962         if (chi2x < chi2xmin) {
0963           chi2xmin = chi2x;
0964           minbinj = j;
0965           minbink = k;
0966         }
0967       }
0968     }
0969     deltaj /= 2;
0970     if (minbinj > fxbin) {
0971       jmin = minbinj - deltaj;
0972     } else {
0973       jmin = fxbin;
0974     }
0975     if (minbinj < lxbin) {
0976       jmax = minbinj + deltaj;
0977     } else {
0978       jmax = lxbin;
0979     }
0980     if (minbink > fxbin) {
0981       kmin = minbink - deltaj;
0982     } else {
0983       kmin = fxbin;
0984     }
0985     if (minbink < lxbin) {
0986       kmax = minbink + deltaj;
0987     } else {
0988       kmax = lxbin;
0989     }
0990   }
0991 
0992   if (theVerboseLevel > 9) {
0993     LOGDEBUG("SiPixelTemplateSplit") << "minbinj/minbink " << minbinj << "/" << minbink << " chi2xmin = " << chi2xmin
0994                                      << ENDL;
0995   }
0996 
0997   // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
0998 
0999   if (logxpx == 1) {
1000     if (nxpix == 1) {
1001       delta = templ.dxone();
1002       sigma = templ.sxone();
1003     } else {
1004       delta = templ.dxtwo();
1005       sigma = templ.sxtwo();
1006     }
1007     xrec1 = 0.5f * (fxpix + lxpix - 2 * shiftx + 2.f * originx) * xsize - delta;
1008     xrec2 = xrec1;
1009     if (sigma <= 0.f) {
1010       sigmax = xsize / sqrt12;
1011     } else {
1012       sigmax = sigma;
1013     }
1014 
1015     // Do probability calculation for one-pixel clusters
1016 
1017     chi21max = fmax(chi21min, (double)templ.chi2xminone());
1018     chi2xmin -= chi21max;
1019     if (chi2xmin < 0.) {
1020       chi2xmin = 0.;
1021     }
1022     meanx = fmax(mean1pix, (double)templ.chi2xavgone());
1023     hchi2 = chi2xmin / 2.;
1024     hndof = meanx / 2.;
1025     prob2x = 1. - TMath::Gamma(hndof, hchi2);
1026 
1027   } else {
1028     // For cluster > 1 pix, use chi^2 minimm to recontruct the two x-positions
1029 
1030     // uncertainty and final correction depend upon charge bin
1031 
1032     bias = templ.xavgc2m(binq);
1033     xrec1 = (0.125f * (minbink - xcbin) + BHX - (float)shiftx + originx) * xsize - bias;
1034     xrec2 = (0.125f * (minbinj - xcbin) + BHX - (float)shiftx + originx) * xsize - bias;
1035     sigmax = sqrt2x * templ.xrmsc2m(binq);
1036 
1037     // Do goodness of fit test in y
1038 
1039     if (chi2xmin < 0.0) {
1040       chi2xmin = 0.0;
1041     }
1042     meanx = templ.chi2xavgc2m(binq);
1043     if (meanx < 0.01) {
1044       meanx = 0.01;
1045     }
1046     hchi2 = chi2xmin / 2.;
1047     hndof = meanx / 2.;
1048     prob2x = 1. - TMath::Gamma(hndof, hchi2);
1049   }
1050 
1051   //  Don't return exact zeros for the probability
1052 
1053   if (prob2y < probmin) {
1054     prob2y = probmin;
1055   }
1056   if (prob2x < probmin) {
1057     prob2x = probmin;
1058   }
1059 
1060   //  New code: resolve the ambiguity if resolve is set to true
1061 
1062   dchisq = 0.;
1063   if (resolve) {
1064     //  First copy the unexpanded cluster into a new working buffer
1065 
1066     for (i = 0; i < BYM2; ++i) {
1067       for (j = 0; j < BXM2; ++j) {
1068         if ((i > 0 && i < BYM3) && (j > 0 && j < BXM3)) {
1069           cluster2[j][i] = qscale * clust[j - 1][i - 1];
1070         } else {
1071           cluster2[j][i] = 0.f;
1072         }
1073       }
1074     }
1075 
1076     //  Next, redefine the local coordinates to start at the lower LH corner of pixel[1][1];
1077 
1078     if (xdouble[0]) {
1079       x1p = xrec1 + xsize;
1080       x2p = xrec2 + xsize;
1081     } else {
1082       x1p = xrec1 + xsize / 2.f;
1083       x2p = xrec2 + xsize / 2.f;
1084     }
1085 
1086     if (ydouble[0]) {
1087       y1p = yrec1 + ysize;
1088       y2p = yrec2 + ysize;
1089     } else {
1090       y1p = yrec1 + ysize / 2.f;
1091       y2p = yrec2 + ysize / 2.f;
1092     }
1093 
1094     // Next, calculate 2-d templates for the (x1,y1)+(x2,y2) and the (x1,y2)+(x2,y1) hypotheses
1095 
1096     //  First zero the two 2-d hypotheses
1097 
1098     for (i = 0; i < BYM2; ++i) {
1099       for (j = 0; j < BXM2; ++j) {
1100         temp2d1[j][i] = 0.f;
1101         temp2d2[j][i] = 0.f;
1102       }
1103     }
1104 
1105     // Add the two hits in the first hypothesis
1106 
1107     any2dfail = templ2D.xytemp(id, cotalpha, cotbeta, x1p, y1p, xdouble, ydouble, temp2d1) &&
1108                 templ2D.xytemp(id, cotalpha, cotbeta, x2p, y2p, xdouble, ydouble, temp2d1);
1109 
1110     // And then the second hypothesis
1111 
1112     any2dfail = any2dfail && templ2D.xytemp(id, cotalpha, cotbeta, x1p, y2p, xdouble, ydouble, temp2d2);
1113 
1114     any2dfail = any2dfail && templ2D.xytemp(id, cotalpha, cotbeta, x2p, y1p, xdouble, ydouble, temp2d2);
1115 
1116     // If any of these have failed, use the simple templates instead
1117 
1118     if (!any2dfail) {
1119       //  Rezero the two 2-d hypotheses
1120 
1121       for (i = 0; i < BYM2; ++i) {
1122         for (j = 0; j < BXM2; ++j) {
1123           temp2d1[j][i] = 0.f;
1124           temp2d2[j][i] = 0.f;
1125         }
1126       }
1127 
1128       // Add the two hits in the first hypothesis
1129 
1130       if (!templ.simpletemplate2D(x1p, y1p, xdouble, ydouble, temp2d1)) {
1131         return 1;
1132       }
1133 
1134       if (!templ.simpletemplate2D(x2p, y2p, xdouble, ydouble, temp2d1)) {
1135         return 1;
1136       }
1137 
1138       // And then the second hypothesis
1139 
1140       if (!templ.simpletemplate2D(x1p, y2p, xdouble, ydouble, temp2d2)) {
1141         return 1;
1142       }
1143 
1144       if (!templ.simpletemplate2D(x2p, y1p, xdouble, ydouble, temp2d2)) {
1145         return 1;
1146       }
1147     }
1148 
1149     // Keep lists of pixels and nearest neighbors
1150 
1151     std::list<std::pair<int, int> > pixellst;
1152 
1153     // Loop through the array and find non-zero elements
1154 
1155     for (i = 0; i < BYM2; ++i) {
1156       for (j = 0; j < BXM2; ++j) {
1157         if (cluster2[j][i] > 0.f || temp2d1[j][i] > 0.f || temp2d2[j][i] > 0.f) {
1158           pixel.first = j;
1159           pixel.second = i;
1160           pixellst.push_back(pixel);
1161         }
1162       }
1163     }
1164 
1165     // Now calculate the product of Landau probabilities (alpha probability)
1166 
1167     templ2D.landau_par(lanpar);
1168     loglike1 = 0.;
1169     loglike2 = 0.;
1170 
1171     // Now, for each neighbor, match it again the pixel list.
1172     // If found, delete it from the neighbor list
1173 
1174     std::list<std::pair<int, int> >::const_iterator pixIter, pixEnd;
1175     pixIter = pixellst.begin();
1176     pixEnd = pixellst.end();
1177     for (; pixIter != pixEnd; ++pixIter) {
1178       j = pixIter->first;
1179       i = pixIter->second;
1180       if ((i < BHY && cotbeta > 0.) || (i >= BHY && cotbeta < 0.)) {
1181         lparm = 0;
1182       } else {
1183         lparm = 1;
1184       }
1185       mpv1 = lanpar[lparm][0] + lanpar[lparm][1] * temp2d1[j][i];
1186       sigmal1 = lanpar[lparm][2] * mpv1;
1187       sigmal1 = sqrt(sigmal1 * sigmal1 + lanpar[lparm][3] * lanpar[lparm][3]);
1188       if (sigmal1 < q05)
1189         sigmal1 = q05;
1190       arg1 = (cluster2[j][i] - mpv1) / sigmal1 - 0.22278;
1191       mpv2 = lanpar[lparm][0] + lanpar[lparm][1] * temp2d2[j][i];
1192       sigmal2 = lanpar[lparm][2] * mpv2;
1193       sigmal2 = sqrt(sigmal2 * sigmal2 + lanpar[lparm][3] * lanpar[lparm][3]);
1194       if (sigmal2 < q05)
1195         sigmal2 = q05;
1196       arg2 = (cluster2[j][i] - mpv2) / sigmal2 - 0.22278;
1197       //            like = ROOT::Math::landau_pdf(arg1)/sigmal1;
1198       like = ROOT::Math::landau_pdf(arg1);
1199       if (like < 1.e-30)
1200         like = 1.e-30;
1201       loglike1 += log(like);
1202       //            like = ROOT::Math::landau_pdf(arg2)/sigmal2;
1203       like = ROOT::Math::landau_pdf(arg2);
1204       if (like < 1.e-30)
1205         like = 1.e-30;
1206       loglike2 += log(like);
1207     }
1208 
1209     // Calculate chisquare difference for the two hypotheses 9don't multiply by 2 for less inconvenient scaling
1210 
1211     deltachi2 = loglike1 - loglike2;
1212 
1213     if (deltachi2 < 0.) {
1214       // Flip the x1 and x2
1215 
1216       x1p = xrec1;
1217       xrec1 = xrec2;
1218       xrec2 = x1p;
1219     }
1220 
1221     // Return a positive definite value
1222 
1223     dchisq = fabs(deltachi2);
1224   }
1225 
1226   return 0;
1227 }  // PixelTempSplit
1228 
1229 // *************************************************************************************************************************************
1230 //! Reconstruct the best estimate of the hit positions for pixel clusters.
1231 //! \param         id - (input) identifier of the template to use
1232 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1233 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1234 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
1235 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
1236 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
1237 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
1238 //! \param      templ - (input) the template used in the reconstruction
1239 //! \param      yrec1 - (output) best estimate of first y-coordinate of hit in microns
1240 //! \param      yrec2 - (output) best estimate of second y-coordinate of hit in microns
1241 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
1242 //! \param      prob2y - (output) probability describing goodness-of-fit for y-reco
1243 //! \param      xrec1 - (output) best estimate of first x-coordinate of hit in microns
1244 //! \param      xrec2 - (output) best estimate of second x-coordinate of hit in microns
1245 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
1246 //! \param      prob2x - (output) probability describing goodness-of-fit for x-reco
1247 //! \param       q2bin - (output) index (0-4) describing the charge of the cluster
1248 //!                     [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]
1249 //! \param      prob2Q - (output) probability that the cluster charge is compatible with a 2-hit merging
1250 //! \param    resolve - (input) use 2-D chisquare to resolve the 2-fold x-y association ambiguity (costs a factor of two in speed)
1251 //! \param     dchisq - (output) the delta chisquare estimator used to break the degeneracy (0 means no discrimination, larger than 0.1 is good)
1252 // *************************************************************************************************************************************
1253 
1254 int SiPixelTemplateSplit::PixelTempSplit(int id,
1255                                          float cotalpha,
1256                                          float cotbeta,
1257                                          array_2d& cluster,
1258                                          std::vector<bool>& ydouble,
1259                                          std::vector<bool>& xdouble,
1260                                          SiPixelTemplate& templ,
1261                                          float& yrec1,
1262                                          float& yrec2,
1263                                          float& sigmay,
1264                                          float& prob2y,
1265                                          float& xrec1,
1266                                          float& xrec2,
1267                                          float& sigmax,
1268                                          float& prob2x,
1269                                          int& q2bin,
1270                                          float& prob2Q,
1271                                          bool resolve,
1272                                          int speed,
1273                                          float& dchisq,
1274                                          SiPixelTemplate2D& templ2D) {
1275   // Local variables
1276   const bool deadpix = false;
1277   std::vector<std::pair<int, int> > zeropix;
1278 
1279   return SiPixelTemplateSplit::PixelTempSplit(id,
1280                                               cotalpha,
1281                                               cotbeta,
1282                                               cluster,
1283                                               ydouble,
1284                                               xdouble,
1285                                               templ,
1286                                               yrec1,
1287                                               yrec2,
1288                                               sigmay,
1289                                               prob2y,
1290                                               xrec1,
1291                                               xrec2,
1292                                               sigmax,
1293                                               prob2x,
1294                                               q2bin,
1295                                               prob2Q,
1296                                               resolve,
1297                                               speed,
1298                                               dchisq,
1299                                               deadpix,
1300                                               zeropix,
1301                                               templ2D);
1302 
1303 }  // PixelTempSplit
1304 
1305 // *************************************************************************************************************************************
1306 //! Reconstruct the best estimate of the hit positions for pixel clusters.
1307 //! \param         id - (input) identifier of the template to use
1308 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1309 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1310 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
1311 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
1312 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
1313 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
1314 //! \param      templ - (input) the template used in the reconstruction
1315 //! \param      yrec1 - (output) best estimate of first y-coordinate of hit in microns
1316 //! \param      yrec2 - (output) best estimate of second y-coordinate of hit in microns
1317 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
1318 //! \param      prob2y - (output) probability describing goodness-of-fit for y-reco
1319 //! \param      xrec1 - (output) best estimate of first x-coordinate of hit in microns
1320 //! \param      xrec2 - (output) best estimate of second x-coordinate of hit in microns
1321 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
1322 //! \param      prob2x - (output) probability describing goodness-of-fit for x-reco
1323 //! \param       q2bin - (output) index (0-4) describing the charge of the cluster
1324 //!                     [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]
1325 //! \param      prob2Q - (output) probability that the cluster charge is compatible with a 2-hit merging
1326 //! \param    resolve - (input) use 2-D chisquare to resolve the 2-fold x-y association ambiguity (costs a factor of two in speed)
1327 //! \param     dchisq - (output) the delta chisquare estimator used to break the degeneracy (0 means no discrimination, larger than 0.1 is good)
1328 // *************************************************************************************************************************************
1329 
1330 int SiPixelTemplateSplit::PixelTempSplit(int id,
1331                                          float cotalpha,
1332                                          float cotbeta,
1333                                          array_2d& cluster,
1334                                          std::vector<bool>& ydouble,
1335                                          std::vector<bool>& xdouble,
1336                                          SiPixelTemplate& templ,
1337                                          float& yrec1,
1338                                          float& yrec2,
1339                                          float& sigmay,
1340                                          float& prob2y,
1341                                          float& xrec1,
1342                                          float& xrec2,
1343                                          float& sigmax,
1344                                          float& prob2x,
1345                                          int& q2bin,
1346                                          float& prob2Q,
1347                                          bool resolve,
1348                                          float& dchisq,
1349                                          SiPixelTemplate2D& templ2D) {
1350   // Local variables
1351   const bool deadpix = false;
1352   std::vector<std::pair<int, int> > zeropix;
1353   const int speed = 1;
1354 
1355   return SiPixelTemplateSplit::PixelTempSplit(id,
1356                                               cotalpha,
1357                                               cotbeta,
1358                                               cluster,
1359                                               ydouble,
1360                                               xdouble,
1361                                               templ,
1362                                               yrec1,
1363                                               yrec2,
1364                                               sigmay,
1365                                               prob2y,
1366                                               xrec1,
1367                                               xrec2,
1368                                               sigmax,
1369                                               prob2x,
1370                                               q2bin,
1371                                               prob2Q,
1372                                               resolve,
1373                                               speed,
1374                                               dchisq,
1375                                               deadpix,
1376                                               zeropix,
1377                                               templ2D);
1378 
1379 }  // PixelTempSplit
1380 
1381 // *************************************************************************************************************************************
1382 //! Reconstruct the best estimate of the hit positions for pixel clusters. Legacy interface to older code.
1383 //! \param         id - (input) identifier of the template to use
1384 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1385 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1386 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
1387 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
1388 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
1389 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
1390 //! \param      templ - (input) the template used in the reconstruction
1391 //! \param      yrec1 - (output) best estimate of first y-coordinate of hit in microns
1392 //! \param      yrec2 - (output) best estimate of second y-coordinate of hit in microns
1393 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
1394 //! \param      prob2y - (output) probability describing goodness-of-fit for y-reco
1395 //! \param      xrec1 - (output) best estimate of first x-coordinate of hit in microns
1396 //! \param      xrec2 - (output) best estimate of second x-coordinate of hit in microns
1397 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
1398 //! \param      prob2x - (output) probability describing goodness-of-fit for x-reco
1399 //! \param       q2bin - (output) index (0-4) describing the charge of the cluster
1400 //!                     [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]
1401 //! \param      prob2Q - (output) probability that the cluster charge is compatible with a 2-hit merging
1402 // *************************************************************************************************************************************
1403 
1404 int SiPixelTemplateSplit::PixelTempSplit(int id,
1405                                          float cotalpha,
1406                                          float cotbeta,
1407                                          array_2d& cluster,
1408                                          std::vector<bool>& ydouble,
1409                                          std::vector<bool>& xdouble,
1410                                          SiPixelTemplate& templ,
1411                                          float& yrec1,
1412                                          float& yrec2,
1413                                          float& sigmay,
1414                                          float& prob2y,
1415                                          float& xrec1,
1416                                          float& xrec2,
1417                                          float& sigmax,
1418                                          float& prob2x,
1419                                          int& q2bin,
1420                                          float& prob2Q,
1421                                          SiPixelTemplate2D& templ2D) {
1422   // Local variables
1423   const bool deadpix = false;
1424   const bool resolve = true;
1425   float dchisq;
1426   std::vector<std::pair<int, int> > zeropix;
1427   const int speed = 1;
1428 
1429   return SiPixelTemplateSplit::PixelTempSplit(id,
1430                                               cotalpha,
1431                                               cotbeta,
1432                                               cluster,
1433                                               ydouble,
1434                                               xdouble,
1435                                               templ,
1436                                               yrec1,
1437                                               yrec2,
1438                                               sigmay,
1439                                               prob2y,
1440                                               xrec1,
1441                                               xrec2,
1442                                               sigmax,
1443                                               prob2x,
1444                                               q2bin,
1445                                               prob2Q,
1446                                               resolve,
1447                                               speed,
1448                                               dchisq,
1449                                               deadpix,
1450                                               zeropix,
1451                                               templ2D);
1452 
1453 }  // PixelTempSplit
1454 
1455 // *************************************************************************************************************************************
1456 //! Reconstruct the best estimate of the hit positions for pixel clusters. Legacy interface to older code.
1457 //! \param         id - (input) identifier of the template to use
1458 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1459 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1460 //! \param    cluster - (input) boost multi_array container of 7x21 array of pixel signals,
1461 //!           origin of local coords (0,0) at center of pixel cluster[0][0].
1462 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel
1463 //! \param    xdouble - (input) STL vector of 7 element array to flag a double-pixel
1464 //! \param      templ - (input) the template used in the reconstruction
1465 //! \param      yrec1 - (output) best estimate of first y-coordinate of hit in microns
1466 //! \param      yrec2 - (output) best estimate of second y-coordinate of hit in microns
1467 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns
1468 //! \param      prob2y - (output) probability describing goodness-of-fit for y-reco
1469 //! \param      xrec1 - (output) best estimate of first x-coordinate of hit in microns
1470 //! \param      xrec2 - (output) best estimate of second x-coordinate of hit in microns
1471 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns
1472 //! \param      prob2x - (output) probability describing goodness-of-fit for x-reco
1473 //! \param       q2bin - (output) index (0-4) describing the charge of the cluster
1474 //!                     [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]
1475 // *************************************************************************************************************************************
1476 
1477 int SiPixelTemplateSplit::PixelTempSplit(int id,
1478                                          float cotalpha,
1479                                          float cotbeta,
1480                                          array_2d& cluster,
1481                                          std::vector<bool>& ydouble,
1482                                          std::vector<bool>& xdouble,
1483                                          SiPixelTemplate& templ,
1484                                          float& yrec1,
1485                                          float& yrec2,
1486                                          float& sigmay,
1487                                          float& prob2y,
1488                                          float& xrec1,
1489                                          float& xrec2,
1490                                          float& sigmax,
1491                                          float& prob2x,
1492                                          int& q2bin,
1493                                          SiPixelTemplate2D& templ2D) {
1494   // Local variables
1495   const bool deadpix = false;
1496   const bool resolve = true;
1497   float dchisq, prob2Q;
1498   std::vector<std::pair<int, int> > zeropix;
1499   const int speed = 1;
1500 
1501   return SiPixelTemplateSplit::PixelTempSplit(id,
1502                                               cotalpha,
1503                                               cotbeta,
1504                                               cluster,
1505                                               ydouble,
1506                                               xdouble,
1507                                               templ,
1508                                               yrec1,
1509                                               yrec2,
1510                                               sigmay,
1511                                               prob2y,
1512                                               xrec1,
1513                                               xrec2,
1514                                               sigmax,
1515                                               prob2x,
1516                                               q2bin,
1517                                               prob2Q,
1518                                               resolve,
1519                                               speed,
1520                                               dchisq,
1521                                               deadpix,
1522                                               zeropix,
1523                                               templ2D);
1524 
1525 }  // PixelTempSplit