/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateSplit.cc

Go to the documentation of this file.

//
//  SiPixelTemplateSplit.cc (Version 2.30)
//
//  Procedure to fit two templates (same angle hypotheses) to a single cluster
//  Return two x- and two y-coordinates for the cluster
//
//  Created by Morris Swartz on 04/10/08.
//
//  Incorporate "cluster repair" to handle dead pixels
//  Take truncation size from new pixmax information
//  Change to allow template sizes to be changed at compile time
//  Move interpolation range error to LogDebug
//  Add q2bin = 5 and change 1-pixel probability to use new template info
//  Add floor for probabilities (no exact zeros)
//  Add ambiguity resolution with crude 2-D templates (v2.00)
//  Pass all containers by alias to prevent excessive cpu-usage (v2.01)
//  Add ambiguity resolution with fancy 2-D templates (v2.10)
//  Make small change to indices for ambiguity resolution (v2.11)
//  Tune x and y errors separately (v2.12)
//  Use template cytemp/cxtemp methods to center the data cluster in the right place when the templates become asymmetric after irradiation (v2.20)
//  Add charge probability to the splitter [tests consistency with a two-hit merged cluster hypothesis]  (v2.20)
//  Improve likelihood normalization slightly (v2.21)
//  Replace hardwired pixel size derived errors with ones from templated pixel sizes (v2.22)
//  Add shape and charge probabilities for the merged cluster hypothesis (v2.23)
//  Incorporate VI-like speed improvements (v2.25)
//  Improve speed by eliminating the triple index boost::multiarray objects and add speed switch to optimize the algorithm (v2.30)
//

#ifndef SI_PIXEL_TEMPLATE_STANDALONE
//#include <cmath.h>
#else
#include <math.h>
#endif
#include <algorithm>
#include <vector>
#include <list>
#include <iostream>
// ROOT::Math has a c++ function that does the probability calc, but only in v5.12 and later
#include "Math/DistFunc.h"
#include "TMath.h"
// Use current version of gsl instead of ROOT::Math
//#include <gsl/gsl_cdf.h>

#ifndef SI_PIXEL_TEMPLATE_STANDALONE
#include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateSplit.h"
#include "RecoLocalTracker/SiPixelRecHits/interface/VVIObj.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#define LOGERROR(x) edm::LogError(x)
#define LOGDEBUG(x) LogDebug(x)
static int theVerboseLevel = 2;
#define ENDL " "
#else
#include "SiPixelTemplateSplit.h"
#include "VVIObj.h"
//#include "SiPixelTemplate2D.h"
//#include "SimpleTemplate2D.h"
//static int theVerboseLevel = {2};
#define LOGERROR(x) std::cout << x << ": "
#define LOGDEBUG(x) std::cout << x << ": "
#define ENDL std::endl
#endif

using namespace SiPixelTemplateSplit;

// *************************************************************************************************************************************
// *************************************************************************************************************************************
int SiPixelTemplateSplit::PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d& clust, 
                    std::vector<bool>& ydouble, std::vector<bool>& xdouble, 
                    SiPixelTemplate& templ, 
                    float& yrec1, float& yrec2, float& sigmay, float& prob2y,
                        float& xrec1, float& xrec2, float& sigmax, float& prob2x, int& q2bin, float& prob2Q, bool resolve, int speed, float& dchisq, bool deadpix, std::vector<std::pair<int, int> >& zeropix, SiPixelTemplate2D& templ2D)
                        
