/data/refman/pasoursint/CMSSW_5_3_9_patch3/src/RecoLocalTracker/SiStripRecHitConverter/src/SiStripTemplateReco.cc

Go to the documentation of this file.

//
//  SiStripTemplateReco.cc Version 2.01 (V1.00 based on SiPixelTemplateReco Version 8.20)
//
//  V1.05 - add VI optimizations from pixel reco
//  V1.07 - Improve cluster centering
//  V2.00 - Change to chi2min estimator and choose barycenter when errors/bias are large
//        - Increase buffer size to avoid truncation
//        - Change pseudo-strip weighting to accommodate asymmetric clusters (deco mode)
//        - Remove obsolete sorting of signal for weighting (truncation makes it worthless)
//  V2.01 - Add barycenter bias correction
//  
//
//
//  Created by Morris Swartz on 10/27/06.
//
//

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

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

using namespace SiStripTemplateReco;

// *************************************************************************************************************************************
// *************************************************************************************************************************************
int SiStripTemplateReco::StripTempReco1D(int id, float cotalpha, float cotbeta, float locBy, std::vector<float>& cluster, 
                    SiStripTemplate& templ, 
                    float& xrec, float& sigmax, float& probx, int& qbin, int speed, float& probQ)
                        
{
    // Local variables 
        int i, j, minbin, binl, binh, binq, midpix;
        int fxpix, nxpix, lxpix, logxpx, shiftx, ftpix, ltpix;
        int nclusx;
        int deltaj, jmin, jmax, fxbin, lxbin, djx;
        float sxthr, rnorm, delta, sigma, pseudopix, qscale, q50, q100;
        float ss2, ssa, sa2, ssba, saba, sba2, rat, fq, qtotal, barycenter, sigmaxbcn;
        float originx, qfx, qlx, bias, biasbcn, maxpix;
        double chi2x, meanx, chi2xmin, chi21max;
        double hchi2, hndof, prvav, mpv, sigmaQ, kappa, xvav, beta2;
        float xtemp[41][BSXSIZE], xsum[BSXSIZE];
        float chi2xbin[41], xsig2[BSXSIZE];
        float xw2[BSXSIZE],  xsw[BSXSIZE];
        bool calc_probQ, use_VVIObj;
        float xsize;
        const float probmin={1.110223e-16};
        const float probQmin={1.e-5};
        
// The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only

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

        if(!templ.interpolate(id, cotalpha, cotbeta, locBy)) {
           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;}        
           return 20;
        }
        
// Check to see if Q probability is selected
        
        calc_probQ = false;
        use_VVIObj = false;
        if(speed < 0) {
                calc_probQ = true;
                if(speed < -1) use_VVIObj = true;
                speed = 0;
        }
        
        if(speed > 3) {
                calc_probQ = true;
                if(speed < 5) use_VVIObj = true;
                speed = 3;
        }
        
// Get pixel dimensions from the template (to allow multiple detectors in the future)
        
        xsize = templ.xsize();
   
// Define size of pseudopixel
        
        q50 = templ.s50();
        q100 = 2.f * q50;
        pseudopix = q50;
        
// Get charge scaling factor

        qscale = templ.qscale();
        
// enforce maximum size 
        
        nclusx = (int)cluster.size();
        
        if(nclusx > TSXSIZE) {nclusx = TSXSIZE;}
        
// First, rescale all strip charges, sum them and trunate the strip charges      

        qtotal = 0.;
        for(i=0; i<BSXSIZE; ++i) {xsum[i] = 0.f;}
        maxpix = templ.sxmax();
        barycenter = 0.f;
        for(j=0; j<nclusx; ++j) {
                xsum[j] = qscale*cluster[j];
                qtotal += xsum[j];
                barycenter += j*xsize*xsum[j];
           if(xsum[j] > maxpix) {xsum[j] = maxpix;}
   }
        
        barycenter = barycenter/qtotal - 0.5f*templ.lorxwidth();
                
// next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx   

        fxpix = -1;
        ftpix = -1;
        nxpix=0;
        lxpix=0;
        ltpix=0;
        logxpx=0;
        for(i=0; i<BSXSIZE; ++i) {
           if(xsum[i] > 0.f) {
              if(fxpix == -1) {fxpix = i;}
                  ++logxpx;
                  ++nxpix;
                  lxpix = i;
                        if(xsum[i] > q100) {
                                if(ftpix == -1) {ftpix = i;}
                                ltpix = i;
         }                              
                }
        }

        
//      dlengthx = (float)nxpix - templ.clslenx();
        
// Make sure cluster is continuous

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

           return 2; 
        }

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

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

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

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

           originx = 0.f;
        
