Reference Manual

Go to the documentation of this file.
00001 //
00002 //  SiPixelTemplateSplit.cc (Version 2.30)
00003 //
00004 //  Procedure to fit two templates (same angle hypotheses) to a single cluster
00005 //  Return two x- and two y-coordinates for the cluster
00006 //
00007 //  Created by Morris Swartz on 04/10/08.
00008 //
00009 //  Incorporate "cluster repair" to handle dead pixels
00010 //  Take truncation size from new pixmax information
00011 //  Change to allow template sizes to be changed at compile time
00012 //  Move interpolation range error to LogDebug
00013 //  Add q2bin = 5 and change 1-pixel probability to use new template info
00014 //  Add floor for probabilities (no exact zeros)
00015 //  Add ambiguity resolution with crude 2-D templates (v2.00)
00016 //  Pass all containers by alias to prevent excessive cpu-usage (v2.01)
00017 //  Add ambiguity resolution with fancy 2-D templates (v2.10)
00018 //  Make small change to indices for ambiguity resolution (v2.11)
00019 //  Tune x and y errors separately (v2.12)
00020 //  Use template cytemp/cxtemp methods to center the data cluster in the right place when the templates become asymmetric after irradiation (v2.20)
00021 //  Add charge probability to the splitter [tests consistency with a two-hit merged cluster hypothesis]  (v2.20)
00022 //  Improve likelihood normalization slightly (v2.21)
00023 //  Replace hardwired pixel size derived errors with ones from templated pixel sizes (v2.22)
00024 //  Add shape and charge probabilities for the merged cluster hypothesis (v2.23)
00025 //  Incorporate VI-like speed improvements (v2.25)
00026 //  Improve speed by eliminating the triple index boost::multiarray objects and add speed switch to optimize the algorithm (v2.30)
00027 //
00028 
00029 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00030 //#include <cmath.h>
00031 #else
00032 #include <math.h>
00033 #endif
00034 #include <algorithm>
00035 #include <vector>
00036 #include <list>
00037 #include <iostream>
00038 // ROOT::Math has a c++ function that does the probability calc, but only in v5.12 and later
00039 #include "Math/DistFunc.h"
00040 #include "TMath.h"
00041 // Use current version of gsl instead of ROOT::Math
00042 //#include <gsl/gsl_cdf.h>
00043 
00044 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00045 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateSplit.h"
00046 #include "RecoLocalTracker/SiPixelRecHits/interface/VVIObj.h"
00047 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00048 #define LOGERROR(x) edm::LogError(x)
00049 #define LOGDEBUG(x) LogDebug(x)
00050 static int theVerboseLevel = 2;
00051 #define ENDL " "
00052 #else
00053 #include "SiPixelTemplateSplit.h"
00054 #include "VVIObj.h"
00055 //#include "SiPixelTemplate2D.h"
00056 //#include "SimpleTemplate2D.h"
00057 //static int theVerboseLevel = {2};
00058 #define LOGERROR(x) std::cout << x << ": "
00059 #define LOGDEBUG(x) std::cout << x << ": "
00060 #define ENDL std::endl
00061 #endif
00062 
00063 using namespace SiPixelTemplateSplit;
00064 
00065 // *************************************************************************************************************************************
00096 // *************************************************************************************************************************************
00097 int SiPixelTemplateSplit::PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d& clust, 
00098                     std::vector<bool>& ydouble, std::vector<bool>& xdouble, 
00099                     SiPixelTemplate& templ, 
00100                     float& yrec1, float& yrec2, float& sigmay, float& prob2y,
00101                         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)
00102                         
00103 {
00104     // Local variables 
00105         int i, j, k, binq, midpix, fypix, nypix, lypix, logypx, lparm;
00106         int fxpix, nxpix, lxpix, logxpx, shifty, shiftx, nyzero[TYSIZE];
00107         int nclusx, nclusy;
00108         int nybin, ycbin, nxbin, xcbin, minbinj, minbink;
00109         int deltaj, jmin, jmax, kmin, kmax, km, fxbin, lxbin, fybin, lybin, djy, djx;
00110         float sythr, sxthr, delta, sigma, sigavg, pseudopix, xsize, ysize, qscale, lanpar[2][5];
00111         float ss2, ssa, sa2, rat, fq, qtotal, qpixel, qavg;
00112         float x1p, x2p, y1p, y2p, deltachi2;
00113         float originx, originy, bias, maxpix, minmax;
00114         double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
00115         double hchi2, hndof, sigmal1, sigmal2, mpv1, mpv2, arg1, arg2, q05, like, loglike1, loglike2;
00116         double prvav, mpv, sigmaQ, kappa, xvav, beta2;
00117         float ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
00118         float ysig2[BYSIZE], xsig2[BXSIZE];
00119         float yw2[BYSIZE], xw2[BXSIZE],  ysw[BYSIZE], xsw[BXSIZE];
00120         float cluster2[BXM2][BYM2], temp2d1[BXM2][BYM2], temp2d2[BXM2][BYM2];
00121         bool yd[BYSIZE], xd[BXSIZE], anyyd, anyxd, any2dfail;
00122         const float sqrt2x={3.0000}, sqrt2y={1.7000};
00123         const float sqrt12={3.4641};
00124         const float probmin={1.110223e-16};
00125         const float prob2Qmin={1.e-5};
00126         std::pair<int, int> pixel;
00127         
00128 //      bool SiPixelTemplateSplit::SimpleTemplate2D(float cotalpha, float cotbeta, float xhit, float yhit, float thick, float lorxwidth, float lorywidth, 
00129 //                                                float qavg, std::vector<bool> ydouble, std::vector<bool> xdouble, float template2d[BXM2][BYM2]);
00130         
00131 // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
00132 
00133         const double mean1pix={0.100}, chi21min={0.160};
00134                       
00135 // First, interpolate the template needed to analyze this cluster     
00136 // check to see of the track direction is in the physical range of the loaded template
00137 
00138         if(!templ.interpolate(id, cotalpha, cotbeta)) {
00139            LOGDEBUG("SiPixelTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha << ", cot(beta) = " << cotbeta 
00140            << " is not within the acceptance of template ID = " << id << ", no reconstruction performed" << ENDL;       
00141            return 20;
00142         }
00143                       
00144 // Get pixel dimensions from the template (to allow multiple detectors in the future)
00145         
00146         xsize = templ.xsize();
00147         ysize = templ.ysize();
00148         
00149 // Define size of pseudopixel
00150 
00151         pseudopix = templ.s50();
00152 //      q05 = 0.28*pseudopix;
00153         q05 = 0.05f*pseudopix;
00154         
00155 // Get charge scaling factor
00156 
00157         qscale = templ.qscale();
00158         
00159 // Make a local copy of the cluster container so that we can muck with it
00160         
00161         array_2d cluster = clust;
00162     
00163 // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
00164 
00165         if(cluster.num_dimensions() != 2) {
00166            LOGERROR("SiPixelTemplateReco") << "input cluster container (BOOST Multiarray) has wrong number of dimensions" << ENDL;      
00167            return 3;
00168         }
00169         nclusx = (int)cluster.shape()[0];
00170         nclusy = (int)cluster.shape()[1];
00171         if(nclusx != (int)xdouble.size()) {
00172            LOGERROR("SiPixelTemplateReco") << "input cluster container x-size is not equal to double pixel flag container size" << ENDL;        
00173            return 4;
00174         }
00175         if(nclusy != (int)ydouble.size()) {
00176            LOGERROR("SiPixelTemplateReco") << "input cluster container y-size is not equal to double pixel flag container size" << ENDL;        
00177            return 5;
00178         }
00179         
00180 // enforce maximum size 
00181         
00182         if(nclusx > TXSIZE) {nclusx = TXSIZE;}
00183         if(nclusy > TYSIZE) {nclusy = TYSIZE;}
00184         
00185 // First, rescale all pixel charges       
00186 
00187     for(i=0; i<nclusy; ++i) {
00188            for(j=0; j<nclusx; ++j) {
00189                   if(cluster[j][i] > 0) {cluster[j][i] *= qscale;}
00190            }
00191         }
00192         
00193 // Next, sum the total charge and "decapitate" big pixels         
00194 
00195         qtotal = 0.f;
00196         minmax = 2.0f*templ.pixmax();
00197         for(i=0; i<nclusy; ++i) {
00198            maxpix = minmax;
00199            if(ydouble[i]) {maxpix *=2.f;}
00200            for(j=0; j<nclusx; ++j) {
00201                   qtotal += cluster[j][i];
00202                   if(cluster[j][i] > maxpix) {cluster[j][i] = maxpix;}
00203            }
00204         }
00205         
00206 // Do the cluster repair here   
00207         
00208     if(deadpix) {
00209            fypix = BYM3; lypix = -1;
00210        for(i=0; i<nclusy; ++i) {
00211               ysum[i] = 0.f; nyzero[i] = 0;
00212 // Do preliminary cluster projection in y
00213               for(j=0; j<nclusx; ++j) {
00214                      ysum[i] += cluster[j][i];
00215                   }
00216                   if(ysum[i] > 0.f) {
00217 // identify ends of cluster to determine what the missing charge should be
00218                      if(i < fypix) {fypix = i;}
00219                          if(i > lypix) {lypix = i;}
00220                   }
00221            }
00222            
00223 // Now loop over dead pixel list and "fix" everything   
00224 
00225 //First see if the cluster ends are redefined and that we have only one dead pixel per column
00226 
00227            std::vector<std::pair<int, int> >::const_iterator zeroIter = zeropix.begin(), zeroEnd = zeropix.end();
00228        for ( ; zeroIter != zeroEnd; ++zeroIter ) {
00229               i = zeroIter->second;
00230                   if(i<0 || i>TYSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column y-index " << i << ", no reconstruction performed" << ENDL;       
00231                                return 11;}
00232                                                    
00233 // count the number of dead pixels in each column
00234                   ++nyzero[i];
00235 // allow them to redefine the cluster ends
00236                   if(i < fypix) {fypix = i;}
00237                   if(i > lypix) {lypix = i;}
00238            }
00239            
00240            nypix = lypix-fypix+1;
00241            
00242 // 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
00243            
00244        for (zeroIter = zeropix.begin(); zeroIter != zeroEnd; ++zeroIter ) {        
00245               i = zeroIter->second; j = zeroIter->first;
00246                   if(j<0 || j>TXSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column x-index " << j << ", no reconstruction performed" << ENDL;       
00247                                return 12;}
00248                   if((i == fypix || i == lypix) && nypix > 1) {maxpix = templ.symax()/2.;} else {maxpix = templ.symax();}
00249                   if(ydouble[i]) {maxpix *=2.;}
00250                   if(nyzero[i] > 0 && nyzero[i] < 3) {qpixel = (maxpix - ysum[i])/(float)nyzero[i];} else {qpixel = 1.;}
00251                   if(qpixel < 1.f) {qpixel = 1.f;}
00252           cluster[j][i] = qpixel;
00253 // Adjust the total cluster charge to reflect the charge of the "repaired" cluster
00254                   qtotal += qpixel;
00255            }
00256 // End of cluster repair section
00257         } 
00258         
00259 // Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container         
00260 
00261     for(i=0; i<BYSIZE; ++i) { ysum[i] = 0.f; yd[i] = false;}
00262         k=0;
00263         anyyd = false;
00264     for(i=0; i<nclusy; ++i) {
00265            for(j=0; j<nclusx; ++j) {
00266                   ysum[k] += cluster[j][i];
00267            }
00268     
00269 // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels  
00270    
00271            if(ydouble[i]) {
00272               ysum[k] /= 2.f;
00273                   ysum[k+1] = ysum[k];
00274                   yd[k] = true;
00275                   yd[k+1] = false;
00276                   k=k+2;
00277                   anyyd = true;
00278            } else {
00279                   yd[k] = false;
00280               ++k;
00281            }
00282            if(k > BYM1) {break;}
00283         }
00284                  
00285 // Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container         
00286 
00287     for(i=0; i<BXSIZE; ++i) { xsum[i] = 0.f; xd[i] = false;}
00288         k=0;
00289         anyxd = false;
00290     for(j=0; j<nclusx; ++j) {
00291            for(i=0; i<nclusy; ++i) {
00292                   xsum[k] += cluster[j][i];
00293            }
00294     
00295 // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels  
00296    
00297            if(xdouble[j]) {
00298               xsum[k] /= 2.f;
00299                   xsum[k+1] = xsum[k];
00300                   xd[k]=true;
00301                   xd[k+1]=false;
00302                   k=k+2;
00303                   anyxd = true;
00304            } else {
00305                   xd[k]=false;
00306               ++k;
00307            }
00308            if(k > BXM1) {break;}
00309         }
00310         
00311 // next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx   
00312 
00313     fypix=-1;
00314         nypix=0;
00315         lypix=0;
00316         logypx=0;
00317         for(i=0; i<BYSIZE; ++i) {
00318            if(ysum[i] > 0.) {
00319               if(fypix == -1) {fypix = i;}
00320                   if(!yd[i]) {
00321                      ysort[logypx] = ysum[i];
00322                          ++logypx;
00323                   }
00324                   ++nypix;
00325                   lypix = i;
00326                 }
00327         }
00328         
00329 // Make sure cluster is continuous
00330 
00331         if((lypix-fypix+1) != nypix || nypix == 0) { 
00332            LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
00333            if (theVerboseLevel > 2) {
00334           LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
00335           for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}           
00336                   LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
00337        }
00338         
00339            return 1; 
00340         }
00341         
00342 // If cluster is longer than max template size, technique fails
00343 
00344         if(nypix > TYSIZE) { 
00345            LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster is larger than maximum template size" << ENDL;
00346            if (theVerboseLevel > 2) {
00347           LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
00348           for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}           
00349                   LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
00350        }
00351         
00352            return 6; 
00353         }
00354         
00355 // next, center the cluster on template center if necessary   
00356 
00357         midpix = (fypix+lypix)/2;
00358 //      shifty = BHY - midpix;
00359         shifty = templ.cytemp() - midpix;
00360         if(shifty > 0) {
00361            for(i=lypix; i>=fypix; --i) {
00362               ysum[i+shifty] = ysum[i];
00363                   ysum[i] = 0.;
00364                   yd[i+shifty] = yd[i];
00365                   yd[i] = false;
00366            }
00367         } else if (shifty < 0) {
00368            for(i=fypix; i<=lypix; ++i) {
00369               ysum[i+shifty] = ysum[i];
00370                   ysum[i] = 0.;
00371                   yd[i+shifty] = yd[i];
00372                   yd[i] = false;
00373            }
00374     }
00375         lypix +=shifty;
00376         fypix +=shifty;
00377         
00378 // If the cluster boundaries are OK, add pesudopixels, otherwise quit
00379         
00380         if(fypix > 1 && fypix < BYM2) {
00381            ysum[fypix-1] = pseudopix;
00382            ysum[fypix-2] = 0.2f*pseudopix;
00383         } else {return 8;}
00384         if(lypix > 1 && lypix < BYM2) {
00385            ysum[lypix+1] = pseudopix;   
00386            ysum[lypix+2] = 0.2f*pseudopix;
00387         } else {return 8;}
00388         
00389 // finally, determine if pixel[0] is a double pixel and make an origin correction if it is   
00390 
00391     if(ydouble[0]) {
00392            originy = -0.5f;
00393         } else {
00394            originy = 0.f;
00395         }
00396         
00397 // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx   
00398 
00399     fxpix=-1;
00400         nxpix=0;
00401         lxpix=0;
00402         logxpx=0;
00403         for(i=0; i<BXSIZE; ++i) {
00404            if(xsum[i] > 0.) {
00405               if(fxpix == -1) {fxpix = i;}
00406                   if(!xd[i]) {
00407                      xsort[logxpx] = xsum[i];
00408                          ++logxpx;
00409                   }
00410                   ++nxpix;
00411                   lxpix = i;
00412                 }
00413         }
00414         
00415 // Make sure cluster is continuous
00416 
00417         if((lxpix-fxpix+1) != nxpix) { 
00418         
00419            LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
00420            if (theVerboseLevel > 2) {
00421           LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
00422           for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}           
00423                   LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
00424        }
00425 
00426            return 2; 
00427         }
00428 
00429 // If cluster is longer than max template size, technique fails
00430 
00431         if(nxpix > TXSIZE) { 
00432         
00433            LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
00434            if (theVerboseLevel > 2) {
00435           LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
00436           for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}           
00437                   LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
00438        }
00439 
00440            return 7; 
00441         }
00442         
00443 // next, center the cluster on template center if necessary   
00444 
00445         midpix = (fxpix+lxpix)/2;
00446 //      shiftx = BHX - midpix;
00447         shiftx = templ.cxtemp() - midpix;
00448         if(shiftx > 0) {
00449            for(i=lxpix; i>=fxpix; --i) {
00450               xsum[i+shiftx] = xsum[i];
00451                   xsum[i] = 0.f;
00452               xd[i+shiftx] = xd[i];
00453                   xd[i] = false;
00454            }
00455         } else if (shiftx < 0) {
00456            for(i=fxpix; i<=lxpix; ++i) {
00457               xsum[i+shiftx] = xsum[i];
00458                   xsum[i] = 0.f;
00459               xd[i+shiftx] = xd[i];
00460                   xd[i] = false;
00461            }
00462     }
00463         lxpix +=shiftx;
00464         fxpix +=shiftx;
00465         
00466 // If the cluster boundaries are OK, add pesudopixels, otherwise quit
00467         
00468         if(fxpix > 1 && fxpix <BXM2) {
00469            xsum[fxpix-1] = pseudopix;
00470            xsum[fxpix-2] = 0.2f*pseudopix;
00471         } else {return 9;}
00472         if(lxpix > 1 && lxpix < BXM2) {
00473            xsum[lxpix+1] = pseudopix;
00474            xsum[lxpix+2] = 0.2f*pseudopix;
00475         } else {return 9;}
00476                         
00477 // finally, determine if pixel[0] is a double pixel and make an origin correction if it is   
00478 
00479     if(xdouble[0]) {
00480            originx = -0.5f;
00481         } else {
00482            originx = 0.f;
00483         }
00484         
00485 // uncertainty and final corrections depend upon total charge bin          
00486            
00487         qavg = templ.qavg();
00488         fq = qtotal/qavg;
00489         if(fq > 3.0f) {
00490            binq=0;
00491         } else {
00492            if(fq > 2.0f) {
00493               binq=1;
00494            } else {
00495                   if(fq > 1.70f) {
00496                          binq=2;
00497                   } else {
00498                          binq=3;
00499                   }
00500            }
00501         }
00502         
00503                 
00504 // Calculate the Vavilov probability that the cluster charge is consistent with a merged cluster
00505         
00506         if(speed < 0) {
00507            templ.vavilov2_pars(mpv, sigmaQ, kappa);
00508 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00509            if((sigmaQ <=0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
00510                         throw cms::Exception("DataCorrupt") << "SiPixelTemplateSplit::Vavilov parameters mpv/sigmaQ/kappa = " << mpv << "/" << sigmaQ << "/" << kappa << std::endl;
00511            }
00512 #else
00513            assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
00514 #endif
00515            xvav = ((double)qtotal-mpv)/sigmaQ;
00516            beta2 = 1.;
00517 //  VVIObj is a private port of CERNLIB VVIDIS
00518            VVIObj vvidist(kappa, beta2, 1);
00519            prvav = vvidist.fcn(xvav);                   
00520            prob2Q = 1. - prvav;
00521            if(prob2Q < prob2Qmin) {prob2Q = prob2Qmin;}
00522         } else {
00523                 prob2Q = -1.f;
00524         }
00525         
00526         
00527 // Return the charge bin via the parameter list unless the charge is too small (then flag it)
00528         
00529         q2bin = binq;
00530         if(!deadpix && qtotal < 1.9f*templ.qmin()) {q2bin = 5;} else {
00531                 if(!deadpix && qtotal < 1.9f*templ.qmin(1)) {q2bin = 4;}
00532         }
00533         
00534         if (theVerboseLevel > 9) {
00535        LOGDEBUG("SiPixelTemplateSplit") <<
00536         "ID = " << id << 
00537          " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta << 
00538          " nclusx = " << nclusx << " nclusy = " << nclusy << ENDL;
00539     }
00540 
00541         
00542 // Next, generate the 3d y- and x-templates
00543    
00544         templ.ytemp3d_int(nypix, nybin);
00545         
00546         ycbin = nybin/2;        
00547    
00548         templ.xtemp3d_int(nxpix, nxbin);
00549         
00550 // retrieve the number of x-bins 
00551         
00552         xcbin = nxbin/2;
00553                 
00554         
00555 // First, decide on chi^2 min search parameters
00556         
00557 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00558         if(speed < -1 || speed > 2) {
00559                 throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::PixelTempReco2D called with illegal speed = " << speed << std::endl;
00560         }
00561 #else
00562         assert(speed >= -1 && speed < 3);
00563 #endif
00564         fybin = 0; lybin = nybin-1; fxbin = 0; lxbin = nxbin-1; djy = 1; djx = 1;
00565         if(speed > 0) {
00566            djy = 2; djx = 2;
00567            if(speed > 1) {
00568               if(!anyyd) {djy = 4;}
00569                         if(!anyxd) {djx = 4;}
00570            }
00571         }
00572         
00573         if (theVerboseLevel > 9) {
00574                 LOGDEBUG("SiPixelTemplateReco") <<
00575                 "fypix " << fypix << " lypix = " << lypix << 
00576                 " fybin = " << fybin << " lybin = " << lybin << 
00577                 " djy = " << djy << " logypx = " << logypx << ENDL;
00578                 LOGDEBUG("SiPixelTemplateReco") <<
00579                 "fxpix " << fxpix << " lxpix = " << lxpix << 
00580                 " fxbin = " << fxbin << " lxbin = " << lxbin << 
00581                 " djx = " << djx << " logxpx = " << logxpx << ENDL;
00582         }
00583         
00584         
00585         
00586 // Do the y-reconstruction first 
00587              
00588 // Define the maximum signal to allow before de-weighting a pixel 
00589 
00590         sythr = 2.1f*(templ.symax());
00591                           
00592 // Make sure that there will be at least two pixels that are not de-weighted 
00593 
00594         std::sort(&ysort[0], &ysort[logypx]);
00595         if(logypx == 1) {sythr = 1.01f*ysort[0];} else {
00596            if (ysort[1] > sythr) { sythr = 1.01f*ysort[1]; }
00597         }
00598         
00599 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
00600 
00601 //      for(i=0; i<BYSIZE; ++i) { ysig2[i] = 0.;}
00602         templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
00603                           
00604 // Find the template bin that minimizes the Chi^2 
00605 
00606         chi2ymin = 1.e15;
00607         ss2 = 0.f;
00608         for(i=fypix-2; i<=lypix+2; ++i) { 
00609                 yw2[i] = 1.f/ysig2[i];
00610                 ysw[i] = ysum[i]*yw2[i];
00611                 ss2 += ysum[i]*ysw[i];
00612         }
00613         minbinj = -1; 
00614         minbink = -1;
00615         deltaj = djy;
00616         jmin = fybin;
00617         jmax = lybin;
00618         kmin = fybin;
00619         kmax = lybin;
00620     std::vector<float> ytemp(BYSIZE);
00621         while(deltaj > 0) {
00622            for(j=jmin; j<jmax; j+=deltaj) {
00623                         km = std::min(kmax, j);
00624                 for(k=kmin; k<=km; k+=deltaj) {
00625                 
00626 // Get the template for this set of indices
00627                 
00628                templ.ytemp3d(j, k, ytemp);
00629                 
00630 // Modify the template if double pixels are present   
00631                 
00632                 if(nypix > logypx) {
00633                     i=fypix;
00634                     while(i < lypix) {
00635                        if(yd[i] && !yd[i+1]) {
00636                                     
00637 // Sum the adjacent cells and put the average signal in both   
00638                                     
00639                            sigavg = (ytemp[i] +  ytemp[i+1])/2.f;
00640                            ytemp[i] = sigavg;
00641                            ytemp[i+1] = sigavg;
00642                            i += 2;
00643                         } else {
00644                             ++i;
00645                         }
00646                     }
00647                 }       
00648                       ssa = 0.f;
00649                       sa2 = 0.f;
00650                       for(i=fypix-2; i<=lypix+2; ++i) {
00651                               ssa += ysw[i]*ytemp[i];
00652                               sa2 += ytemp[i]*ytemp[i]*yw2[i];
00653                       }
00654                       rat=ssa/ss2;
00655                       if(rat <= 0.) {LOGERROR("SiPixelTemplateSplit") << "illegal chi2ymin normalization = " << rat << ENDL; rat = 1.;}
00656                       chi2y=ss2-2.f*ssa/rat+sa2/(rat*rat);
00657                       if(chi2y < chi2ymin) {
00658                               chi2ymin = chi2y;
00659                               minbinj = j;
00660                               minbink = k;
00661                       }
00662               } 
00663            }
00664        deltaj /= 2;
00665            if(minbinj > fybin) {jmin = minbinj - deltaj;} else {jmin = fybin;}
00666            if(minbinj < lybin) {jmax = minbinj + deltaj;} else {jmax = lybin;}
00667            if(minbink > fybin) {kmin = minbink - deltaj;} else {kmin = fybin;}
00668            if(minbink < lybin) {kmax = minbink + deltaj;} else {kmax = lybin;}          
00669         }
00670         
00671         if (theVerboseLevel > 9) {
00672        LOGDEBUG("SiPixelTemplateReco") <<
00673         "minbins " << minbinj << "," << minbink << " chi2ymin = " << chi2ymin << ENDL;
00674     }
00675         
00676 // Do not apply final template pass to 1-pixel clusters (use calibrated offset) 
00677         
00678         if(logypx == 1) {
00679         
00680            if(nypix ==1) {
00681               delta = templ.dyone();
00682                   sigma = templ.syone();
00683            } else {
00684               delta = templ.dytwo();
00685                   sigma = templ.sytwo();
00686            }
00687            
00688            yrec1 = 0.5f*(fypix+lypix-2*shifty+2.f*originy)*ysize-delta;
00689            yrec2 = yrec1;
00690            
00691            if(sigma <= 0.) {
00692               sigmay = ysize/sqrt12;
00693            } else {
00694           sigmay = sigma;
00695            }
00696            
00697 // Do probability calculation for one-pixel clusters
00698 
00699                 chi21max = fmax(chi21min, (double)templ.chi2yminone());
00700                 chi2ymin -=chi21max;
00701                 if(chi2ymin < 0.) {chi2ymin = 0.;}
00702                 //         prob2y = gsl_cdf_chisq_Q(chi2ymin, mean1pix);
00703                 meany = fmax(mean1pix, (double)templ.chi2yavgone());
00704                 hchi2 = chi2ymin/2.; hndof = meany/2.;
00705                 prob2y = 1. - TMath::Gamma(hndof, hchi2);
00706            
00707         } else {
00708            
00709 // For cluster > 1 pix, use chi^2 minimm to recontruct the two y-positions
00710 
00711 // at small eta, the templates won't actually work on two pixel y-clusters so just return the pixel centers        
00712 
00713            if(logypx == 2 && fabsf(cotbeta) < 0.25f) {
00714                    switch(nypix) {
00715                       case 2:
00716 //  Both pixels are small
00717                      yrec1 = (fypix-shifty+originy)*ysize;
00718                          yrec2 = (lypix-shifty+originy)*ysize;
00719                          sigmay = ysize/sqrt12;
00720                                  break;
00721                           case 3:
00722 //  One big pixel and one small pixel
00723                             if(yd[fypix]) {
00724                                    yrec1 = (fypix+0.5f-shifty+originy)*ysize;
00725                                    yrec2 = (lypix-shifty+originy)*ysize;
00726                                    sigmay = ysize/sqrt12;
00727                                 } else {
00728                                    yrec1 = (fypix-shifty+originy)*ysize;
00729                                    yrec2 = (lypix-0.5f-shifty+originy)*ysize;
00730                                    sigmay = 1.5f*ysize/sqrt12;
00731                                 }
00732                                 break;
00733                          case 4:
00734 //  Two big pixels
00735                                 yrec1 = (fypix+0.5f-shifty+originy)*ysize;
00736                                 yrec2 = (lypix-0.5f-shifty+originy)*ysize;
00737                                 sigmay = 2.f*ysize/sqrt12;
00738                             break;
00739                          default:
00740 //  Something is screwy ...
00741                     LOGERROR("SiPixelTemplateReco") << "weird problem: logical y-pixels = " << logypx << ", total ysize in normal pixels = " << nypix << ENDL;  
00742                     return 10;
00743               }
00744            } else {
00745            
00746 // uncertainty and final correction depend upon charge bin       
00747   
00748               bias = templ.yavgc2m(binq);
00749               yrec1 = (0.125f*(minbink-ycbin)+BHY-(float)shifty+originy)*ysize - bias;
00750               yrec2 = (0.125f*(minbinj-ycbin)+BHY-(float)shifty+originy)*ysize - bias;
00751               sigmay = sqrt2y*templ.yrmsc2m(binq);
00752                   
00753            }
00754            
00755 // Do goodness of fit test in y  
00756            
00757            if(chi2ymin < 0.0) {chi2ymin = 0.0;}
00758            meany = templ.chi2yavgc2m(binq);
00759            if(meany < 0.01) {meany = 0.01;}
00760 // gsl function that calculates the chi^2 tail prob for non-integral dof
00761 //         prob2y = gsl_cdf_chisq_Q(chi2y, meany);
00762 //         prob2y = ROOT::Math::chisquared_cdf_c(chi2y, meany);
00763        hchi2 = chi2ymin/2.; hndof = meany/2.;
00764            prob2y = 1. - TMath::Gamma(hndof, hchi2);
00765         }
00766         
00767 // Do the x-reconstruction next 
00768                           
00769                                   
00770 // Define the maximum signal to allow before de-weighting a pixel 
00771 
00772         sxthr = 2.1f*templ.sxmax();
00773                           
00774 // Make sure that there will be at least two pixels that are not de-weighted 
00775         std::sort(&xsort[0], &xsort[logxpx]);
00776         if(logxpx == 1) {sxthr = 1.01f*xsort[0];} else {
00777            if (xsort[1] > sxthr) { sxthr = 1.01f*xsort[1]; }
00778         }
00779            
00780 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
00781 
00782 //      for(i=0; i<BYSIZE; ++i) { xsig2[i] = 0.; }
00783         templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
00784                           
00785 // Find the template bin that minimizes the Chi^2 
00786 
00787         chi2xmin = 1.e15;
00788         ss2 = 0.f;
00789         for(i=fxpix-2; i<=lxpix+2; ++i) {
00790                 xw2[i] = 1.f/xsig2[i];
00791                 xsw[i] = xsum[i]*xw2[i];
00792                 ss2 += xsum[i]*xsw[i];
00793         }
00794         minbinj = -1; 
00795         minbink = -1;
00796         deltaj = djx;
00797         jmin = fxbin;
00798         jmax = lxbin;
00799         kmin = fxbin;
00800         kmax = lxbin;
00801     std::vector<float> xtemp(BXSIZE);
00802         while(deltaj > 0) {
00803            for(j=jmin; j<jmax; j+=deltaj) {
00804                         km = std::min(kmax, j);
00805               for(k=kmin; k<=km; k+=deltaj) {
00806               
00807 // Get the template for this set of indices
00808               
00809               templ.xtemp3d(j, k, xtemp);
00810               
00811 // Modify the template if double pixels are present   
00812               
00813               if(nxpix > logxpx) {
00814                   i=fxpix;
00815                   while(i < lxpix) {
00816                      if(xd[i] && !xd[i+1]) {
00817                                   
00818 // Sum the adjacent cells and put the average signal in both   
00819                                   
00820                         sigavg = (xtemp[i] +  xtemp[i+1])/2.f;
00821                         xtemp[i] = sigavg;
00822                         xtemp[i+1] = sigavg;
00823                           i += 2;
00824                       } else {
00825                           ++i;
00826                       }
00827                   }
00828               } 
00829                       ssa = 0.f;
00830                       sa2 = 0.f;
00831                       for(i=fxpix-2; i<=lxpix+2; ++i) {
00832                               ssa += xsw[i]*xtemp[i];
00833                               sa2 += xtemp[i]*xtemp[i]*xw2[i];
00834                       }
00835                       rat=ssa/ss2;
00836                       if(rat <= 0.f) {LOGERROR("SiPixelTemplateSplit") << "illegal chi2xmin normalization = " << rat << ENDL; rat = 1.;}
00837                       chi2x=ss2-2.f*ssa/rat+sa2/(rat*rat);
00838                       if(chi2x < chi2xmin) {
00839                               chi2xmin = chi2x;
00840                               minbinj = j;
00841                               minbink = k;
00842                                 }
00843                    }
00844            } 
00845        deltaj /= 2;
00846            if(minbinj > fxbin) {jmin = minbinj - deltaj;} else {jmin = fxbin;}
00847            if(minbinj < lxbin) {jmax = minbinj + deltaj;} else {jmax = lxbin;}
00848            if(minbink > fxbin) {kmin = minbink - deltaj;} else {kmin = fxbin;}
00849            if(minbink < lxbin) {kmax = minbink + deltaj;} else {kmax = lxbin;}          
00850         }
00851         
00852         if (theVerboseLevel > 9) {
00853        LOGDEBUG("SiPixelTemplateSplit") <<
00854         "minbinj/minbink " << minbinj<< "/" << minbink << " chi2xmin = " << chi2xmin << ENDL;
00855     }
00856 
00857 // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
00858         
00859         if(logxpx == 1) {
00860         
00861            if(nxpix ==1) {
00862               delta = templ.dxone();
00863                   sigma = templ.sxone();
00864            } else {
00865               delta = templ.dxtwo();
00866                   sigma = templ.sxtwo();
00867            }
00868            xrec1 = 0.5f*(fxpix+lxpix-2*shiftx+2.f*originx)*xsize-delta;
00869            xrec2 = xrec1;
00870            if(sigma <= 0.f) {
00871               sigmax = xsize/sqrt12;
00872            } else {
00873           sigmax = sigma;
00874            }
00875            
00876 // Do probability calculation for one-pixel clusters
00877 
00878                 chi21max = fmax(chi21min, (double)templ.chi2xminone());
00879                 chi2xmin -=chi21max;
00880                 if(chi2xmin < 0.) {chi2xmin = 0.;}
00881                 meanx = fmax(mean1pix, (double)templ.chi2xavgone());
00882                 hchi2 = chi2xmin/2.; hndof = meanx/2.;
00883                 prob2x = 1. - TMath::Gamma(hndof, hchi2);
00884            
00885         } else {
00886            
00887 // For cluster > 1 pix, use chi^2 minimm to recontruct the two x-positions
00888 
00889 // uncertainty and final correction depend upon charge bin         
00890            
00891            bias = templ.xavgc2m(binq);
00892            k = std::min(minbink, minbinj);
00893            j = std::max(minbink, minbinj);
00894        xrec1 = (0.125f*(minbink-xcbin)+BHX-(float)shiftx+originx)*xsize - bias;
00895        xrec2 = (0.125f*(minbinj-xcbin)+BHX-(float)shiftx+originx)*xsize - bias;
00896            sigmax = sqrt2x*templ.xrmsc2m(binq);
00897            
00898 // Do goodness of fit test in y  
00899            
00900            if(chi2xmin < 0.0) {chi2xmin = 0.0;}
00901            meanx = templ.chi2xavgc2m(binq);
00902            if(meanx < 0.01) {meanx = 0.01;}
00903        hchi2 = chi2xmin/2.; hndof = meanx/2.;
00904            prob2x = 1. - TMath::Gamma(hndof, hchi2);
00905         }
00906         
00907         //  Don't return exact zeros for the probability
00908         
00909         if(prob2y < probmin) {prob2y = probmin;}
00910         if(prob2x < probmin) {prob2x = probmin;}
00911         
00912 //  New code: resolve the ambiguity if resolve is set to true
00913         
00914         dchisq = 0.;
00915         if(resolve) {
00916                 
00917 //  First copy the unexpanded cluster into a new working buffer
00918                 
00919                 for(i=0; i<BYM2; ++i) {
00920                         for(j=0; j<BXM2; ++j) {
00921                                 if((i>0 && i<BYM3)&&(j>0 && j<BXM3)) {
00922                                         cluster2[j][i] = qscale*clust[j-1][i-1];
00923                                 } else {
00924                                         cluster2[j][i] = 0.f;
00925                                 }
00926                         }
00927                 }
00928                 
00929 //  Next, redefine the local coordinates to start at the lower LH corner of pixel[1][1];
00930                 
00931                 if(xdouble[0]) {
00932                         x1p = xrec1+xsize;
00933                         x2p = xrec2+xsize;
00934                 } else {
00935                         x1p = xrec1+xsize/2.f;
00936                         x2p = xrec2+xsize/2.f;  
00937                 }
00938                 
00939                 if(ydouble[0]) {
00940                         y1p = yrec1+ysize;
00941                         y2p = yrec2+ysize;
00942                 } else {
00943                         y1p = yrec1+ysize/2.f;
00944                         y2p = yrec2+ysize/2.f;  
00945                 }
00946                 
00947 // Next, calculate 2-d templates for the (x1,y1)+(x2,y2) and the (x1,y2)+(x2,y1) hypotheses
00948                 
00949 //  First zero the two 2-d hypotheses
00950                 
00951                 for(i=0; i<BYM2; ++i) {
00952                         for(j=0; j<BXM2; ++j) {
00953                                 temp2d1[j][i] = 0.f;
00954                                 temp2d2[j][i] = 0.f;
00955                         }
00956                 }
00957                 
00958                 
00959 // Add the two hits in the first hypothesis
00960                 
00961                 any2dfail = templ2D.xytemp(id, cotalpha, cotbeta, x1p, y1p, xdouble, ydouble, temp2d1) && templ2D.xytemp(id, cotalpha, cotbeta, x2p, y2p, xdouble, ydouble, temp2d1);
00962                 
00963 // And then the second hypothesis
00964                 
00965                 any2dfail = any2dfail && templ2D.xytemp(id, cotalpha, cotbeta, x1p, y2p, xdouble, ydouble, temp2d2);
00966                 
00967                 any2dfail = any2dfail && templ2D.xytemp(id, cotalpha, cotbeta, x2p, y1p, xdouble, ydouble, temp2d2);
00968                 
00969 // If any of these have failed, use the simple templates instead
00970                 
00971                 if(!any2dfail) {
00972                         
00973 //  Rezero the two 2-d hypotheses
00974                         
00975                         for(i=0; i<BYM2; ++i) {
00976                                 for(j=0; j<BXM2; ++j) {
00977                                         temp2d1[j][i] = 0.f;
00978                                         temp2d2[j][i] = 0.f;
00979                                 }
00980                         }
00981                         
00982 // Add the two hits in the first hypothesis
00983                         
00984                         if(!templ.simpletemplate2D(x1p, y1p, xdouble, ydouble, temp2d1)) {return 1;}
00985                         
00986                         if(!templ.simpletemplate2D(x2p, y2p, xdouble, ydouble, temp2d1)) {return 1;}
00987                         
00988 // And then the second hypothesis
00989                         
00990                         if(!templ.simpletemplate2D(x1p, y2p, xdouble, ydouble, temp2d2)) {return 1;}
00991                         
00992                         if(!templ.simpletemplate2D(x2p, y1p, xdouble, ydouble, temp2d2)) {return 1;}
00993                 }
00994                 
00995 // Keep lists of pixels and nearest neighbors   
00996                 
00997                 std::list<std::pair<int, int> > pixellst;
00998                 
00999 // Loop through the array and find non-zero elements
01000                 
01001                 for(i=0; i<BYM2; ++i) {
01002                         for(j=0; j<BXM2; ++j) {
01003                                 if(cluster2[j][i] > 0.f || temp2d1[j][i] > 0.f || temp2d2[j][i] > 0.f) {
01004                                         pixel.first = j; pixel.second = i;
01005                                         pixellst.push_back(pixel);
01006                                 }
01007                         }
01008                 }
01009                 
01010                 // Now calculate the product of Landau probabilities (alpha probability)
01011                 
01012       templ2D.landau_par(lanpar);               
01013                 loglike1 = 0.; loglike2 = 0.;
01014                 
01015                 // Now, for each neighbor, match it again the pixel list.
01016                 // If found, delete it from the neighbor list
01017                 
01018                 std::list<std::pair<int, int> >::const_iterator pixIter, pixEnd;
01019                 pixIter = pixellst.begin();
01020                 pixEnd = pixellst.end();
01021                 for( ; pixIter != pixEnd; ++pixIter ) {
01022                         j = pixIter->first; i = pixIter->second;
01023                         if((i < BHY && cotbeta > 0.) || (i >= BHY && cotbeta < 0.)) {lparm = 0;} else {lparm = 1;}
01024                         mpv1 = lanpar[lparm][0] + lanpar[lparm][1]*temp2d1[j][i];
01025                         sigmal1 = lanpar[lparm][2]*mpv1;
01026                         sigmal1 = sqrt(sigmal1*sigmal1+lanpar[lparm][3]*lanpar[lparm][3]);
01027                         if(sigmal1 < q05) sigmal1 = q05;
01028                         arg1 = (cluster2[j][i]-mpv1)/sigmal1-0.22278;
01029                         mpv2 = lanpar[lparm][0] + lanpar[lparm][1]*temp2d2[j][i];
01030                         sigmal2 = lanpar[lparm][2]*mpv2;
01031                         sigmal2 = sqrt(sigmal2*sigmal2+lanpar[lparm][3]*lanpar[lparm][3]);
01032                         if(sigmal2 < q05) sigmal2 = q05;
01033                         arg2 = (cluster2[j][i]-mpv2)/sigmal2-0.22278;
01034 //                      like = ROOT::Math::landau_pdf(arg1)/sigmal1;
01035                         like = ROOT::Math::landau_pdf(arg1);
01036                         if(like < 1.e-30) like = 1.e-30;
01037                         loglike1 += log(like);
01038 //                      like = ROOT::Math::landau_pdf(arg2)/sigmal2;
01039                         like = ROOT::Math::landau_pdf(arg2);
01040                         if(like < 1.e-30) like = 1.e-30;
01041                         loglike2 += log(like);
01042                 }
01043                 
01044 // Calculate chisquare difference for the two hypotheses 9don't multiply by 2 for less inconvenient scaling
01045                 
01046                 deltachi2 = loglike1 - loglike2;
01047                                 
01048                 if(deltachi2 < 0.) {
01049                         
01050                         // Flip the x1 and x2
01051                         
01052                         x1p = xrec1;
01053                         xrec1 = xrec2;
01054                         xrec2 = x1p;
01055                 }
01056                 
01057                 // Return a positive definite value
01058                 
01059                 dchisq = fabs(deltachi2);               
01060         }
01061         
01062     return 0;
01063 } // PixelTempSplit 
01064 
01065 
01066 // *************************************************************************************************************************************
01089 // *************************************************************************************************************************************
01090 
01091 int SiPixelTemplateSplit::PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d& cluster, 
01092                     std::vector<bool>& ydouble, std::vector<bool>& xdouble, 
01093                     SiPixelTemplate& templ, 
01094                     float& yrec1, float& yrec2, float& sigmay, float& prob2y,
01095                         float& xrec1, float& xrec2, float& sigmax, float& prob2x, int& q2bin, float& prob2Q, bool resolve, int speed, float& dchisq, SiPixelTemplate2D& templ2D)
01096 {
01097     // Local variables 
01098         const bool deadpix = false;
01099         std::vector<std::pair<int, int> > zeropix;
01100     
01101         return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ, 
01102                     yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);
01103 
01104 } // PixelTempSplit
01105 
01106 
01107 // *************************************************************************************************************************************
01130 // *************************************************************************************************************************************
01131 
01132 int SiPixelTemplateSplit::PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d& cluster, 
01133                                          std::vector<bool>& ydouble, std::vector<bool>& xdouble, 
01134                                          SiPixelTemplate& templ, 
01135                                          float& yrec1, float& yrec2, float& sigmay, float& prob2y,
01136                                          float& xrec1, float& xrec2, float& sigmax, float& prob2x, int& q2bin, float& prob2Q, bool resolve, float& dchisq, SiPixelTemplate2D& templ2D)
01137 {
01138     // Local variables 
01139         const bool deadpix = false;
01140         std::vector<std::pair<int, int> > zeropix;
01141     const int speed = 1;
01142     
01143         return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ, 
01144                                                 yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);
01145     
01146 } // PixelTempSplit
01147 
01148 // *************************************************************************************************************************************
01169 // *************************************************************************************************************************************
01170 
01171 int SiPixelTemplateSplit::PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d& cluster, 
01172                                                                                 std::vector<bool>& ydouble, std::vector<bool>& xdouble, 
01173                                                                                 SiPixelTemplate& templ, 
01174                                                                                 float& yrec1, float& yrec2, float& sigmay, float& prob2y,
01175                                                                                 float& xrec1, float& xrec2, float& sigmax, float& prob2x, int& q2bin, float& prob2Q, SiPixelTemplate2D& templ2D)
01176 {
01177     // Local variables 
01178         const bool deadpix = false;
01179         const bool resolve = true;
01180         float dchisq;
01181         std::vector<std::pair<int, int> > zeropix;
01182         const int speed = 1;
01183     
01184         return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ, 
01185                                                                                            yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);
01186         
01187 } // PixelTempSplit
01188 
01189 
01190 
01191 // *************************************************************************************************************************************
01211 // *************************************************************************************************************************************
01212 
01213 int SiPixelTemplateSplit::PixelTempSplit(int id, float cotalpha, float cotbeta, array_2d& cluster, 
01214                                                                                                          std::vector<bool>& ydouble, std::vector<bool>& xdouble, 
01215                                                                                                          SiPixelTemplate& templ, 
01216                                                                                                          float& yrec1, float& yrec2, float& sigmay, float& prob2y,
01217                                                                                                          float& xrec1, float& xrec2, float& sigmax, float& prob2x, int& q2bin, SiPixelTemplate2D& templ2D)
01218 {
01219         // Local variables 
01220         const bool deadpix = false;
01221         const bool resolve = true;
01222         float dchisq, prob2Q;
01223         std::vector<std::pair<int, int> > zeropix;
01224         const int speed = 1;
01225         
01226         return SiPixelTemplateSplit::PixelTempSplit(id, cotalpha, cotbeta, cluster, ydouble, xdouble, templ, 
01227                                                                                                                          yrec1, yrec2, sigmay, prob2y, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q, resolve, speed, dchisq, deadpix, zeropix, templ2D);
01228         
01229 } // PixelTempSplit
01230 
01231
/data/refman/pasoursint/CMSSW_5_3_9_patch3/src/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateSplit.cc