{
    // Local variables 
        int i, j, k, binq, midpix, fypix, nypix, lypix, logypx, lparm;
        int fxpix, nxpix, lxpix, logxpx, shifty, shiftx, nyzero[TYSIZE];
        int nclusx, nclusy;
        int nybin, ycbin, nxbin, xcbin, minbinj, minbink;
        int deltaj, jmin, jmax, kmin, kmax, km, fxbin, lxbin, fybin, lybin, djy, djx;
        float sythr, sxthr, delta, sigma, sigavg, pseudopix, xsize, ysize, qscale, lanpar[2][5];
        float ss2, ssa, sa2, rat, fq, qtotal, qpixel, qavg;
        float x1p, x2p, y1p, y2p, deltachi2;
        float originx, originy, bias, maxpix, minmax;
        double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
        double hchi2, hndof, sigmal1, sigmal2, mpv1, mpv2, arg1, arg2, q05, like, loglike1, loglike2;
        double prvav, mpv, sigmaQ, kappa, xvav, beta2;
        float ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
        float ysig2[BYSIZE], xsig2[BXSIZE];
        float yw2[BYSIZE], xw2[BXSIZE],  ysw[BYSIZE], xsw[BXSIZE];
        float cluster2[BXM2][BYM2], temp2d1[BXM2][BYM2], temp2d2[BXM2][BYM2];
        bool yd[BYSIZE], xd[BXSIZE], anyyd, anyxd, any2dfail;
        const float sqrt2x={3.0000}, sqrt2y={1.7000};
        const float sqrt12={3.4641};
        const float probmin={1.110223e-16};
        const float prob2Qmin={1.e-5};
        std::pair<int, int> pixel;
        
//      bool SiPixelTemplateSplit::SimpleTemplate2D(float cotalpha, float cotbeta, float xhit, float yhit, float thick, float lorxwidth, float lorywidth, 
//                                                float qavg, std::vector<bool> ydouble, std::vector<bool> xdouble, float template2d[BXM2][BYM2]);
        
// The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only

        const double mean1pix={0.100}, chi21min={0.160};
                      
// First, interpolate the template needed to analyze this cluster     
// check to see of the track direction is in the physical range of the loaded template

        if(!templ.interpolate(id, cotalpha, cotbeta)) {
           LOGDEBUG("SiPixelTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha << ", cot(beta) = " << cotbeta 
           << " is not within the acceptance of template ID = " << id << ", no reconstruction performed" << ENDL;       
           return 20;
        }
                      
// Get pixel dimensions from the template (to allow multiple detectors in the future)
        
        xsize = templ.xsize();
        ysize = templ.ysize();
        
// Define size of pseudopixel

        pseudopix = templ.s50();
//      q05 = 0.28*pseudopix;
        q05 = 0.05f*pseudopix;
        
// Get charge scaling factor

        qscale = templ.qscale();
        
// Make a local copy of the cluster container so that we can muck with it
        
        array_2d cluster = clust;
    
// Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags

        if(cluster.num_dimensions() != 2) {
           LOGERROR("SiPixelTemplateReco") << "input cluster container (BOOST Multiarray) has wrong number of dimensions" << ENDL;      
           return 3;
        }
        nclusx = (int)cluster.shape()[0];
        nclusy = (int)cluster.shape()[1];
        if(nclusx != (int)xdouble.size()) {
           LOGERROR("SiPixelTemplateReco") << "input cluster container x-size is not equal to double pixel flag container size" << ENDL;        
           return 4;
        }
        if(nclusy != (int)ydouble.size()) {
           LOGERROR("SiPixelTemplateReco") << "input cluster container y-size is not equal to double pixel flag container size" << ENDL;        
           return 5;
        }
        
// enforce maximum size 
        
        if(nclusx > TXSIZE) {nclusx = TXSIZE;}
        if(nclusy > TYSIZE) {nclusy = TYSIZE;}
        
// First, rescale all pixel charges       

    for(i=0; i<nclusy; ++i) {
           for(j=0; j<nclusx; ++j) {
                  if(cluster[j][i] > 0) {cluster[j][i] *= qscale;}
           }
        }
        
// Next, sum the total charge and "decapitate" big pixels         

        qtotal = 0.f;
        minmax = 2.0f*templ.pixmax();
        for(i=0; i<nclusy; ++i) {
           maxpix = minmax;
           if(ydouble[i]) {maxpix *=2.f;}
           for(j=0; j<nclusx; ++j) {
                  qtotal += cluster[j][i];
                  if(cluster[j][i] > maxpix) {cluster[j][i] = maxpix;}
           }
        }
        
// Do the cluster repair here   
        
    if(deadpix) {
           fypix = BYM3; lypix = -1;
       for(i=0; i<nclusy; ++i) {
              ysum[i] = 0.f; nyzero[i] = 0;
// Do preliminary cluster projection in y
              for(j=0; j<nclusx; ++j) {
                     ysum[i] += cluster[j][i];
                  }
                  if(ysum[i] > 0.f) {
// identify ends of cluster to determine what the missing charge should be
                     if(i < fypix) {fypix = i;}
                         if(i > lypix) {lypix = i;}
                  }
           }
           
// Now loop over dead pixel list and "fix" everything   

//First see if the cluster ends are redefined and that we have only one dead pixel per column

           std::vector<std::pair<int, int> >::const_iterator zeroIter = zeropix.begin(), zeroEnd = zeropix.end();
       for ( ; zeroIter != zeroEnd; ++zeroIter ) {
              i = zeroIter->second;
                  if(i<0 || i>TYSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column y-index " << i << ", no reconstruction performed" << ENDL;       
                               return 11;}
                                                   
// count the number of dead pixels in each column
                  ++nyzero[i];
// allow them to redefine the cluster ends
                  if(i < fypix) {fypix = i;}
                  if(i > lypix) {lypix = i;}
           }
           
           nypix = lypix-fypix+1;
           
// 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
           
       for (zeroIter = zeropix.begin(); zeroIter != zeroEnd; ++zeroIter ) {        
              i = zeroIter->second; j = zeroIter->first;
                  if(j<0 || j>TXSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column x-index " << j << ", no reconstruction performed" << ENDL;       
                               return 12;}
                  if((i == fypix || i == lypix) && nypix > 1) {maxpix = templ.symax()/2.;} else {maxpix = templ.symax();}
                  if(ydouble[i]) {maxpix *=2.;}
                  if(nyzero[i] > 0 && nyzero[i] < 3) {qpixel = (maxpix - ysum[i])/(float)nyzero[i];} else {qpixel = 1.;}
                  if(qpixel < 1.f) {qpixel = 1.f;}
          cluster[j][i] = qpixel;
// Adjust the total cluster charge to reflect the charge of the "repaired" cluster
                  qtotal += qpixel;
           }
// End of cluster repair section
        } 
        
// Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container         

    for(i=0; i<BYSIZE; ++i) { ysum[i] = 0.f; yd[i] = false;}
        k=0;
        anyyd = false;
    for(i=0; i<nclusy; ++i) {
           for(j=0; j<nclusx; ++j) {
                  ysum[k] += cluster[j][i];
           }
    
// If this is a double pixel, put 1/2 of the charge in 2 consective single pixels  
   
           if(ydouble[i]) {
              ysum[k] /= 2.f;
                  ysum[k+1] = ysum[k];
                  yd[k] = true;
                  yd[k+1] = false;
                  k=k+2;
                  anyyd = true;
           } else {
                  yd[k] = false;
              ++k;
           }
           if(k > BYM1) {break;}
        }
                 
// Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container         

    for(i=0; i<BXSIZE; ++i) { xsum[i] = 0.f; xd[i] = false;}
        k=0;
        anyxd = false;
    for(j=0; j<nclusx; ++j) {
           for(i=0; i<nclusy; ++i) {
                  xsum[k] += cluster[j][i];
           }
    
// If this is a double pixel, put 1/2 of the charge in 2 consective single pixels  
   
           if(xdouble[j]) {
              xsum[k] /= 2.f;
                  xsum[k+1] = xsum[k];
                  xd[k]=true;
                  xd[k+1]=false;
                  k=k+2;
                  anyxd = true;
           } else {
                  xd[k]=false;
              ++k;
           }
           if(k > BXM1) {break;}
        }
        
// next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx   

    fypix=-1;
        nypix=0;
        lypix=0;
        logypx=0;
        for(i=0; i<BYSIZE; ++i) {
           if(ysum[i] > 0.) {
              if(fypix == -1) {fypix = i;}
                  if(!yd[i]) {
                     ysort[logypx] = ysum[i];
                         ++logypx;
                  }
                  ++nypix;
                  lypix = i;
                }
        }
        
// Make sure cluster is continuous

        if((lypix-fypix+1) != nypix || nypix == 0) { 
           LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
           if (theVerboseLevel > 2) {
          LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
          for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}           
                  LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
       }
        
           return 1; 
        }
        
// If cluster is longer than max template size, technique fails

        if(nypix > TYSIZE) { 
           LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster is larger than maximum template size" << ENDL;
           if (theVerboseLevel > 2) {
          LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
          for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}           
                  LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
       }
        
           return 6; 
        }
        
