Reference Manual

Go to the documentation of this file.
00001 //
00002 //  SiStripTemplateSplit.cc
00003 //
00004 //  Procedure to fit two templates (same angle hypotheses) to a single cluster
00005 //
00006 //  Version 1.00 [based on SiPixelTemplateSplit.cc Version 2.30]
00007 //  Version 1.01 [improve error estimation for qbin=3 events]
00008 //  Version 1.05 [Incorporate VI-like speed improvements]
00009 //  Version 1.06 Clean-up irrelevant (since truncation) sorting
00010 //  Version 2.10 Clone speed improvements from the pixels (eliminate 3-d multi-arays, improve seach algorithm)
00011 //
00012 //  Created by Morris Swartz on 04/10/08.
00013 //
00014 //
00015 
00016 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00017 //#include <cmath.h>
00018 #else
00019 #include <math.h>
00020 #endif
00021 #include <algorithm>
00022 #include <vector>
00023 #include <list>
00024 #include <iostream>
00025 // ROOT::Math has a c++ function that does the probability calc, but only in v5.12 and later
00026 #include "Math/DistFunc.h"
00027 #include "TMath.h"
00028 // Use current version of gsl instead of ROOT::Math
00029 //#include <gsl/gsl_cdf.h>
00030 
00031 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00032 #include "RecoLocalTracker/SiStripRecHitConverter/interface/SiStripTemplateSplit.h"
00033 #include "RecoLocalTracker/SiStripRecHitConverter/interface/VVIObj.h"
00034 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00035 #define LOGERROR(x) edm::LogError(x)
00036 #define LOGDEBUG(x) LogDebug(x)
00037 static int theVerboseLevel = 2;
00038 #define ENDL " "
00039 #else
00040 #include "SiStripTemplateSplit.h"
00041 #include "VVIObj.h"
00042 //#include "SiStripTemplate2D.h"
00043 //#include "SimpleTemplate2D.h"
00044 //static int theVerboseLevel = {2};
00045 #define LOGERROR(x) std::cout << x << ": "
00046 #define LOGDEBUG(x) std::cout << x << ": "
00047 #define ENDL std::endl
00048 #endif
00049 
00050 using namespace SiStripTemplateSplit;
00051 
00052 // *************************************************************************************************************************************
00074 // *************************************************************************************************************************************
00075 int SiStripTemplateSplit::StripTempSplit(int id, float cotalpha, float cotbeta, float locBy, int speed, std::vector<float>& cluster, 
00076                     SiStripTemplate& templ, 
00077                         float& xrec1, float& xrec2, float& sigmax, float& prob2x, int& q2bin, float& prob2Q)
00078                         
00079 {
00080     // Local variables 
00081         int i, j, k, binq, binqerr, midpix;
00082         int fxpix, nxpix, lxpix, logxpx, shiftx;
00083         int nclusx;
00084         int nxbin, xcbin, minbinj, minbink;
00085         int deltaj, jmin, jmax, kmin, kmax, km, fxbin, lxbin, djx;
00086         float sxthr, delta, sigma, pseudopix, xsize, qscale;
00087         float ss2, ssa, sa2, rat, fq, qtotal, qavg;
00088         float originx, bias, maxpix;
00089         double chi2x, meanx, chi2xmin, chi21max;
00090         double hchi2, hndof;
00091         double prvav, mpv, sigmaQ, kappa, xvav, beta2;
00092         float xsum[BSXSIZE];
00093         float xsig2[BSXSIZE];
00094         float xw2[BSXSIZE], xsw[BSXSIZE];
00095         const float sqrt2x={2.00000};
00096         const float sqrt12={3.4641};
00097         const float probmin={1.110223e-16};
00098         const float prob2Qmin={1.e-5};
00099         
00100 //      bool SiStripTemplateSplit::SimpleTemplate2D(float cotalpha, float cotbeta, float xhit, float yhit, float thick, float lorxwidth, float lorywidth, 
00101 //                                                float qavg, std::vector<bool> ydouble, std::vector<bool> xdouble, float template2d[BSXM2][BYM2]);
00102         
00103 // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
00104 
00105         const double mean1pix={0.100}, chi21min={0.160};
00106                       
00107 // First, interpolate the template needed to analyze this cluster     
00108 // check to see of the track direction is in the physical range of the loaded template
00109 
00110         if(!templ.interpolate(id, cotalpha, cotbeta, locBy)) {
00111            LOGDEBUG("SiStripTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha << ", cot(beta) = " << cotbeta 
00112            << " is not within the acceptance of template ID = " << id << ", no reconstruction performed" << ENDL;       
00113            return 20;
00114         }
00115                       
00116         // Get pixel dimensions from the template (to allow multiple detectors in the future)
00117         
00118         xsize = templ.xsize();
00119    
00120         // Define size of pseudopixel
00121         
00122         pseudopix = templ.s50();
00123         
00124         // Get charge scaling factor
00125         
00126         qscale = templ.qscale();
00127         
00128         // enforce maximum size 
00129         
00130         nclusx = (int)cluster.size();
00131         
00132         if(nclusx > TSXSIZE) {nclusx = TSXSIZE;}
00133         
00134         // First, rescale all strip charges, sum them and trunate the strip charges      
00135         
00136         qtotal = 0.f;
00137         for(i=0; i<BSXSIZE; ++i) {xsum[i] = 0.f;}
00138         maxpix = 2.f*templ.sxmax();
00139         for(j=0; j<nclusx; ++j) {
00140                 xsum[j] = qscale*cluster[j];
00141                 qtotal += xsum[j];
00142            if(xsum[j] > maxpix) {xsum[j] = maxpix;}
00143    }
00144         
00145         // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx   
00146         
00147         fxpix = -1;
00148         nxpix=0;
00149         lxpix=0;
00150         logxpx=0;
00151         for(i=0; i<BSXSIZE; ++i) {
00152            if(xsum[i] > 0.f) {
00153               if(fxpix == -1) {fxpix = i;}
00154                         ++logxpx;
00155                         ++nxpix;
00156                         lxpix = i;
00157                 }
00158         }
00159         
00160         
00161         //      dlengthx = (float)nxpix - templ.clslenx();
00162         
00163         // Make sure cluster is continuous
00164         
00165         if((lxpix-fxpix+1) != nxpix) { 
00166                 
00167            LOGDEBUG("SiStripTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
00168            if (theVerboseLevel > 2) {
00169                         LOGDEBUG("SiStripTemplateReco") << "xsum[] = ";
00170                         for(i=0; i<BSXSIZE-1; ++i) {LOGDEBUG("SiStripTemplateReco") << xsum[i] << ", ";}           
00171                         LOGDEBUG("SiStripTemplateReco") << ENDL;
00172                 }
00173                 
00174            return 2; 
00175         }
00176         
00177         // If cluster is longer than max template size, technique fails
00178         
00179         if(nxpix > TSXSIZE) { 
00180                 
00181            LOGDEBUG("SiStripTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
00182            if (theVerboseLevel > 2) {
00183                         LOGDEBUG("SiStripTemplateReco") << "xsum[] = ";
00184                         for(i=0; i<BSXSIZE-1; ++i) {LOGDEBUG("SiStripTemplateReco") << xsum[i] << ", ";}           
00185                         LOGDEBUG("SiStripTemplateReco") << ENDL;
00186                 }
00187                 
00188            return 7; 
00189         }
00190         
00191         // next, center the cluster on template center if necessary   
00192         
00193         midpix = (fxpix+lxpix)/2;
00194         shiftx = templ.cxtemp() - midpix;
00195         if(shiftx > 0) {
00196            for(i=lxpix; i>=fxpix; --i) {
00197                    xsum[i+shiftx] = xsum[i];
00198                    xsum[i] = 0.f;
00199            }
00200         } else if (shiftx < 0) {
00201            for(i=fxpix; i<=lxpix; ++i) {
00202               xsum[i+shiftx] = xsum[i];
00203                    xsum[i] = 0.f;
00204            }
00205         }
00206         lxpix +=shiftx;
00207         fxpix +=shiftx;
00208         
00209         // If the cluster boundaries are OK, add pesudopixels, otherwise quit
00210         
00211         if(fxpix > 1 && fxpix <BSXM2) {
00212            xsum[fxpix-1] = pseudopix;
00213            xsum[fxpix-2] = 0.2f*pseudopix;
00214         } else {return 9;}
00215         if(lxpix > 1 && lxpix < BSXM2) {
00216            xsum[lxpix+1] = pseudopix;
00217            xsum[lxpix+2] = 0.2f*pseudopix;
00218         } else {return 9;}
00219         
00220         // finally, determine if pixel[0] is a double pixel and make an origin correction if it is   
00221         
00222         originx = 0.f;
00223                 
00224 // uncertainty and final corrections depend upon total charge bin          
00225            
00226         qavg = templ.qavg();
00227         fq = qtotal/qavg;
00228         if(fq > 3.0f) {
00229            binq=0;
00230         } else {
00231            if(fq > 2.0f) {
00232               binq=1;
00233            } else {
00234                   if(fq > 1.70f) {
00235                          binq=2;
00236                   } else {
00237                          binq=3;
00238                   }
00239            }
00240         }
00241         
00242         // Return the charge bin via the parameter list unless the charge is too small (then flag it)
00243         
00244         q2bin = binq;
00245         if(qtotal < 1.9f*templ.qmin()) {q2bin = 5;} else {
00246                 if(qtotal < 1.9f*templ.qmin(1)) {q2bin = 4;}
00247         }
00248         if (theVerboseLevel > 9) {
00249                 LOGDEBUG("SiStripTemplateReco") <<
00250                 "ID = " << id <<  
00251                 " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta << 
00252                 " nclusx = " << nclusx << ENDL;
00253         }
00254         
00255 // binqerr is the charge bin for error estimation
00256         
00257         binqerr = binq;
00258         if(binqerr > 2) binqerr = 2;    
00259                 
00260 // Calculate the Vavilov probability that the cluster charge is consistent with a merged cluster
00261                 
00262         if(speed < 0) {
00263            templ.vavilov2_pars(mpv, sigmaQ, kappa);
00264 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00265            if((sigmaQ <=0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
00266                         throw cms::Exception("DataCorrupt") << "SiStripTemplateSplit::Vavilov parameters mpv/sigmaQ/kappa = " << mpv << "/" << sigmaQ << "/" << kappa << std::endl;
00267            }
00268 #else
00269            assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
00270 #endif
00271            xvav = ((double)qtotal-mpv)/sigmaQ;
00272            beta2 = 1.;
00273 //  VVIObj is a private port of CERNLIB VVIDIS
00274            VVIObj vvidist(kappa, beta2, 1);
00275            prvav = vvidist.fcn(xvav);                   
00276            prob2Q = 1. - prvav;
00277            if(prob2Q < prob2Qmin) {prob2Q = prob2Qmin;}
00278         } else {
00279                 prob2Q = -1.f;
00280         }
00281         
00282         
00283 // Next, generate the 3d x-template
00284    
00285         templ.xtemp3d_int(nxpix, nxbin);
00286         
00287 // retrieve the number of x-bins 
00288         
00289         xcbin = nxbin/2;
00290         
00291 // Next, decide on chi^2 min search parameters
00292         
00293 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00294         if(speed < -1 || speed > 2) {
00295                 throw cms::Exception("DataCorrupt") << "SiStripTemplateReco::PixelTempReco2D called with illegal speed = " << speed << std::endl;
00296         }
00297 #else
00298         assert(speed >= -1 && speed < 3);
00299 #endif
00300         fxbin = 0; lxbin = nxbin-1; djx = 1;
00301         if(speed > 0) {
00302         djx = 2;
00303         if(speed > 1) {djx = 4;}
00304         }
00305                                                           
00306 // Define the maximum signal to allow before de-weighting a pixel 
00307 
00308         sxthr = 1.1f*maxpix;
00309                                    
00310 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
00311 
00312 //      for(i=0; i<BYSIZE; ++i) { xsig2[i] = 0.; }
00313         templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
00314                           
00315 // Find the template bin that minimizes the Chi^2 
00316 
00317         chi2xmin = 1.e15;
00318         ss2 = 0.f;
00319         for(i=fxpix-2; i<=lxpix+2; ++i) {
00320                 xw2[i] = 1.f/xsig2[i];
00321                 xsw[i] = xsum[i]*xw2[i];
00322                 ss2 += xsum[i]*xsw[i];
00323         }
00324         minbinj = -1; 
00325         minbink = -1;
00326         deltaj = djx;
00327         jmin = fxbin;
00328         jmax = lxbin;
00329         kmin = fxbin;
00330         kmax = lxbin;
00331     std::vector<float> xtemp(BSXSIZE);
00332         while(deltaj > 0) {
00333         for(j=jmin; j<jmax; j+=deltaj) {
00334                         km = std::min(kmax, j);
00335             for(k=kmin; k<=km; k+=deltaj) {
00336                 
00337                 // Get the template for this set of indices
00338                 
00339                 templ.xtemp3d(j, k, xtemp);
00340                 
00341                 ssa = 0.f;
00342                 sa2 = 0.f;
00343                 for(i=fxpix-2; i<=lxpix+2; ++i) {
00344                     ssa += xsw[i]*xtemp[i];
00345                     sa2 += xtemp[i]*xtemp[i]*xw2[i];
00346                 }
00347                 rat=ssa/ss2;
00348                 if(rat <= 0.f) {LOGERROR("SiStripTemplateSplit") << "illegal chi2xmin normalization = " << rat << ENDL; rat = 1.;}
00349                 chi2x=ss2-2.f*ssa/rat+sa2/(rat*rat);
00350                 if(chi2x < chi2xmin) {
00351                     chi2xmin = chi2x;
00352                     minbinj = j;
00353                     minbink = k;
00354                                 }
00355             }
00356         } 
00357         deltaj /= 2;
00358         if(minbinj > fxbin) {jmin = minbinj - deltaj;} else {jmin = fxbin;}
00359         if(minbinj < lxbin) {jmax = minbinj + deltaj;} else {jmax = lxbin;}
00360         if(minbink > fxbin) {kmin = minbink - deltaj;} else {kmin = fxbin;}
00361         if(minbink < lxbin) {kmax = minbink + deltaj;} else {kmax = lxbin;}             
00362         }
00363         
00364         if (theVerboseLevel > 9) {
00365                 LOGDEBUG("SiStripTemplateReco") <<
00366                 "minbinj/minbink " << minbinj<< "/" << minbink << " chi2xmin = " << chi2xmin << ENDL;
00367         }
00368 
00369 // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
00370         
00371         if(logxpx == 1) {
00372         
00373                 delta = templ.dxone();
00374                 sigma = templ.sxone();
00375                 xrec1 = 0.5f*(fxpix+lxpix-2*shiftx+2.f*originx)*xsize-delta;
00376            xrec2 = xrec1;
00377            if(sigma <= 0.f) {
00378               sigmax = xsize/sqrt12;
00379            } else {
00380           sigmax = sigma;
00381            }
00382            
00383 // Do probability calculation for one-pixel clusters
00384 
00385                 chi21max = fmax(chi21min, (double)templ.chi2xminone());
00386                 chi2xmin -=chi21max;
00387                 if(chi2xmin < 0.) {chi2xmin = 0.;}
00388                 meanx = fmax(mean1pix, (double)templ.chi2xavgone());
00389                 hchi2 = chi2xmin/2.; hndof = meanx/2.;
00390                 prob2x = 1. - TMath::Gamma(hndof, hchi2);
00391            
00392         } else {
00393            
00394                         
00395 // uncertainty and final correction depend upon charge bin         
00396            
00397                 bias = templ.xavgc2m(binq);
00398                 k = std::min(minbink, minbinj);
00399                 j = std::max(minbink, minbinj);
00400                 xrec1 = (0.125f*(minbink-xcbin)+BSHX-(float)shiftx+originx)*xsize - bias;
00401                 xrec2 = (0.125f*(minbinj-xcbin)+BSHX-(float)shiftx+originx)*xsize - bias;
00402                 sigmax = sqrt2x*templ.xrmsc2m(binqerr);
00403                         
00404            
00405 // Do goodness of fit test in y  
00406            
00407            if(chi2xmin < 0.0) {chi2xmin = 0.0;}
00408            meanx = templ.chi2xavgc2m(binq);
00409            if(meanx < 0.01) {meanx = 0.01;}
00410        hchi2 = chi2xmin/2.; hndof = meanx/2.;
00411            prob2x = 1. - TMath::Gamma(hndof, hchi2);
00412    }
00413         
00414         //  Don't return exact zeros for the probability
00415         
00416         if(prob2x < probmin) {prob2x = probmin;}
00417         
00418         
00419     return 0;
00420 } // StripTempSplit 
00421 
00422 
00423 // *************************************************************************************************************************************
00439 // *************************************************************************************************************************************
00440 int SiStripTemplateSplit::StripTempSplit(int id, float cotalpha, float cotbeta, float locBy, std::vector<float>& cluster, 
00441                                          SiStripTemplate& templ, 
00442                                          float& xrec1, float& xrec2, float& sigmax, float& prob2x, int& q2bin, float& prob2Q)
00443 
00444 {
00445     // Local variables 
00446         const int speed = 1;
00447 
00448     return SiStripTemplateSplit::StripTempSplit(id, cotalpha, cotbeta, locBy, speed, cluster, templ, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q);
00449 } // StripTempSplit
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoLocalTracker/SiStripRecHitConverter/src/SiStripTemplateSplit.cc
