Reference Manual

Go to the documentation of this file.
00001 //
00002 //  SiStripTemplateReco.cc Version 2.01 (V1.00 based on SiPixelTemplateReco Version 8.20)
00003 //
00004 //  V1.05 - add VI optimizations from pixel reco
00005 //  V1.07 - Improve cluster centering
00006 //  V2.00 - Change to chi2min estimator and choose barycenter when errors/bias are large
00007 //        - Increase buffer size to avoid truncation
00008 //        - Change pseudo-strip weighting to accommodate asymmetric clusters (deco mode)
00009 //        - Remove obsolete sorting of signal for weighting (truncation makes it worthless)
00010 //  V2.01 - Add barycenter bias correction
00011 //  
00012 //
00013 //
00014 //  Created by Morris Swartz on 10/27/06.
00015 //
00016 //
00017 
00018 #include <math.h>
00019 #include <algorithm>
00020 #include <vector>
00021 #include <utility>
00022 #include <iostream>
00023 // ROOT::Math has a c++ function that does the probability calc, but only in v5.12 and later
00024 #include "TMath.h"
00025 #include "Math/DistFunc.h"
00026 // Use current version of gsl instead of ROOT::Math
00027 //#include <gsl/gsl_cdf.h>
00028 
00029 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00030 #include "RecoLocalTracker/SiStripRecHitConverter/interface/SiStripTemplateReco.h"
00031 #include "RecoLocalTracker/SiStripRecHitConverter/interface/VVIObj.h"
00032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00033 #define LOGERROR(x) edm::LogError(x)
00034 #define LOGDEBUG(x) LogDebug(x)
00035 static int theVerboseLevel = 2;
00036 #define ENDL " "
00037 #include "FWCore/Utilities/interface/Exception.h"
00038 #else
00039 #include "SiStripTemplateReco.h"
00040 #include "VVIObj.h"
00041 //static int theVerboseLevel = {2};
00042 #define LOGERROR(x) std::cout << x << ": "
00043 #define LOGDEBUG(x) std::cout << x << ": "
00044 #define ENDL std::endl
00045 #endif
00046 
00047 using namespace SiStripTemplateReco;
00048 
00049 // *************************************************************************************************************************************
00083 // *************************************************************************************************************************************
00084 int SiStripTemplateReco::StripTempReco1D(int id, float cotalpha, float cotbeta, float locBy, std::vector<float>& cluster, 
00085                     SiStripTemplate& templ, 
00086                     float& xrec, float& sigmax, float& probx, int& qbin, int speed, float& probQ)
00087                         
00088 {
00089     // Local variables 
00090         int i, j, minbin, binl, binh, binq, midpix;
00091         int fxpix, nxpix, lxpix, logxpx, shiftx, ftpix, ltpix;
00092         int nclusx;
00093         int deltaj, jmin, jmax, fxbin, lxbin, djx;
00094         float sxthr, rnorm, delta, sigma, pseudopix, qscale, q50, q100;
00095         float ss2, ssa, sa2, ssba, saba, sba2, rat, fq, qtotal, barycenter, sigmaxbcn;
00096         float originx, qfx, qlx, bias, biasbcn, maxpix;
00097         double chi2x, meanx, chi2xmin, chi21max;
00098         double hchi2, hndof, prvav, mpv, sigmaQ, kappa, xvav, beta2;
00099         float xtemp[41][BSXSIZE], xsum[BSXSIZE];
00100         float chi2xbin[41], xsig2[BSXSIZE];
00101         float xw2[BSXSIZE],  xsw[BSXSIZE];
00102         bool calc_probQ, use_VVIObj;
00103         float xsize;
00104         const float probmin={1.110223e-16};
00105         const float probQmin={1.e-5};
00106         
00107 // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
00108 
00109         const double mean1pix={0.100}, chi21min={0.160};
00110                       
00111 // First, interpolate the template needed to analyze this cluster     
00112 // check to see of the track direction is in the physical range of the loaded template
00113 
00114         if(!templ.interpolate(id, cotalpha, cotbeta, locBy)) {
00115            if (theVerboseLevel > 2) {LOGDEBUG("SiStripTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha << ", cot(beta) = " << cotbeta << ", local B_y = " << locBy << ", template ID = " << id << ", no reconstruction performed" << ENDL;}        
00116            return 20;
00117         }
00118         
00119 // Check to see if Q probability is selected
00120         
00121         calc_probQ = false;
00122         use_VVIObj = false;
00123         if(speed < 0) {
00124                 calc_probQ = true;
00125                 if(speed < -1) use_VVIObj = true;
00126                 speed = 0;
00127         }
00128         
00129         if(speed > 3) {
00130                 calc_probQ = true;
00131                 if(speed < 5) use_VVIObj = true;
00132                 speed = 3;
00133         }
00134         
00135 // Get pixel dimensions from the template (to allow multiple detectors in the future)
00136         
00137         xsize = templ.xsize();
00138    
00139 // Define size of pseudopixel
00140         
00141         q50 = templ.s50();
00142         q100 = 2.f * q50;
00143         pseudopix = q50;
00144         
00145 // Get charge scaling factor
00146 
00147         qscale = templ.qscale();
00148         
00149 // enforce maximum size 
00150         
00151         nclusx = (int)cluster.size();
00152         
00153         if(nclusx > TSXSIZE) {nclusx = TSXSIZE;}
00154         
00155 // First, rescale all strip charges, sum them and trunate the strip charges      
00156 
00157         qtotal = 0.;
00158         for(i=0; i<BSXSIZE; ++i) {xsum[i] = 0.f;}
00159         maxpix = templ.sxmax();
00160         barycenter = 0.f;
00161         for(j=0; j<nclusx; ++j) {
00162                 xsum[j] = qscale*cluster[j];
00163                 qtotal += xsum[j];
00164                 barycenter += j*xsize*xsum[j];
00165            if(xsum[j] > maxpix) {xsum[j] = maxpix;}
00166    }
00167         
00168         barycenter = barycenter/qtotal - 0.5f*templ.lorxwidth();
00169                 
00170 // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx   
00171 
00172         fxpix = -1;
00173         ftpix = -1;
00174         nxpix=0;
00175         lxpix=0;
00176         ltpix=0;
00177         logxpx=0;
00178         for(i=0; i<BSXSIZE; ++i) {
00179            if(xsum[i] > 0.f) {
00180               if(fxpix == -1) {fxpix = i;}
00181                   ++logxpx;
00182                   ++nxpix;
00183                   lxpix = i;
00184                         if(xsum[i] > q100) {
00185                                 if(ftpix == -1) {ftpix = i;}
00186                                 ltpix = i;
00187          }                              
00188                 }
00189         }
00190 
00191         
00192 //      dlengthx = (float)nxpix - templ.clslenx();
00193         
00194 // Make sure cluster is continuous
00195 
00196         if((lxpix-fxpix+1) != nxpix) { 
00197         
00198            LOGDEBUG("SiStripTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
00199            if (theVerboseLevel > 2) {
00200           LOGDEBUG("SiStripTemplateReco") << "xsum[] = ";
00201           for(i=0; i<BSXSIZE-1; ++i) {LOGDEBUG("SiStripTemplateReco") << xsum[i] << ", ";}           
00202                   LOGDEBUG("SiStripTemplateReco") << ENDL;
00203        }
00204 
00205            return 2; 
00206         }
00207 
00208 // If cluster is longer than max template size, technique fails
00209 
00210         if(nxpix > TSXSIZE) { 
00211         
00212            LOGDEBUG("SiStripTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
00213            if (theVerboseLevel > 2) {
00214           LOGDEBUG("SiStripTemplateReco") << "xsum[] = ";
00215           for(i=0; i<BSXSIZE-1; ++i) {LOGDEBUG("SiStripTemplateReco") << xsum[i] << ", ";}           
00216                   LOGDEBUG("SiStripTemplateReco") << ENDL;
00217        }
00218 
00219            return 7; 
00220         }
00221         
00222 // next, center the cluster on template center if necessary   
00223 
00224         midpix = (ftpix+ltpix)/2;
00225         shiftx = templ.cxtemp() - midpix;
00226                 
00227         if(shiftx > 0) {
00228            for(i=lxpix; i>=fxpix; --i) {
00229                    xsum[i+shiftx] = xsum[i];
00230                    xsum[i] = 0.f;
00231            }
00232         } else if (shiftx < 0) {
00233            for(i=fxpix; i<=lxpix; ++i) {
00234               xsum[i+shiftx] = xsum[i];
00235                    xsum[i] = 0.f;
00236            }
00237         }
00238         lxpix +=shiftx;
00239         fxpix +=shiftx;
00240         
00241 // If the cluster boundaries are OK, add pesudopixels, otherwise quit
00242         
00243         if(fxpix > 1 && fxpix < BSXM2) {
00244            xsum[fxpix-1] = pseudopix;
00245            xsum[fxpix-2] = 0.2f*pseudopix;
00246         } else {return 9;}
00247         if(lxpix > 1 && lxpix < BSXM2) {
00248            xsum[lxpix+1] = pseudopix;
00249            xsum[lxpix+2] = 0.2f*pseudopix;
00250         } else {return 9;}
00251                         
00252 // finally, determine if pixel[0] is a double pixel and make an origin correction if it is   
00253 
00254            originx = 0.f;
00255         
00256 // uncertainty and final corrections depend upon total charge bin          
00257            
00258         fq = qtotal/templ.qavg();
00259         if(fq > 1.5f) {
00260            binq=0;
00261         } else {
00262            if(fq > 1.0f) {
00263               binq=1;
00264            } else {
00265                   if(fq > 0.85f) {
00266                          binq=2;
00267                   } else {
00268                          binq=3;
00269                   }
00270            }
00271         }
00272         
00273 // Return the charge bin via the parameter list unless the charge is too small (then flag it)
00274         
00275         qbin = binq;
00276         if(qtotal < 0.95f*templ.qmin()) {qbin = 5;} else {
00277                 if(qtotal < 0.95f*templ.qmin(1)) {qbin = 4;}
00278         }
00279         if (theVerboseLevel > 9) {
00280        LOGDEBUG("SiStripTemplateReco") <<
00281         "ID = " << id <<  
00282          " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta << 
00283          " nclusx = " << nclusx << ENDL;
00284     }
00285 
00286         
00287 // Next, copy the y- and x-templates to local arrays
00288    
00289 // First, decide on chi^2 min search parameters
00290     
00291 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00292     if(speed < 0 || speed > 3) {
00293                 throw cms::Exception("DataCorrupt") << "SiStripTemplateReco::StripTempReco2D called with illegal speed = " << speed << std::endl;
00294         }
00295 #else
00296     assert(speed >= 0 && speed < 4);
00297 #endif
00298         fxbin = 2; lxbin = 38; djx = 1;
00299     if(speed > 0) {
00300        fxbin = 8; lxbin = 32;
00301          }
00302         
00303         if(speed > 1) { 
00304            djx = 2;
00305            if(speed > 2) {
00306                   djx = 4;
00307            }
00308         }
00309         
00310         if (theVerboseLevel > 9) {
00311        LOGDEBUG("SiStripTemplateReco") <<
00312         "fxpix " << fxpix << " lxpix = " << lxpix << 
00313          " fxbin = " << fxbin << " lxbin = " << lxbin << 
00314          " djx = " << djx << " logxpx = " << logxpx << ENDL;
00315     }
00316         
00317 // Now do the copies
00318                 
00319         templ.xtemp(fxbin, lxbin, xtemp);
00320         
00321 // Do the x-reconstruction next 
00322                           
00323 // Apply the first-pass template algorithm to all clusters
00324 
00325 // Modify the template if double pixels are present 
00326                                   
00327 // Define the maximum signal to allow before de-weighting a pixel 
00328 
00329         sxthr = 1.1f*maxpix;
00330                           
00331 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
00332 
00333 //      for(i=0; i<BSXSIZE; ++i) { xsig2[i] = 0.; }
00334         templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
00335                           
00336 // Find the template bin that minimizes the Chi^2 
00337 
00338         chi2xmin = 1.e15;
00339         for(i=fxbin; i<=lxbin; ++i) { chi2xbin[i] = -1.e15f;}
00340         ss2 = 0.f;
00341         for(i=fxpix-2; i<=lxpix+2; ++i) {
00342                 xw2[i] = 1.f/xsig2[i];
00343                 xsw[i] = xsum[i]*xw2[i];
00344                 ss2 += xsum[i]*xsw[i];
00345         }
00346         minbin = -1;
00347         deltaj = djx;
00348         jmin = fxbin;
00349         jmax = lxbin;
00350         while(deltaj > 0) {
00351            for(j=jmin; j<=jmax; j+=deltaj) {
00352               if(chi2xbin[j] < -100.f) {
00353                                 ssa = 0.f;
00354                                 sa2 = 0.f;
00355                                 for(i=fxpix-2; i<=lxpix+2; ++i) {
00356                                         ssa += xsw[i]*xtemp[j][i];
00357                                         sa2 += xtemp[j][i]*xtemp[j][i]*xw2[i];
00358                                 }
00359                      rat=ssa/ss2;
00360                      if(rat <= 0.f) {LOGERROR("SiStripTemplateReco") << "illegal chi2xmin normalization (1) = " << rat << ENDL; rat = 1.;}
00361                      chi2xbin[j]=ss2-2.*ssa/rat+sa2/(rat*rat);
00362                   }
00363                   if(chi2xbin[j] < chi2xmin) {
00364                           chi2xmin = chi2xbin[j];
00365                           minbin = j;
00366                   }
00367            } 
00368            deltaj /= 2;
00369            if(minbin > fxbin) {jmin = minbin - deltaj;} else {jmin = fxbin;}
00370            if(minbin < lxbin) {jmax = minbin + deltaj;} else {jmax = lxbin;}
00371         }
00372         
00373         if (theVerboseLevel > 9) {
00374        LOGDEBUG("SiStripTemplateReco") <<
00375         "minbin " << minbin << " chi2xmin = " << chi2xmin << ENDL;
00376     }
00377 
00378 // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
00379         
00380         if(nxpix == 1) {
00381         
00382                 delta = templ.dxone();
00383                 sigma = templ.sxone();
00384            xrec = 0.5f*(fxpix+lxpix-2*shiftx+2.f*originx)*xsize-delta;
00385            if(sigma <= 0.f) {
00386               sigmax = 28.9f;
00387            } else {
00388           sigmax = sigma;
00389            }
00390            
00391 // Do probability calculation for one-pixel clusters
00392 
00393                 chi21max = fmax(chi21min, (double)templ.chi2xminone());
00394                 chi2xmin -=chi21max;
00395                 if(chi2xmin < 0.) {chi2xmin = 0.;}
00396                 meanx = fmax(mean1pix, (double)templ.chi2xavgone());
00397                 hchi2 = chi2xmin/2.; hndof = meanx/2.;
00398                 probx = 1. - TMath::Gamma(hndof, hchi2);
00399            
00400         } else {
00401                 
00402 // Now make the second, interpolating pass with the templates 
00403         
00404                 binl = minbin - 1;
00405         binh = binl + 2;
00406         if(binl < fxbin) { binl = fxbin;}
00407         if(binh > lxbin) { binh = lxbin;}         
00408         ssa = 0.;
00409         sa2 = 0.;
00410         ssba = 0.;
00411         saba = 0.;
00412         sba2 = 0.;
00413         for(i=fxpix-2; i<=lxpix+2; ++i) {
00414             ssa += xsw[i]*xtemp[binl][i];
00415             sa2 += xtemp[binl][i]*xtemp[binl][i]*xw2[i];
00416             ssba += xsw[i]*(xtemp[binh][i] - xtemp[binl][i]);
00417                         saba += xtemp[binl][i]*(xtemp[binh][i] - xtemp[binl][i])*xw2[i];
00418                         sba2 += (xtemp[binh][i] - xtemp[binl][i])*(xtemp[binh][i] - xtemp[binl][i])*xw2[i];
00419         }
00420         
00421 // rat is the fraction of the "distance" from template a to template b     
00422         
00423         rat=(ssba*ssa-ss2*saba)/(ss2*sba2-ssba*ssba);
00424         if(rat < 0.f) {rat=0.f;}
00425         if(rat > 1.f) {rat=1.0f;}
00426         rnorm = (ssa+rat*ssba)/ss2;
00427         
00428 // Calculate the charges in the first and last pixels
00429         
00430         qfx = xsum[fxpix];
00431         if(logxpx > 1) {
00432             qlx=xsum[lxpix];
00433             } else {
00434             qlx = qfx;
00435             }
00436                 
00437 //  Now calculate the mean bias correction and uncertainties
00438         
00439         float qxfrac = (qfx-qlx)/(qfx+qlx);
00440                 bias = templ.xflcorr(binq,qxfrac)+templ.xavg(binq);
00441         
00442 // uncertainty and final correction depend upon charge bin         
00443         
00444         xrec = (0.125f*binl+BSHX-2.5f+rat*(binh-binl)*0.125f-(float)shiftx+originx)*xsize - bias;
00445         sigmax = templ.xrms(binq);
00446         
00447 // Do goodness of fit test in x  
00448         
00449         if(rnorm <= 0.f) {LOGERROR("SiStripTemplateReco") << "illegal chi2x normalization (2) = " << rnorm << ENDL; rnorm = 1.;}
00450         chi2x=ss2-2.f/rnorm*ssa-2.f/rnorm*rat*ssba+(sa2+2.f*rat*saba+rat*rat*sba2)/(rnorm*rnorm)-templ.chi2xmin(binq);
00451         if(chi2x < 0.0) {chi2x = 0.0;}
00452         meanx = templ.chi2xavg(binq);
00453         if(meanx < 0.01) {meanx = 0.01;}
00454         // gsl function that calculates the chi^2 tail prob for non-integral dof
00455         //         probx = gsl_cdf_chisq_Q(chi2x, meanx);
00456         //         probx = ROOT::Math::chisquared_cdf_c(chi2x, meanx, trx0);
00457         hchi2 = chi2x/2.; hndof = meanx/2.;
00458         probx = 1. - TMath::Gamma(hndof, hchi2);
00459                            
00460 // Now choose the better result
00461 
00462         bias = templ.xavg(binq);
00463             biasbcn = templ.xavgbcn(binq);
00464             sigmaxbcn = templ.xrmsbcn(binq);
00465                 
00466                 if((bias*bias+sigmax*sigmax) > (biasbcn*biasbcn+sigmaxbcn*sigmaxbcn)) {
00467                                 
00468                         xrec = barycenter - biasbcn;
00469                         sigmax = sigmaxbcn;
00470             
00471                 }       
00472            
00473         }
00474         
00475 //  Don't return exact zeros for the probability
00476         
00477         if(probx < probmin) {probx = probmin;}
00478         
00479 //  Decide whether to generate a cluster charge probability
00480         
00481         if(calc_probQ) {
00482                 
00483 // Calculate the Vavilov probability that the cluster charge is OK
00484         
00485                 templ.vavilov_pars(mpv, sigmaQ, kappa);
00486 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00487                 if((sigmaQ <=0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
00488                         throw cms::Exception("DataCorrupt") << "SiStripTemplateReco::Vavilov parameters mpv/sigmaQ/kappa = " << mpv << "/" << sigmaQ << "/" << kappa << std::endl;
00489                 }
00490 #else
00491                 assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
00492 #endif
00493                 xvav = ((double)qtotal-mpv)/sigmaQ;
00494                 beta2 = 1.;
00495                 if(use_VVIObj) {                        
00496 //  VVIObj is a private port of CERNLIB VVIDIS
00497                    VVIObj vvidist(kappa, beta2, 1);
00498                    prvav = vvidist.fcn(xvav);                   
00499                 } else {
00500 //  Use faster but less accurate TMath Vavilov distribution function
00501                         prvav = TMath::VavilovI(xvav, kappa, beta2);
00502       }
00503 //  Change to upper tail probability
00504 //              if(prvav > 0.5) prvav = 1. - prvav;
00505 //              probQ = (float)(2.*prvav);
00506                 probQ = 1. - prvav;
00507                 if(probQ < probQmin) {probQ = probQmin;}
00508         } else {
00509                 probQ = -1;
00510         }
00511         
00512         return 0;
00513 } // StripTempReco2D 
00514
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoLocalTracker/SiStripRecHitConverter/src/SiStripTemplateReco.cc