// next, center the cluster on template center if necessary   

        midpix = (fypix+lypix)/2;
//      shifty = BHY - midpix;
        shifty = templ.cytemp() - midpix;
        if(shifty > 0) {
           for(i=lypix; i>=fypix; --i) {
              ysum[i+shifty] = ysum[i];
                  ysum[i] = 0.;
                  yd[i+shifty] = yd[i];
                  yd[i] = false;
           }
        } else if (shifty < 0) {
           for(i=fypix; i<=lypix; ++i) {
              ysum[i+shifty] = ysum[i];
                  ysum[i] = 0.;
                  yd[i+shifty] = yd[i];
                  yd[i] = false;
           }
    }
        lypix +=shifty;
        fypix +=shifty;
        
// If the cluster boundaries are OK, add pesudopixels, otherwise quit
        
        if(fypix > 1 && fypix < BYM2) {
           ysum[fypix-1] = pseudopix;
           ysum[fypix-2] = 0.2f*pseudopix;
        } else {return 8;}
        if(lypix > 1 && lypix < BYM2) {
           ysum[lypix+1] = pseudopix;   
           ysum[lypix+2] = 0.2f*pseudopix;
        } else {return 8;}
        
// finally, determine if pixel[0] is a double pixel and make an origin correction if it is   

    if(ydouble[0]) {
           originy = -0.5f;
        } else {
           originy = 0.f;
        }
        
// next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx   

    fxpix=-1;
        nxpix=0;
        lxpix=0;
        logxpx=0;
        for(i=0; i<BXSIZE; ++i) {
           if(xsum[i] > 0.) {
              if(fxpix == -1) {fxpix = i;}
                  if(!xd[i]) {
                     xsort[logxpx] = xsum[i];
                         ++logxpx;
                  }
                  ++nxpix;
                  lxpix = i;
                }
        }
        
// Make sure cluster is continuous

        if((lxpix-fxpix+1) != nxpix) { 
        
           LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
           if (theVerboseLevel > 2) {
          LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
          for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}           
                  LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
       }

           return 2; 
        }

// If cluster is longer than max template size, technique fails

        if(nxpix > TXSIZE) { 
        
           LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
           if (theVerboseLevel > 2) {
          LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
          for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}           
                  LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
       }

           return 7; 
        }
        
// next, center the cluster on template center if necessary   

        midpix = (fxpix+lxpix)/2;
//      shiftx = BHX - midpix;
        shiftx = templ.cxtemp() - midpix;
        if(shiftx > 0) {
           for(i=lxpix; i>=fxpix; --i) {
              xsum[i+shiftx] = xsum[i];
                  xsum[i] = 0.f;
              xd[i+shiftx] = xd[i];
                  xd[i] = false;
           }
        } else if (shiftx < 0) {
           for(i=fxpix; i<=lxpix; ++i) {
              xsum[i+shiftx] = xsum[i];
                  xsum[i] = 0.f;
              xd[i+shiftx] = xd[i];
                  xd[i] = false;
           }
    }
        lxpix +=shiftx;
        fxpix +=shiftx;
        