// uncertainty and final corrections depend upon total charge bin          
           
        fq = qtotal/templ.qavg();
        if(fq > 1.5f) {
           binq=0;
        } else {
           if(fq > 1.0f) {
              binq=1;
           } else {
                  if(fq > 0.85f) {
                         binq=2;
                  } else {
                         binq=3;
                  }
           }
        }
        
// Return the charge bin via the parameter list unless the charge is too small (then flag it)
        
        qbin = binq;
        if(qtotal < 0.95f*templ.qmin()) {qbin = 5;} else {
                if(qtotal < 0.95f*templ.qmin(1)) {qbin = 4;}
        }
        if (theVerboseLevel > 9) {
       LOGDEBUG("SiStripTemplateReco") <<
        "ID = " << id <<  
         " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta << 
         " nclusx = " << nclusx << ENDL;
    }

        
// Next, copy the y- and x-templates to local arrays
   
// First, decide on chi^2 min search parameters
    
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
    if(speed < 0 || speed > 3) {
                throw cms::Exception("DataCorrupt") << "SiStripTemplateReco::StripTempReco2D called with illegal speed = " << speed << std::endl;
        }
#else
    assert(speed >= 0 && speed < 4);
#endif
        fxbin = 2; lxbin = 38; djx = 1;
    if(speed > 0) {
       fxbin = 8; lxbin = 32;
         }
        
        if(speed > 1) { 
           djx = 2;
           if(speed > 2) {
                  djx = 4;
           }
        }
        
        if (theVerboseLevel > 9) {
       LOGDEBUG("SiStripTemplateReco") <<
        "fxpix " << fxpix << " lxpix = " << lxpix << 
         " fxbin = " << fxbin << " lxbin = " << lxbin << 
         " djx = " << djx << " logxpx = " << logxpx << ENDL;
    }
        
// Now do the copies
                
        templ.xtemp(fxbin, lxbin, xtemp);
        
// Do the x-reconstruction next 
                          
// Apply the first-pass template algorithm to all clusters

// Modify the template if double pixels are present 
                                  
// Define the maximum signal to allow before de-weighting a pixel 

        sxthr = 1.1f*maxpix;
                          
// Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 

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

        chi2xmin = 1.e15;
        for(i=fxbin; i<=lxbin; ++i) { chi2xbin[i] = -1.e15f;}
        ss2 = 0.f;
        for(i=fxpix-2; i<=lxpix+2; ++i) {
                xw2[i] = 1.f/xsig2[i];
                xsw[i] = xsum[i]*xw2[i];
                ss2 += xsum[i]*xsw[i];
        }
        minbin = -1;
        deltaj = djx;
        jmin = fxbin;
        jmax = lxbin;
        while(deltaj > 0) {
           for(j=jmin; j<=jmax; j+=deltaj) {
              if(chi2xbin[j] < -100.f) {
                                ssa = 0.f;
                                sa2 = 0.f;
                                for(i=fxpix-2; i<=lxpix+2; ++i) {
                                        ssa += xsw[i]*xtemp[j][i];
                                        sa2 += xtemp[j][i]*xtemp[j][i]*xw2[i];
                                }
                     rat=ssa/ss2;
                     if(rat <= 0.f) {LOGERROR("SiStripTemplateReco") << "illegal chi2xmin normalization (1) = " << rat << ENDL; rat = 1.;}
                     chi2xbin[j]=ss2-2.*ssa/rat+sa2/(rat*rat);
                  }
                  if(chi2xbin[j] < chi2xmin) {
                          chi2xmin = chi2xbin[j];
                          minbin = j;
                  }
           } 
           deltaj /= 2;
           if(minbin > fxbin) {jmin = minbin - deltaj;} else {jmin = fxbin;}
           if(minbin < lxbin) {jmax = minbin + deltaj;} else {jmax = lxbin;}
        }
        
        if (theVerboseLevel > 9) {
       LOGDEBUG("SiStripTemplateReco") <<
        "minbin " << minbin << " chi2xmin = " << chi2xmin << ENDL;
    }

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

                chi21max = fmax(chi21min, (double)templ.chi2xminone());
                chi2xmin -=chi21max;
                if(chi2xmin < 0.) {chi2xmin = 0.;}
                meanx = fmax(mean1pix, (double)templ.chi2xavgone());
                hchi2 = chi2xmin/2.; hndof = meanx/2.;
                probx = 1. - TMath::Gamma(hndof, hchi2);
           
        } else {
                
// Now make the second, interpolating pass with the templates 
        
                binl = minbin - 1;
        binh = binl + 2;
        if(binl < fxbin) { binl = fxbin;}
        if(binh > lxbin) { binh = lxbin;}         
        ssa = 0.;
        sa2 = 0.;
        ssba = 0.;
        saba = 0.;
        sba2 = 0.;
        for(i=fxpix-2; i<=lxpix+2; ++i) {
            ssa += xsw[i]*xtemp[binl][i];
            sa2 += xtemp[binl][i]*xtemp[binl][i]*xw2[i];
            ssba += xsw[i]*(xtemp[binh][i] - xtemp[binl][i]);
                        saba += xtemp[binl][i]*(xtemp[binh][i] - xtemp[binl][i])*xw2[i];
                        sba2 += (xtemp[binh][i] - xtemp[binl][i])*(xtemp[binh][i] - xtemp[binl][i])*xw2[i];
        }
        
// rat is the fraction of the "distance" from template a to template b     
        
        rat=(ssba*ssa-ss2*saba)/(ss2*sba2-ssba*ssba);
        if(rat < 0.f) {rat=0.f;}
        if(rat > 1.f) {rat=1.0f;}
        rnorm = (ssa+rat*ssba)/ss2;
        
// Calculate the charges in the first and last pixels
        
        qfx = xsum[fxpix];
        if(logxpx > 1) {
            qlx=xsum[lxpix];
            } else {
            qlx = qfx;
            }
                
//  Now calculate the mean bias correction and uncertainties
        
        float qxfrac = (qfx-qlx)/(qfx+qlx);
                bias = templ.xflcorr(binq,qxfrac)+templ.xavg(binq);
        
// uncertainty and final correction depend upon charge bin         
        
        xrec = (0.125f*binl+BSHX-2.5f+rat*(binh-binl)*0.125f-(float)shiftx+originx)*xsize - bias;
        sigmax = templ.xrms(binq);
        
// Do goodness of fit test in x  
        
        if(rnorm <= 0.f) {LOGERROR("SiStripTemplateReco") << "illegal chi2x normalization (2) = " << rnorm << ENDL; rnorm = 1.;}
        chi2x=ss2-2.f/rnorm*ssa-2.f/rnorm*rat*ssba+(sa2+2.f*rat*saba+rat*rat*sba2)/(rnorm*rnorm)-templ.chi2xmin(binq);
        if(chi2x < 0.0) {chi2x = 0.0;}
        meanx = templ.chi2xavg(binq);
        if(meanx < 0.01) {meanx = 0.01;}
        // gsl function that calculates the chi^2 tail prob for non-integral dof
        //         probx = gsl_cdf_chisq_Q(chi2x, meanx);
        //         probx = ROOT::Math::chisquared_cdf_c(chi2x, meanx, trx0);
        hchi2 = chi2x/2.; hndof = meanx/2.;
        probx = 1. - TMath::Gamma(hndof, hchi2);
                           
// Now choose the better result

        bias = templ.xavg(binq);
            biasbcn = templ.xavgbcn(binq);
            sigmaxbcn = templ.xrmsbcn(binq);
                
                if((bias*bias+sigmax*sigmax) > (biasbcn*biasbcn+sigmaxbcn*sigmaxbcn)) {
                                
                        xrec = barycenter - biasbcn;
                        sigmax = sigmaxbcn;
            
                }       
           
        }
        
//  Don't return exact zeros for the probability
        
        if(probx < probmin) {probx = probmin;}
        
//  Decide whether to generate a cluster charge probability
        
        if(calc_probQ) {
                
// Calculate the Vavilov probability that the cluster charge is OK
        
                templ.vavilov_pars(mpv, sigmaQ, kappa);
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
                if((sigmaQ <=0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
                        throw cms::Exception("DataCorrupt") << "SiStripTemplateReco::Vavilov parameters mpv/sigmaQ/kappa = " << mpv << "/" << sigmaQ << "/" << kappa << std::endl;
                }
#else
                assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
#endif
                xvav = ((double)qtotal-mpv)/sigmaQ;
                beta2 = 1.;
                if(use_VVIObj) {                        
//  VVIObj is a private port of CERNLIB VVIDIS
                   VVIObj vvidist(kappa, beta2, 1);
                   prvav = vvidist.fcn(xvav);                   
                } else {
//  Use faster but less accurate TMath Vavilov distribution function
                        prvav = TMath::VavilovI(xvav, kappa, beta2);
      }
//  Change to upper tail probability
//              if(prvav > 0.5) prvav = 1. - prvav;
//              probQ = (float)(2.*prvav);
                probQ = 1. - prvav;
                if(probQ < probQmin) {probQ = probQmin;}
        } else {
                probQ = -1;
        }
        
        return 0;
} // StripTempReco2D