Reference Manual

Go to the documentation of this file.
00001 //
00002 //  SiPixelTemplateSplit.cc (Version 1.05)
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 //  Copyright 2008 __TheJohnsHopkinsUniversity__. All rights reserved.
00009 //
00010 //  Incorporate "cluster repair" to handle dead pixels
00011 //  Take truncation size from new pixmax information
00012 //  Change to allow template sizes to be changed at compile time
00013 //  Move interpolation range error to LogDebug
00014 //  Add qbin = 5 and change 1-pixel probability to use new template info
00015 //  Add floor for probabilities (no exact zeros)
00016 //
00017 
00018 #include <math.h>
00019 #include <algorithm>
00020 #include <vector>
00021 #include <iostream>
00022 // ROOT::Math has a c++ function that does the probability calc, but only in v5.12 and later
00023 //#include "Math/DistFunc.h"
00024 #include "TMath.h"
00025 // Use current version of gsl instead of ROOT::Math
00026 //#include <gsl/gsl_cdf.h>
00027 
00028 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00029 #include "SiPixelTemplateSplit.h"
00030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00031 #define LOGERROR(x) edm::LogError(x)
00032 #define LOGDEBUG(x) LogDebug(x)
00033 static int theVerboseLevel = 2;
00034 #define ENDL " "
00035 #else
00036 #include "SiPixelTemplateSplit.h"
00037 //static int theVerboseLevel = {2};
00038 #define LOGERROR(x) std::cout << x << ": "
00039 #define LOGDEBUG(x) std::cout << x << ": "
00040 #define ENDL std::endl
00041 #endif
00042 
00043 using namespace SiPixelTemplateReco;
00044 
00045 // *************************************************************************************************************************************
00069 // *************************************************************************************************************************************
00070 int SiPixelTemplateReco::PixelTempSplit(int id, bool fpix, float cotalpha, float cotbeta, array_2d cluster, 
00071                     std::vector<bool> ydouble, std::vector<bool> xdouble, 
00072                     SiPixelTemplate& templ, 
00073                     float& yrec1, float& yrec2, float& sigmay, float& proby,
00074                         float& xrec1, float& xrec2, float& sigmax, float& probx, int& qbin, bool deadpix, std::vector<std::pair<int, int> > zeropix)
00075                         
00076 {
00077     // Local variables 
00078         static int i, j, k, minbin, binq, midpix, fypix, nypix, lypix, logypx;
00079     static int fxpix, nxpix, lxpix, logxpx, shifty, shiftx, nyzero[TYSIZE];
00080         static int nclusx, nclusy;
00081         static int nybin, ycbin, nxbin, xcbin, minbinj, minbink;
00082         static float sythr, sxthr, delta, sigma, sigavg, pseudopix, qscale;
00083         static float ss2, ssa, sa2, rat, fq, qtotal, qpixel;
00084         static float originx, originy, bias, maxpix, minmax;
00085         static double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
00086         static double hchi2, hndof;
00087         static float ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
00088         static float ysig2[BYSIZE], xsig2[BXSIZE];
00089         static bool yd[BYSIZE], xd[BXSIZE], anyyd, anyxd;
00090         const float ysize={150.}, xsize={100.}, sqrt2={2.};
00091         const float probmin={1.110223e-16};
00092 //  const float sqrt2={1.41421356};
00093         
00094 // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
00095 
00096         const double mean1pix={0.100}, chi21min={0.160};
00097                       
00098 // First, interpolate the template needed to analyze this cluster     
00099 // check to see of the track direction is in the physical range of the loaded template
00100 
00101         if(!templ.interpolate(id, fpix, cotalpha, cotbeta)) {
00102            LOGDEBUG("SiPixelTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha << ", cot(beta) = " << cotbeta << " is not within the acceptance of fpix = "
00103            << fpix << ", template ID = " << id << ", no reconstruction performed" << ENDL;      
00104            return 20;
00105         }
00106                       
00107 // Define size of pseudopixel
00108 
00109     pseudopix = templ.s50();
00110         
00111 // Get charge scaling factor
00112 
00113         qscale = templ.qscale();
00114     
00115 // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
00116 
00117         if(cluster.num_dimensions() != 2) {
00118            LOGERROR("SiPixelTemplateReco") << "input cluster container (BOOST Multiarray) has wrong number of dimensions" << ENDL;      
00119            return 3;
00120         }
00121         nclusx = (int)cluster.shape()[0];
00122         nclusy = (int)cluster.shape()[1];
00123         if(nclusx != (int)xdouble.size()) {
00124            LOGERROR("SiPixelTemplateReco") << "input cluster container x-size is not equal to double pixel flag container size" << ENDL;        
00125            return 4;
00126         }
00127         if(nclusy != (int)ydouble.size()) {
00128            LOGERROR("SiPixelTemplateReco") << "input cluster container y-size is not equal to double pixel flag container size" << ENDL;        
00129            return 5;
00130         }
00131         
00132 // enforce maximum size 
00133         
00134         if(nclusx > TXSIZE) {nclusx = TXSIZE;}
00135         if(nclusy > TYSIZE) {nclusy = TYSIZE;}
00136         
00137 // First, rescale all pixel charges       
00138 
00139     for(i=0; i<nclusy; ++i) {
00140            for(j=0; j<nclusx; ++j) {
00141                   if(cluster[j][i] > 0) {cluster[j][i] *= qscale;}
00142            }
00143         }
00144         
00145 // Next, sum the total charge and "decapitate" big pixels         
00146 
00147         qtotal = 0.;
00148         minmax = 2.*templ.pixmax();
00149         for(i=0; i<nclusy; ++i) {
00150            maxpix = minmax;
00151            if(ydouble[i]) {maxpix *=2.;}
00152            for(j=0; j<nclusx; ++j) {
00153                   qtotal += cluster[j][i];
00154                   if(cluster[j][i] > maxpix) {cluster[j][i] = maxpix;}
00155            }
00156         }
00157         
00158 // Do the cluster repair here   
00159         
00160     if(deadpix) {
00161            fypix = BYM3; lypix = -1;
00162        for(i=0; i<nclusy; ++i) {
00163               ysum[i] = 0.; nyzero[i] = 0;
00164 // Do preliminary cluster projection in y
00165               for(j=0; j<nclusx; ++j) {
00166                      ysum[i] += cluster[j][i];
00167                   }
00168                   if(ysum[i] > 0.) {
00169 // identify ends of cluster to determine what the missing charge should be
00170                      if(i < fypix) {fypix = i;}
00171                          if(i > lypix) {lypix = i;}
00172                   }
00173            }
00174            
00175 // Now loop over dead pixel list and "fix" everything   
00176 
00177 //First see if the cluster ends are redefined and that we have only one dead pixel per column
00178 
00179            std::vector<std::pair<int, int> >::const_iterator zeroIter = zeropix.begin(), zeroEnd = zeropix.end();
00180        for ( ; zeroIter != zeroEnd; ++zeroIter ) {
00181               i = zeroIter->second;
00182                   if(i<0 || i>TYSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column y-index " << i << ", no reconstruction performed" << ENDL;       
00183                                return 11;}
00184                                                    
00185 // count the number of dead pixels in each column
00186                   ++nyzero[i];
00187 // allow them to redefine the cluster ends
00188                   if(i < fypix) {fypix = i;}
00189                   if(i > lypix) {lypix = i;}
00190            }
00191            
00192            nypix = lypix-fypix+1;
00193            
00194 // 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
00195            
00196        for (zeroIter = zeropix.begin(); zeroIter != zeroEnd; ++zeroIter ) {        
00197               i = zeroIter->second; j = zeroIter->first;
00198                   if(j<0 || j>TXSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column x-index " << j << ", no reconstruction performed" << ENDL;       
00199                                return 12;}
00200                   if((i == fypix || i == lypix) && nypix > 1) {maxpix = templ.symax()/2.;} else {maxpix = templ.symax();}
00201                   if(ydouble[i]) {maxpix *=2.;}
00202                   if(nyzero[i] > 0 && nyzero[i] < 3) {qpixel = (maxpix - ysum[i])/(float)nyzero[i];} else {qpixel = 1.;}
00203                   if(qpixel < 1.) {qpixel = 1.;}
00204           cluster[j][i] = qpixel;
00205 // Adjust the total cluster charge to reflect the charge of the "repaired" cluster
00206                   qtotal += qpixel;
00207            }
00208 // End of cluster repair section
00209         } 
00210         
00211 // Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container         
00212 
00213     for(i=0; i<BYSIZE; ++i) { ysum[i] = 0.; yd[i] = false;}
00214         k=0;
00215         anyyd = false;
00216     for(i=0; i<nclusy; ++i) {
00217            for(j=0; j<nclusx; ++j) {
00218                   ysum[k] += cluster[j][i];
00219            }
00220     
00221 // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels  
00222    
00223            if(ydouble[i]) {
00224               ysum[k] /= 2.;
00225                   ysum[k+1] = ysum[k];
00226                   yd[k] = true;
00227                   yd[k+1] = false;
00228                   k=k+2;
00229                   anyyd = true;
00230            } else {
00231                   yd[k] = false;
00232               ++k;
00233            }
00234            if(k > BYM1) {break;}
00235         }
00236                  
00237 // Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container         
00238 
00239     for(i=0; i<BXSIZE; ++i) { xsum[i] = 0.; xd[i] = false;}
00240         k=0;
00241         anyxd = false;
00242     for(j=0; j<nclusx; ++j) {
00243            for(i=0; i<nclusy; ++i) {
00244                   xsum[k] += cluster[j][i];
00245            }
00246     
00247 // If this is a double pixel, put 1/2 of the charge in 2 consective single pixels  
00248    
00249            if(xdouble[j]) {
00250               xsum[k] /= 2.;
00251                   xsum[k+1] = xsum[k];
00252                   xd[k]=true;
00253                   xd[k+1]=false;
00254                   k=k+2;
00255                   anyxd = true;
00256            } else {
00257                   xd[k]=false;
00258               ++k;
00259            }
00260            if(k > BXM1) {break;}
00261         }
00262         
00263 // next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx   
00264 
00265     fypix=-1;
00266         nypix=0;
00267         lypix=0;
00268         logypx=0;
00269         for(i=0; i<BYSIZE; ++i) {
00270            if(ysum[i] > 0.) {
00271               if(fypix == -1) {fypix = i;}
00272                   if(!yd[i]) {
00273                      ysort[logypx] = ysum[i];
00274                          ++logypx;
00275                   }
00276                   ++nypix;
00277                   lypix = i;
00278                 }
00279         }
00280         
00281 // Make sure cluster is continuous
00282 
00283         if((lypix-fypix+1) != nypix || nypix == 0) { 
00284            LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
00285            if (theVerboseLevel > 2) {
00286           LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
00287           for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}           
00288                   LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
00289        }
00290         
00291            return 1; 
00292         }
00293         
00294 // If cluster is longer than max template size, technique fails
00295 
00296         if(nypix > TYSIZE) { 
00297            LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster is larger than maximum template size" << ENDL;
00298            if (theVerboseLevel > 2) {
00299           LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
00300           for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}           
00301                   LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
00302        }
00303         
00304            return 6; 
00305         }
00306         
00307 // next, center the cluster on pixel 12 if necessary   
00308 
00309         midpix = (fypix+lypix)/2;
00310         shifty = BHY - midpix;
00311         if(shifty > 0) {
00312            for(i=lypix; i>=fypix; --i) {
00313               ysum[i+shifty] = ysum[i];
00314                   ysum[i] = 0.;
00315                   yd[i+shifty] = yd[i];
00316                   yd[i] = false;
00317            }
00318         } else if (shifty < 0) {
00319            for(i=fypix; i<=lypix; ++i) {
00320               ysum[i+shifty] = ysum[i];
00321                   ysum[i] = 0.;
00322                   yd[i+shifty] = yd[i];
00323                   yd[i] = false;
00324            }
00325     }
00326         lypix +=shifty;
00327         fypix +=shifty;
00328         
00329 // If the cluster boundaries are OK, add pesudopixels, otherwise quit
00330         
00331         if(fypix > 1 && fypix < BYM2) {
00332            ysum[fypix-1] = pseudopix;
00333            ysum[fypix-2] = 0.2*pseudopix;
00334         } else {return 8;}
00335         if(lypix > 1 && lypix < BYM2) {
00336            ysum[lypix+1] = pseudopix;   
00337            ysum[lypix+2] = 0.2*pseudopix;
00338         } else {return 8;}
00339         
00340 // finally, determine if pixel[0] is a double pixel and make an origin correction if it is   
00341 
00342     if(ydouble[0]) {
00343            originy = -0.5;
00344         } else {
00345            originy = 0.;
00346         }
00347         
00348 // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx   
00349 
00350     fxpix=-1;
00351         nxpix=0;
00352         lxpix=0;
00353         logxpx=0;
00354         for(i=0; i<BXSIZE; ++i) {
00355            if(xsum[i] > 0.) {
00356               if(fxpix == -1) {fxpix = i;}
00357                   if(!xd[i]) {
00358                      xsort[logxpx] = xsum[i];
00359                          ++logxpx;
00360                   }
00361                   ++nxpix;
00362                   lxpix = i;
00363                 }
00364         }
00365         
00366 // Make sure cluster is continuous
00367 
00368         if((lxpix-fxpix+1) != nxpix) { 
00369         
00370            LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
00371            if (theVerboseLevel > 2) {
00372           LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
00373           for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}           
00374                   LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
00375        }
00376 
00377            return 2; 
00378         }
00379 
00380 // If cluster is longer than max template size, technique fails
00381 
00382         if(nxpix > TXSIZE) { 
00383         
00384            LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
00385            if (theVerboseLevel > 2) {
00386           LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
00387           for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}           
00388                   LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
00389        }
00390 
00391            return 7; 
00392         }
00393         
00394 // next, center the cluster on pixel 5 if necessary   
00395 
00396         midpix = (fxpix+lxpix)/2;
00397         shiftx = BHX - midpix;
00398         if(shiftx > 0) {
00399            for(i=lxpix; i>=fxpix; --i) {
00400               xsum[i+shiftx] = xsum[i];
00401                   xsum[i] = 0.;
00402               xd[i+shiftx] = xd[i];
00403                   xd[i] = false;
00404            }
00405         } else if (shiftx < 0) {
00406            for(i=fxpix; i<=lxpix; ++i) {
00407               xsum[i+shiftx] = xsum[i];
00408                   xsum[i] = 0.;
00409               xd[i+shiftx] = xd[i];
00410                   xd[i] = false;
00411            }
00412     }
00413         lxpix +=shiftx;
00414         fxpix +=shiftx;
00415         
00416 // If the cluster boundaries are OK, add pesudopixels, otherwise quit
00417         
00418         if(fxpix > 1 && fxpix <BXM2) {
00419            xsum[fxpix-1] = pseudopix;
00420            xsum[fxpix-2] = 0.2*pseudopix;
00421         } else {return 9;}
00422         if(lxpix > 1 && lxpix < BXM2) {
00423            xsum[lxpix+1] = pseudopix;
00424            xsum[lxpix+2] = 0.2*pseudopix;
00425         } else {return 9;}
00426                         
00427 // finally, determine if pixel[0] is a double pixel and make an origin correction if it is   
00428 
00429     if(xdouble[0]) {
00430            originx = -0.5;
00431         } else {
00432            originx = 0.;
00433         }
00434         
00435 // uncertainty and final corrections depend upon total charge bin          
00436            
00437         fq = qtotal/templ.qavg();
00438         if(fq > 3.0) {
00439            binq=0;
00440         } else {
00441            if(fq > 2.0) {
00442               binq=1;
00443            } else {
00444                   if(fq > 1.70) {
00445                          binq=2;
00446                   } else {
00447                          binq=3;
00448                   }
00449            }
00450         }
00451         
00452 // Return the charge bin via the parameter list unless the charge is too small (then flag it)
00453         
00454         qbin = binq;
00455         if(!deadpix && qtotal < 1.9*templ.qmin()) {qbin = 5;} else {
00456                 if(!deadpix && qtotal < 1.9*templ.qmin(1)) {qbin = 4;}
00457         }
00458         
00459         if (theVerboseLevel > 9) {
00460        LOGDEBUG("SiPixelTemplateReco") <<
00461         "ID = " << id << " FPix = " << fpix << 
00462          " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta << 
00463          " nclusx = " << nclusx << " nclusy = " << nclusy << ENDL;
00464     }
00465 
00466         
00467 // Next, generate the 3d y- and x-templates
00468    
00469     array_3d ytemp;     
00470         templ.ytemp3d(nypix, ytemp);
00471    
00472     array_3d xtemp;
00473         templ.xtemp3d(nxpix, xtemp);
00474         
00475 // Do the y-reconstruction first 
00476                                         
00477 // retrieve the number of y-bins 
00478 
00479     nybin = ytemp.shape()[0]; ycbin = nybin/2;
00480                           
00481 // Modify the template if double pixels are present   
00482         
00483         if(nypix > logypx) {
00484                 i=fypix;
00485                 while(i < lypix) {
00486                    if(yd[i] && !yd[i+1]) {
00487                           for(j=0; j<nybin; ++j) {
00488                              for(k=0; k<=j; ++k) {
00489                 
00490 // Sum the adjacent cells and put the average signal in both   
00491 
00492                                     sigavg = (ytemp[j][k][i] +  ytemp[j][k][i+1])/2.;
00493                                     ytemp[j][k][i] = sigavg;
00494                                     ytemp[j][k][i+1] = sigavg;
00495                                   }
00496                            }
00497                            i += 2;
00498                         } else {
00499                            ++i;
00500                         }
00501                  }
00502         }       
00503              
00504 // Define the maximum signal to allow before de-weighting a pixel 
00505 
00506         sythr = 2.1*(templ.symax());
00507                           
00508 // Make sure that there will be at least two pixels that are not de-weighted 
00509 
00510         std::sort(&ysort[0], &ysort[logypx]);
00511         if(logypx == 1) {sythr = 1.01*ysort[0];} else {
00512            if (ysort[1] > sythr) { sythr = 1.01*ysort[1]; }
00513         }
00514         
00515 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
00516 
00517         for(i=0; i<BYSIZE; ++i) { ysig2[i] = 0.;}
00518         templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
00519                           
00520 // Find the template bin that minimizes the Chi^2 
00521 
00522         chi2ymin = 1.e15;
00523         minbinj = -1; 
00524         minbink = -1;
00525         for(j=0; j<nybin; ++j) {
00526            for(k=0; k<=j; ++k) {
00527               ss2 = 0.;
00528                   ssa = 0.;
00529                   sa2 = 0.;
00530                   for(i=fypix-2; i<=lypix+2; ++i) {
00531                          ss2 += ysum[i]*ysum[i]/ysig2[i];
00532                          ssa += ysum[i]*ytemp[j][k][i]/ysig2[i];
00533                          sa2 += ytemp[j][k][i]*ytemp[j][k][i]/ysig2[i];
00534                   }
00535                   rat=ssa/ss2;
00536                   if(rat <= 0.) {LOGERROR("SiPixelTemplateReco") << "illegal chi2ymin normalization = " << rat << ENDL; rat = 1.;}
00537                   chi2y=ss2-2.*ssa/rat+sa2/(rat*rat);
00538                   if(chi2y < chi2ymin) {
00539                           chi2ymin = chi2y;
00540                           minbinj = j;
00541                           minbink = k;
00542                   }
00543            } 
00544         }
00545         
00546         if (theVerboseLevel > 9) {
00547        LOGDEBUG("SiPixelTemplateReco") <<
00548         "minbins " << minbinj << "," << minbink << " chi2ymin = " << chi2ymin << ENDL;
00549     }
00550         
00551 // Do not apply final template pass to 1-pixel clusters (use calibrated offset) 
00552         
00553         if(logypx == 1) {
00554         
00555            if(nypix ==1) {
00556               delta = templ.dyone();
00557                   sigma = templ.syone();
00558            } else {
00559               delta = templ.dytwo();
00560                   sigma = templ.sytwo();
00561            }
00562            
00563            yrec1 = 0.5*(fypix+lypix-2*shifty+2.*originy)*ysize-delta;
00564            yrec2 = yrec1;
00565            
00566            if(sigma <= 0.) {
00567               sigmay = 43.3;
00568            } else {
00569           sigmay = sigma;
00570            }
00571            
00572 // Do probability calculation for one-pixel clusters
00573 
00574                 chi21max = fmax(chi21min, (double)templ.chi2yminone());
00575                 chi2ymin -=chi21max;
00576                 if(chi2ymin < 0.) {chi2ymin = 0.;}
00577                 //         proby = gsl_cdf_chisq_Q(chi2ymin, mean1pix);
00578                 meany = fmax(mean1pix, (double)templ.chi2yavgone());
00579                 hchi2 = chi2ymin/2.; hndof = meany/2.;
00580                 proby = 1. - TMath::Gamma(hndof, hchi2);
00581            
00582         } else {
00583            
00584 // For cluster > 1 pix, use chi^2 minimm to recontruct the two y-positions
00585 
00586 // at small eta, the templates won't actually work on two pixel y-clusters so just return the pixel centers        
00587 
00588            if(logypx == 2 && fabsf(cotbeta) < 0.25) {
00589                    switch(nypix) {
00590                       case 2:
00591 //  Both pixels are small
00592                      yrec1 = (fypix-shifty+originy)*ysize;
00593                          yrec2 = (lypix-shifty+originy)*ysize;
00594                          sigmay = 43.3;
00595                                  break;
00596                           case 3:
00597 //  One big pixel and one small pixel
00598                             if(yd[fypix]) {
00599                                    yrec1 = (fypix+0.5-shifty+originy)*ysize;
00600                                    yrec2 = (lypix-shifty+originy)*ysize;
00601                                    sigmay = 43.3;
00602                                 } else {
00603                                    yrec1 = (fypix-shifty+originy)*ysize;
00604                                    yrec2 = (lypix-0.5-shifty+originy)*ysize;
00605                                    sigmay = 65.;
00606                                 }
00607                                 break;
00608                          case 4:
00609 //  Two big pixels
00610                                 yrec1 = (fypix+0.5-shifty+originy)*ysize;
00611                                 yrec2 = (lypix-0.5-shifty+originy)*ysize;
00612                                 sigmay = 86.6;
00613                             break;
00614                          default:
00615 //  Something is screwy ...
00616                     LOGERROR("SiPixelTemplateReco") << "weird problem: logical y-pixels = " << logypx << ", total ysize in normal pixels = " << nypix << ENDL;  
00617                     return 10;
00618               }
00619            } else {
00620            
00621 // uncertainty and final correction depend upon charge bin       
00622   
00623               bias = templ.yavgc2m(binq);
00624               yrec1 = (0.125*(minbink-ycbin)+BHY-(float)shifty+originy)*ysize - bias;
00625               yrec2 = (0.125*(minbinj-ycbin)+BHY-(float)shifty+originy)*ysize - bias;
00626               sigmay = sqrt2*templ.yrmsc2m(binq);
00627                   
00628            }
00629            
00630 // Do goodness of fit test in y  
00631            
00632            if(chi2ymin < 0.0) {chi2ymin = 0.0;}
00633            meany = 2.*templ.chi2yavg(binq);
00634            if(meany < 0.01) {meany = 0.01;}
00635 // gsl function that calculates the chi^2 tail prob for non-integral dof
00636 //         proby = gsl_cdf_chisq_Q(chi2y, meany);
00637 //         proby = ROOT::Math::chisquared_cdf_c(chi2y, meany);
00638        hchi2 = chi2ymin/2.; hndof = meany/2.;
00639            proby = 1. - TMath::Gamma(hndof, hchi2);
00640         }
00641         
00642 // Do the x-reconstruction next 
00643                           
00644 // retrieve the number of x-bins 
00645 
00646     nxbin = xtemp.shape()[0]; xcbin = nxbin/2;
00647                           
00648 // Modify the template if double pixels are present   
00649         
00650         if(nxpix > logxpx) {
00651                 i=fxpix;
00652                 while(i < lxpix) {
00653                    if(xd[i] && !xd[i+1]) {
00654                           for(j=0; j<nxbin; ++j) {
00655                              for(k=0; k<=j; ++k) {
00656                 
00657 // Sum the adjacent cells and put the average signal in both   
00658 
00659                                     sigavg = (xtemp[j][k][i] +  xtemp[j][k][i+1])/2.;
00660                                     xtemp[j][k][i] = sigavg;
00661                                     xtemp[j][k][i+1] = sigavg;
00662                                   }
00663                            }
00664                            i += 2;
00665                         } else {
00666                            ++i;
00667                         }
00668                  }
00669         }       
00670                                   
00671 // Define the maximum signal to allow before de-weighting a pixel 
00672 
00673         sxthr = 2.1*templ.sxmax();
00674                           
00675 // Make sure that there will be at least two pixels that are not de-weighted 
00676         std::sort(&xsort[0], &xsort[logxpx]);
00677         if(logxpx == 1) {sxthr = 1.01*xsort[0];} else {
00678            if (xsort[1] > sxthr) { sxthr = 1.01*xsort[1]; }
00679         }
00680            
00681 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
00682 
00683         for(i=0; i<BYSIZE; ++i) { xsig2[i] = 0.; }
00684         templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
00685                           
00686 // Find the template bin that minimizes the Chi^2 
00687 
00688         chi2xmin = 1.e15;
00689         minbinj = -1; 
00690         minbink = -1;
00691         for(j=0; j<nxbin; ++j) {
00692            for(k=0; k<=j; ++k) {
00693               ss2 = 0.;
00694                   ssa = 0.;
00695                   sa2 = 0.;
00696                   for(i=fxpix-2; i<=lxpix+2; ++i) {
00697                          ss2 += xsum[i]*xsum[i]/xsig2[i];
00698                          ssa += xsum[i]*xtemp[j][k][i]/xsig2[i];
00699                          sa2 += xtemp[j][k][i]*xtemp[j][k][i]/xsig2[i];
00700                   }
00701                   rat=ssa/ss2;
00702                   if(rat <= 0.) {LOGERROR("SiPixelTemplateReco") << "illegal chi2ymin normalization = " << rat << ENDL; rat = 1.;}
00703                   chi2x=ss2-2.*ssa/rat+sa2/(rat*rat);
00704                   if(chi2x < chi2xmin) {
00705                           chi2xmin = chi2x;
00706                           minbinj = j;
00707                           minbink = k;
00708                   }
00709            } 
00710         }
00711         
00712         if (theVerboseLevel > 9) {
00713        LOGDEBUG("SiPixelTemplateReco") <<
00714         "minbin " << minbin << " chi2xmin = " << chi2xmin << ENDL;
00715     }
00716 
00717 // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
00718         
00719         if(logxpx == 1) {
00720         
00721            if(nxpix ==1) {
00722               delta = templ.dxone();
00723                   sigma = templ.sxone();
00724            } else {
00725               delta = templ.dxtwo();
00726                   sigma = templ.sxtwo();
00727            }
00728            xrec1 = 0.5*(fxpix+lxpix-2*shiftx+2.*originx)*xsize-delta;
00729            xrec2 = xrec1;
00730            if(sigma <= 0.) {
00731               sigmax = 28.9;
00732            } else {
00733           sigmax = sigma;
00734            }
00735            
00736 // Do probability calculation for one-pixel clusters
00737 
00738                 chi21max = fmax(chi21min, (double)templ.chi2xminone());
00739                 chi2xmin -=chi21max;
00740                 if(chi2xmin < 0.) {chi2xmin = 0.;}
00741                 meanx = fmax(mean1pix, (double)templ.chi2xavgone());
00742                 hchi2 = chi2xmin/2.; hndof = meanx/2.;
00743                 probx = 1. - TMath::Gamma(hndof, hchi2);
00744            
00745         } else {
00746            
00747 // For cluster > 1 pix, use chi^2 minimm to recontruct the two x-positions
00748 
00749 // uncertainty and final correction depend upon charge bin         
00750            
00751            bias = templ.xavgc2m(binq);
00752            k = std::min(minbink, minbinj);
00753            j = std::max(minbink, minbinj);
00754        xrec1 = (0.125*(minbink-xcbin)+BHX-(float)shiftx+originx)*xsize - bias;
00755        xrec2 = (0.125*(minbinj-xcbin)+BHX-(float)shiftx+originx)*xsize - bias;
00756            sigmax = sqrt2*templ.xrmsc2m(binq);
00757            
00758 // Do goodness of fit test in y  
00759            
00760            if(chi2xmin < 0.0) {chi2xmin = 0.0;}
00761            meanx = 2.*templ.chi2xavg(binq);
00762            if(meanx < 0.01) {meanx = 0.01;}
00763        hchi2 = chi2xmin/2.; hndof = meanx/2.;
00764            probx = 1. - TMath::Gamma(hndof, hchi2);
00765         }
00766         
00767         //  Don't return exact zeros for the probability
00768         
00769         if(proby < probmin) {proby = probmin;}
00770         if(probx < probmin) {probx = probmin;}
00771         
00772     return 0;
00773 } // PixelTempSplit 
00774 
00775 
00776 // *************************************************************************************************************************************
00798 // *************************************************************************************************************************************
00799 
00800 int SiPixelTemplateReco::PixelTempSplit(int id, bool fpix, float cotalpha, float cotbeta, array_2d cluster, 
00801                     std::vector<bool> ydouble, std::vector<bool> xdouble, 
00802                     SiPixelTemplate& templ, 
00803                     float& yrec1, float& yrec2, float& sigmay, float& proby,
00804                         float& xrec1, float& xrec2, float& sigmax, float& probx, int& qbin)
00805 {
00806     // Local variables 
00807         const bool deadpix = false;
00808         std::vector<std::pair<int, int> > zeropix;
00809     
00810         return SiPixelTemplateReco::PixelTempSplit(id, fpix, cotalpha, cotbeta, cluster, ydouble, xdouble, templ, 
00811                     yrec1, yrec2, sigmay, proby, xrec1, xrec2, sigmax, probx, qbin, deadpix, zeropix);
00812 
00813 } // PixelTempSplit
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/RecoLocalTracker/SiPixelRecHits/plugins/SiPixelTemplateSplit.cc