// If the cluster boundaries are OK, add pesudopixels, otherwise quit
        
        if(fxpix > 1 && fxpix <BXM2) {
           xsum[fxpix-1] = pseudopix;
           xsum[fxpix-2] = 0.2f*pseudopix;
        } else {return 9;}
        if(lxpix > 1 && lxpix < BXM2) {
           xsum[lxpix+1] = pseudopix;
           xsum[lxpix+2] = 0.2f*pseudopix;
        } else {return 9;}
                        
// finally, determine if pixel[0] is a double pixel and make an origin correction if it is   

    if(xdouble[0]) {
           originx = -0.5f;
        } else {
           originx = 0.f;
        }
        
// uncertainty and final corrections depend upon total charge bin          
           
        qavg = templ.qavg();
        fq = qtotal/qavg;
        if(fq > 3.0f) {
           binq=0;
        } else {
           if(fq > 2.0f) {
              binq=1;
           } else {
                  if(fq > 1.70f) {
                         binq=2;
                  } else {
                         binq=3;
                  }
           }
        }
        
                
// Calculate the Vavilov probability that the cluster charge is consistent with a merged cluster
        
        if(speed < 0) {
           templ.vavilov2_pars(mpv, sigmaQ, kappa);
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
           if((sigmaQ <=0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
                        throw cms::Exception("DataCorrupt") << "SiPixelTemplateSplit::Vavilov parameters mpv/sigmaQ/kappa = " << mpv << "/" << sigmaQ << "/" << kappa << std::endl;
           }
#else
           assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
#endif
           xvav = ((double)qtotal-mpv)/sigmaQ;
           beta2 = 1.;
//  VVIObj is a private port of CERNLIB VVIDIS
           VVIObj vvidist(kappa, beta2, 1);
           prvav = vvidist.fcn(xvav);                   
           prob2Q = 1. - prvav;
           if(prob2Q < prob2Qmin) {prob2Q = prob2Qmin;}
        } else {
                prob2Q = -1.f;
        }
        
        
// Return the charge bin via the parameter list unless the charge is too small (then flag it)
        
        q2bin = binq;
        if(!deadpix && qtotal < 1.9f*templ.qmin()) {q2bin = 5;} else {
                if(!deadpix && qtotal < 1.9f*templ.qmin(1)) {q2bin = 4;}
        }
        
        if (theVerboseLevel > 9) {
       LOGDEBUG("SiPixelTemplateSplit") <<
        "ID = " << id << 
         " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta << 
         " nclusx = " << nclusx << " nclusy = " << nclusy << ENDL;
    }

        
// Next, generate the 3d y- and x-templates
   
        templ.ytemp3d_int(nypix, nybin);
        
        ycbin = nybin/2;        
   
        templ.xtemp3d_int(nxpix, nxbin);
        
// retrieve the number of x-bins 
        
        xcbin = nxbin/2;
                
        
// First, decide on chi^2 min search parameters
        
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
        if(speed < -1 || speed > 2) {
                throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::PixelTempReco2D called with illegal speed = " << speed << std::endl;
        }
#else
        assert(speed >= -1 && speed < 3);
#endif
        fybin = 0; lybin = nybin-1; fxbin = 0; lxbin = nxbin-1; djy = 1; djx = 1;
        if(speed > 0) {
           djy = 2; djx = 2;
           if(speed > 1) {
              if(!anyyd) {djy = 4;}
                        if(!anyxd) {djx = 4;}
           }
        }
        
        if (theVerboseLevel > 9) {
                LOGDEBUG("SiPixelTemplateReco") <<
                "fypix " << fypix << " lypix = " << lypix << 
                " fybin = " << fybin << " lybin = " << lybin << 
                " djy = " << djy << " logypx = " << logypx << ENDL;
                LOGDEBUG("SiPixelTemplateReco") <<
                "fxpix " << fxpix << " lxpix = " << lxpix << 
                " fxbin = " << fxbin << " lxbin = " << lxbin << 
                " djx = " << djx << " logxpx = " << logxpx << ENDL;
        }
        
        
        
// Do the y-reconstruction first 
             
// Define the maximum signal to allow before de-weighting a pixel 

        sythr = 2.1f*(templ.symax());
                          
// Make sure that there will be at least two pixels that are not de-weighted 

        std::sort(&ysort[0], &ysort[logypx]);
        if(logypx == 1) {sythr = 1.01f*ysort[0];} else {
           if (ysort[1] > sythr) { sythr = 1.01f*ysort[1]; }
        }
        
// Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 

//      for(i=0; i<BYSIZE; ++i) { ysig2[i] = 0.;}
        templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
                          
// Find the template bin that minimizes the Chi^2 

        chi2ymin = 1.e15;
        ss2 = 0.f;
        for(i=fypix-2; i<=lypix+2; ++i) { 
                yw2[i] = 1.f/ysig2[i];
                ysw[i] = ysum[i]*yw2[i];
                ss2 += ysum[i]*ysw[i];
        }
        minbinj = -1; 
        minbink = -1;
        deltaj = djy;
        jmin = fybin;
        jmax = lybin;
        kmin = fybin;
        kmax = lybin;
    std::vector<float> ytemp(BYSIZE);
        while(deltaj > 0) {
           for(j=jmin; j<jmax; j+=deltaj) {
                        km = std::min(kmax, j);
                for(k=kmin; k<=km; k+=deltaj) {
                
// Get the template for this set of indices
                
               templ.ytemp3d(j, k, ytemp);
                
// Modify the template if double pixels are present   
                
                if(nypix > logypx) {
                    i=fypix;
                    while(i < lypix) {
                       if(yd[i] && !yd[i+1]) {
                                    
// Sum the adjacent cells and put the average signal in both   
                                    
                           sigavg = (ytemp[i] +  ytemp[i+1])/2.f;
                           ytemp[i] = sigavg;
                           ytemp[i+1] = sigavg;
                           i += 2;
                        } else {
                            ++i;
                        }
                    }
                }       
                      ssa = 0.f;
                      sa2 = 0.f;
                      for(i=fypix-2; i<=lypix+2; ++i) {
                              ssa += ysw[i]*ytemp[i];
                              sa2 += ytemp[i]*ytemp[i]*yw2[i];
                      }
                      rat=ssa/ss2;
                      if(rat <= 0.) {LOGERROR("SiPixelTemplateSplit") << "illegal chi2ymin normalization = " << rat << ENDL; rat = 1.;}
                      chi2y=ss2-2.f*ssa/rat+sa2/(rat*rat);
                      if(chi2y < chi2ymin) {
                              chi2ymin = chi2y;
                              minbinj = j;
                              minbink = k;
                      }
              } 
           }
       deltaj /= 2;
           if(minbinj > fybin) {jmin = minbinj - deltaj;} else {jmin = fybin;}
           if(minbinj < lybin) {jmax = minbinj + deltaj;} else {jmax = lybin;}
           if(minbink > fybin) {kmin = minbink - deltaj;} else {kmin = fybin;}
           if(minbink < lybin) {kmax = minbink + deltaj;} else {kmax = lybin;}          
        }
        
        if (theVerboseLevel > 9) {
       LOGDEBUG("SiPixelTemplateReco") <<
        "minbins " << minbinj << "," << minbink << " chi2ymin = " << chi2ymin << ENDL;
    }
        
// Do not apply final template pass to 1-pixel clusters (use calibrated offset) 
        
        if(logypx == 1) {
        
           if(nypix ==1) {
              delta = templ.dyone();
                  sigma = templ.syone();
           } else {
              delta = templ.dytwo();
                  sigma = templ.sytwo();
           }
           
           yrec1 = 0.5f*(fypix+lypix-2*shifty+2.f*originy)*ysize-delta;
           yrec2 = yrec1;
           
           if(sigma <= 0.) {
              sigmay = ysize/sqrt12;
           } else {
          sigmay = sigma;
           }
           
// Do probability calculation for one-pixel clusters

                chi21max = fmax(chi21min, (double)templ.chi2yminone());
                chi2ymin -=chi21max;
                if(chi2ymin < 0.) {chi2ymin = 0.;}
                //         prob2y = gsl_cdf_chisq_Q(chi2ymin, mean1pix);
                meany = fmax(mean1pix, (double)templ.chi2yavgone());
                hchi2 = chi2ymin/2.; hndof = meany/2.;
                prob2y = 1. - TMath::Gamma(hndof, hchi2);
           
        } else {
           
// For cluster > 1 pix, use chi^2 minimm to recontruct the two y-positions

// at small eta, the templates won't actually work on two pixel y-clusters so just return the pixel centers        

           if(logypx == 2 && fabsf(cotbeta) < 0.25f) {
                   switch(nypix) {
                      case 2:
//  Both pixels are small
                     yrec1 = (fypix-shifty+originy)*ysize;
                         yrec2 = (lypix-shifty+originy)*ysize;
                         sigmay = ysize/sqrt12;
                                 break;
                          case 3:
//  One big pixel and one small pixel
                            if(yd[fypix]) {
                                   yrec1 = (fypix+0.5f-shifty+originy)*ysize;
                                   yrec2 = (lypix-shifty+originy)*ysize;
                                   sigmay = ysize/sqrt12;
                                } else {
                                   yrec1 = (fypix-shifty+originy)*ysize;
                                   yrec2 = (lypix-0.5f-shifty+originy)*ysize;
                                   sigmay = 1.5f*ysize/sqrt12;
                                }
                                break;
                         case 4:
//  Two big pixels
                                yrec1 = (fypix+0.5f-shifty+originy)*ysize;
                                yrec2 = (lypix-0.5f-shifty+originy)*ysize;
                                sigmay = 2.f*ysize/sqrt12;
                            break;
                         default:
//  Something is screwy ...
                    LOGERROR("SiPixelTemplateReco") << "weird problem: logical y-pixels = " << logypx << ", total ysize in normal pixels = " << nypix << ENDL;  
                    return 10;
              }
           } else {
           
// uncertainty and final correction depend upon charge bin       
  
              bias = templ.yavgc2m(binq);
              yrec1 = (0.125f*(minbink-ycbin)+BHY-(float)shifty+originy)*ysize - bias;
              yrec2 = (0.125f*(minbinj-ycbin)+BHY-(float)shifty+originy)*ysize - bias;
              sigmay = sqrt2y*templ.yrmsc2m(binq);
                  
           }
           
// Do goodness of fit test in y  
           
           if(chi2ymin < 0.0) {chi2ymin = 0.0;}
           meany = templ.chi2yavgc2m(binq);
           if(meany < 0.01) {meany = 0.01;}
// gsl function that calculates the chi^2 tail prob for non-integral dof
//         prob2y = gsl_cdf_chisq_Q(chi2y, meany);
//         prob2y = ROOT::Math::chisquared_cdf_c(chi2y, meany);
       hchi2 = chi2ymin/2.; hndof = meany/2.;
           prob2y = 1. - TMath::Gamma(hndof, hchi2);
        }
        
// Do the x-reconstruction next 
                          
                                  
// Define the maximum signal to allow before de-weighting a pixel 

        sxthr = 2.1f*templ.sxmax();
                          
// Make sure that there will be at least two pixels that are not de-weighted 
        std::sort(&xsort[0], &xsort[logxpx]);
        if(logxpx == 1) {sxthr = 1.01f*xsort[0];} else {
           if (xsort[1] > sxthr) { sxthr = 1.01f*xsort[1]; }
        }
           
// Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 

//      for(i=0; i<BYSIZE; ++i) { xsig2[i] = 0.; }
        templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
                          
// Find the template bin that minimizes the Chi^2 

        chi2xmin = 1.e15;
        ss2 = 0.f;
        for(i=fxpix-2; i<=lxpix+2; ++i) {
                xw2[i] = 1.f/xsig2[i];
                xsw[i] = xsum[i]*xw2[i];
                ss2 += xsum[i]*xsw[i];
        }
        minbinj = -1; 
        minbink = -1;
        deltaj = djx;
        jmin = fxbin;
        jmax = lxbin;
        kmin = fxbin;
        kmax = lxbin;
    std::vector<float> xtemp(BXSIZE);
        while(deltaj > 0) {
           for(j=jmin; j<jmax; j+=deltaj) {
                        km = std::min(kmax, j);
              for(k=kmin; k<=km; k+=deltaj) {
              
// Get the template for this set of indices
              
              templ.xtemp3d(j, k, xtemp);
              
// Modify the template if double pixels are present   
              
              if(nxpix > logxpx) {
                  i=fxpix;
                  while(i < lxpix) {
                     if(xd[i] && !xd[i+1]) {
                                  
// Sum the adjacent cells and put the average signal in both   
                                  
                        sigavg = (xtemp[i] +  xtemp[i+1])/2.f;
                        xtemp[i] = sigavg;
                        xtemp[i+1] = sigavg;
                          i += 2;
                      } else {
                          ++i;
                      }
                  }
              } 
                      ssa = 0.f;
                      sa2 = 0.f;
                      for(i=fxpix-2; i<=lxpix+2; ++i) {
                              ssa += xsw[i]*xtemp[i];
                              sa2 += xtemp[i]*xtemp[i]*xw2[i];
                      }
                      rat=ssa/ss2;
                      if(rat <= 0.f) {LOGERROR("SiPixelTemplateSplit") << "illegal chi2xmin normalization = " << rat << ENDL; rat = 1.;}
                      chi2x=ss2-2.f*ssa/rat+sa2/(rat*rat);
                      if(chi2x < chi2xmin) {
                              chi2xmin = chi2x;
                              minbinj = j;
                              minbink = k;
                                }
                   }
           } 
       deltaj /= 2;
           if(minbinj > fxbin) {jmin = minbinj - deltaj;} else {jmin = fxbin;}
           if(minbinj < lxbin) {jmax = minbinj + deltaj;} else {jmax = lxbin;}
           if(minbink > fxbin) {kmin = minbink - deltaj;} else {kmin = fxbin;}
           if(minbink < lxbin) {kmax = minbink + deltaj;} else {kmax = lxbin;}          
        }
        
        if (theVerboseLevel > 9) {
       LOGDEBUG("SiPixelTemplateSplit") <<
        "minbinj/minbink " << minbinj<< "/" << minbink << " chi2xmin = " << chi2xmin << ENDL;
    }

// Do not apply final template pass to 1-pixel clusters (use calibrated offset)
        
        if(logxpx == 1) {
        
           if(nxpix ==1) {
              delta = templ.dxone();
                  sigma = templ.sxone();
           } else {
              delta = templ.dxtwo();
                  sigma = templ.sxtwo();
           }
           xrec1 = 0.5f*(fxpix+lxpix-2*shiftx+2.f*originx)*xsize-delta;
           xrec2 = xrec1;
           if(sigma <= 0.f) {
              sigmax = xsize/sqrt12;
           } else {
          sigmax = sigma;
           }
           
// Do probability calculation for one-pixel clusters

                chi21max = fmax(chi21min, (double)templ.chi2xminone());
                chi2xmin -=chi21max;
                if(chi2xmin < 0.) {chi2xmin = 0.;}
                meanx = fmax(mean1pix, (double)templ.chi2xavgone());
                hchi2 = chi2xmin/2.; hndof = meanx/2.;
                prob2x = 1. - TMath::Gamma(hndof, hchi2);
           
        } else {
           
// For cluster > 1 pix, use chi^2 minimm to recontruct the two x-positions

// uncertainty and final correction depend upon charge bin         
           
           bias = templ.xavgc2m(binq);
           k = std::min(minbink, minbinj);
           j = std::max(minbink, minbinj);
       xrec1 = (0.125f*(minbink-xcbin)+BHX-(float)shiftx+originx)*xsize - bias;
       xrec2 = (0.125f*(minbinj-xcbin)+BHX-(float)shiftx+originx)*xsize - bias;
           sigmax = sqrt2x*templ.xrmsc2m(binq);
           
// Do goodness of fit test in y  
           
           if(chi2xmin < 0.0) {chi2xmin = 0.0;}
           meanx = templ.chi2xavgc2m(binq);
           if(meanx < 0.01) {meanx = 0.01;}
       hchi2 = chi2xmin/2.; hndof = meanx/2.;
           prob2x = 1. - TMath::Gamma(hndof, hchi2);
        }
        
        //  Don't return exact zeros for the probability
        
        if(prob2y < probmin) {prob2y = probmin;}
        if(prob2x < probmin) {prob2x = probmin;}
        
//  New code: resolve the ambiguity if resolve is set to true
        
        dchisq = 0.;
        if(resolve) {
                
//  First copy the unexpanded cluster into a new working buffer
                
                for(i=0; i<BYM2; ++i) {
                        for(j=0; j<BXM2; ++j) {
                                if((i>0 && i<BYM3)&&(j>0 && j<BXM3)) {
                                        cluster2[j][i] = qscale*clust[j-1][i-1];
                                } else {
                                        cluster2[j][i] = 0.f;
                                }
                        }
                }
                
//  Next, redefine the local coordinates to start at the lower LH corner of pixel[1][1];
                
                if(xdouble[0]) {
                        x1p = xrec1+xsize;
                        x2p = xrec2+xsize;
                } else {
                        x1p = xrec1+xsize/2.f;
                        x2p = xrec2+xsize/2.f;  
                }
                
                if(ydouble[0]) {
                        y1p = yrec1+ysize;
                        y2p = yrec2+ysize;
                } else {
                        y1p = yrec1+ysize/2.f;
                        y2p = yrec2+ysize/2.f;  
                }
                
// Next, calculate 2-d templates for the (x1,y1)+(x2,y2) and the (x1,y2)+(x2,y1) hypotheses
                
//  First zero the two 2-d hypotheses
                
                for(i=0; i<BYM2; ++i) {
                        for(j=0; j<BXM2; ++j) {
                                temp2d1[j][i] = 0.f;
                                temp2d2[j][i] = 0.f;
                        }
                }
                
                
// Add the two hits in the first hypothesis
                
                any2dfail = templ2D.xytemp(id, cotalpha, cotbeta, x1p, y1p, xdouble, ydouble, temp2d1) && templ2D.xytemp(id, cotalpha, cotbeta, x2p, y2p, xdouble, ydouble, temp2d1);
                
// And then the second hypothesis
                
                any2dfail = any2dfail && templ2D.xytemp(id, cotalpha, cotbeta, x1p, y2p, xdouble, ydouble, temp2d2);
                
                any2dfail = any2dfail && templ2D.xytemp(id, cotalpha, cotbeta, x2p, y1p, xdouble, ydouble, temp2d2);
                
// If any of these have failed, use the simple templates instead
                
                if(!any2dfail) {
                        
//  Rezero the two 2-d hypotheses
                        
                        for(i=0; i<BYM2; ++i) {
                                for(j=0; j<BXM2; ++j) {
                                        temp2d1[j][i] = 0.f;
                                        temp2d2[j][i] = 0.f;
                                }
                        }
                        
// Add the two hits in the first hypothesis
                        
                        if(!templ.simpletemplate2D(x1p, y1p, xdouble, ydouble, temp2d1)) {return 1;}
                        
                        if(!templ.simpletemplate2D(x2p, y2p, xdouble, ydouble, temp2d1)) {return 1;}
                        
// And then the second hypothesis
                        
                        if(!templ.simpletemplate2D(x1p, y2p, xdouble, ydouble, temp2d2)) {return 1;}
                        
                        if(!templ.simpletemplate2D(x2p, y1p, xdouble, ydouble, temp2d2)) {return 1;}
                }
                
// Keep lists of pixels and nearest neighbors   
                
                std::list<std::pair<int, int> > pixellst;
                
// Loop through the array and find non-zero elements
                
                for(i=0; i<BYM2; ++i) {
                        for(j=0; j<BXM2; ++j) {
                                if(cluster2[j][i] > 0.f || temp2d1[j][i] > 0.f || temp2d2[j][i] > 0.f) {
                                        pixel.first = j; pixel.second = i;
                                        pixellst.push_back(pixel);
                                }
                        }
                }
                
                // Now calculate the product of Landau probabilities (alpha probability)
                
      templ2D.landau_par(lanpar);               
                loglike1 = 0.; loglike2 = 0.;
                
                // Now, for each neighbor, match it again the pixel list.
                // If found, delete it from the neighbor list
                
                std::list<std::pair<int, int> >::const_iterator pixIter, pixEnd;
                pixIter = pixellst.begin();
                pixEnd = pixellst.end();
                for( ; pixIter != pixEnd; ++pixIter ) {
                        j = pixIter->first; i = pixIter->second;
                        if((i < BHY && cotbeta > 0.) || (i >= BHY && cotbeta < 0.)) {lparm = 0;} else {lparm = 1;}
                        mpv1 = lanpar[lparm][0] + lanpar[lparm][1]*temp2d1[j][i];
                        sigmal1 = lanpar[lparm][2]*mpv1;
                        sigmal1 = sqrt(sigmal1*sigmal1+lanpar[lparm][3]*lanpar[lparm][3]);
                        if(sigmal1 < q05) sigmal1 = q05;
                        arg1 = (cluster2[j][i]-mpv1)/sigmal1-0.22278;
                        mpv2 = lanpar[lparm][0] + lanpar[lparm][1]*temp2d2[j][i];
                        sigmal2 = lanpar[lparm][2]*mpv2;
                        sigmal2 = sqrt(sigmal2*sigmal2+lanpar[lparm][3]*lanpar[lparm][3]);
                        if(sigmal2 < q05) sigmal2 = q05;
                        arg2 = (cluster2[j][i]-mpv2)/sigmal2-0.22278;
//                      like = ROOT::Math::landau_pdf(arg1)/sigmal1;
                        like = ROOT::Math::landau_pdf(arg1);
                        if(like < 1.e-30) like = 1.e-30;
                        loglike1 += log(like);
//                      like = ROOT::Math::landau_pdf(arg2)/sigmal2;
                        like = ROOT::Math::landau_pdf(arg2);
                        if(like < 1.e-30) like = 1.e-30;
                        loglike2 += log(like);
                }
                
// Calculate chisquare difference for the two hypotheses 9don't multiply by 2 for less inconvenient scaling
                
                deltachi2 = loglike1 - loglike2;
                                
                if(deltachi2 < 0.) {
                        
                        // Flip the x1 and x2
                        
                        x1p = xrec1;
                        xrec1 = xrec2;
                        xrec2 = x1p;
                }
                
                // Return a positive definite value
                
                dchisq = fabs(deltachi2);               
        }
        
    return 0;
} // PixelTempSplit 


// *************************************************************************************************************************************
// *************************************************************************************************************************************

int SiPixelTemplateSplit::PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d& cluster, 
                    std::vector<bool>& ydouble, std::vector<bool>& xdouble, 
                    SiPixelTemplate& templ, 
                    float& yrec1, float& yrec2, float& sigmay, float& prob2y,
                        float& xrec1, float& xrec2, float& sigmax, float& prob2x, int& q2bin, float& prob2Q, bool resolve, int speed, float& dchisq, SiPixelTemplate2D& templ2D)
{
    // Local variables 
        const bool deadpix = false;
        std::vector<std::pair<int, int> > zeropix;
    
        return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ, 
                    yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);

} // PixelTempSplit


// *************************************************************************************************************************************
// *************************************************************************************************************************************

int SiPixelTemplateSplit::PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d& cluster, 
                                         std::vector<bool>& ydouble, std::vector<bool>& xdouble, 
                                         SiPixelTemplate& templ, 
                                         float& yrec1, float& yrec2, float& sigmay, float& prob2y,
                                         float& xrec1, float& xrec2, float& sigmax, float& prob2x, int& q2bin, float& prob2Q, bool resolve, float& dchisq, SiPixelTemplate2D& templ2D)
{
    // Local variables 
        const bool deadpix = false;
        std::vector<std::pair<int, int> > zeropix;
    const int speed = 1;
    
        return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ, 
                                                yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);
    
} // PixelTempSplit

// *************************************************************************************************************************************
// *************************************************************************************************************************************

int SiPixelTemplateSplit::PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d& cluster, 
                                                                                std::vector<bool>& ydouble, std::vector<bool>& xdouble, 
                                                                                SiPixelTemplate& templ, 
                                                                                float& yrec1, float& yrec2, float& sigmay, float& prob2y,
                                                                                float& xrec1, float& xrec2, float& sigmax, float& prob2x, int& q2bin, float& prob2Q, SiPixelTemplate2D& templ2D)
{
    // Local variables 
        const bool deadpix = false;
        const bool resolve = true;
        float dchisq;
        std::vector<std::pair<int, int> > zeropix;
        const int speed = 1;
    
        return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ, 
                                                                                           yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);
        
} // PixelTempSplit



// *************************************************************************************************************************************
// *************************************************************************************************************************************

int SiPixelTemplateSplit::PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d& cluster, 
                                                                                                         std::vector<bool>& ydouble, std::vector<bool>& xdouble, 
                                                                                                         SiPixelTemplate& templ, 
                                                                                                         float& yrec1, float& yrec2, float& sigmay, float& prob2y,
                                                                                                         float& xrec1, float& xrec2, float& sigmax, float& prob2x, int& q2bin, SiPixelTemplate2D& templ2D)
{
        // Local variables 
        const bool deadpix = false;
        const bool resolve = true;
        float dchisq, prob2Q;
        std::vector<std::pair<int, int> > zeropix;
        const int speed = 1;
        
        return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ, 
                                                                                                                         yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);
        
} // PixelTempSplit