FcnBeamSpotFitPV.cc

Go to the documentation of this file.

#include "RecoVertex/BeamSpotProducer/interface/FcnBeamSpotFitPV.h"
#include "Math/SVector.h"
#include "Math/SMatrix.h"
// #include "Math/SMatrixDfwd.h"
//#include "FWCore/MessageLogger/interface/MessageLogger.h"

using namespace std;

//
// constructor from vertex data
//
FcnBeamSpotFitPV::FcnBeamSpotFitPV(const vector<BeamSpotFitPVData>& data) : 
  data_(data), errorDef_(1.) { 
  //
  // default: no additional cuts
  //
  lowerLimits_[0] = lowerLimits_[1] = lowerLimits_[2] = -1.e30;
  upperLimits_[0] = upperLimits_[1] = upperLimits_[2] =  1.e30;
}

void 
FcnBeamSpotFitPV::setLimits (float xmin, float xmax,
              float ymin, float ymax,
              float zmin, float zmax) 
{
  lowerLimits_[0] = xmin;
  lowerLimits_[1] = ymin;
  lowerLimits_[2] = zmin;
  upperLimits_[0] = xmax;
  upperLimits_[1] = ymax;
  upperLimits_[2] = zmax;
}

unsigned int
FcnBeamSpotFitPV::nrOfVerticesUsed () const
{
  //
  // count vertices imposing the current limits
  //
  unsigned int nVtx = 0;
  double v1(0);
  double v2(0);
  double v3(0);
  for ( vector<BeamSpotFitPVData>::const_iterator ipv=data_.begin();
    ipv!=data_.end(); ++ipv ) {
    v1 = (*ipv).position[0];
    if ( v1<lowerLimits_[0] || v1>upperLimits_[0] )  continue;
    v2 = (*ipv).position[1];
    if ( v2<lowerLimits_[1] || v2>upperLimits_[1] )  continue;
    v3 = (*ipv).position[2];
    if ( v3<lowerLimits_[2] || v3>upperLimits_[2] )  continue;

    ++nVtx;
  }

  return nVtx;
}

double
FcnBeamSpotFitPV::operator() (const std::vector<double>& pars) const
{
  //
  // fit parameters
  //
  double vb1 = pars[0];
  double vb2 = pars[1];
  double vb3 = pars[2];
  double sigb1 = pars[3]; 
  double corrb12 = pars[4];
  double sigb2 = pars[5];
  double dxdz = pars[6];
  double dydz = pars[7];
  double sigb3 = pars[8];
  double escale = pars[9];
  //
  // covariance matrix of the beamspot distribution
  //
  typedef ROOT::Math::SVector<double,3> Vector3D;
  typedef ROOT::Math::SMatrix<double,3,3,ROOT::Math::MatRepSym<double,3> > Matrix3D;
  Matrix3D covb;
  double varb1 = sigb1*sigb1;
  double varb2 = sigb2*sigb2;
  double varb3 = sigb3*sigb3;
// parametrisation: rotation (dx/dz, dy/dz); covxy
  covb(0,0) = varb1;
  covb(1,0) = covb(0,1) = corrb12*sigb1*sigb2;
  covb(1,1) = varb2;
  covb(2,0) = covb(0,2) = dxdz*(varb3-varb1)-dydz*covb(1,0);
  covb(2,1) = covb(1,2) = dydz*(varb3-varb2)-dxdz*covb(1,0);
  covb(2,2) = varb3;

  //
  // calculation of the likelihood function
  //
  double sumLL(0);
  double v1(0);
  double v2(0);
  double v3(0);
  double ev1(0);
  double corr12(0);
  double ev2(0);
  double corr13(0);
  double corr23(0);
  double ev3(0);
  //
  // vertex - beamspot difference and total covariance matrix
  //
  Vector3D dv;
  Matrix3D cov;
  Matrix3D wgt;
  //
  // iteration over vertices
  //
  for ( vector<BeamSpotFitPVData>::const_iterator ipv=data_.begin();
    ipv!=data_.end(); ++ipv ) {
    //
    // additional selection
    //
    v1 = (*ipv).position[0];
    if ( v1<lowerLimits_[0] || v1>upperLimits_[0] )  continue;
    v2 = (*ipv).position[1];
    if ( v2<lowerLimits_[1] || v2>upperLimits_[1] )  continue;
    v3 = (*ipv).position[2];
    if ( v3<lowerLimits_[2] || v3>upperLimits_[2] )  continue;
    //
    // vertex errors (after scaling) and correlations
    //
    ev1 = (*ipv).posError[0];
    corr12 = (*ipv).posCorr[0];
    ev2 = (*ipv).posError[1];
    corr13 = (*ipv).posCorr[1];
    corr23 = (*ipv).posCorr[2];
    ev3 = (*ipv).posError[2];
    ev1 *= escale;
    ev2 *= escale;
    ev3 *= escale;
    //
    // vertex covariance matrix
    //
    cov(0,0) = ev1*ev1;
    cov(1,0) = cov(0,1) = ev1*ev2*corr12;
    cov(1,1) = ev2*ev2;
    cov(2,0) = cov(0,2) = ev1*ev3*corr13;
    cov(2,1) = cov(1,2) = ev2*ev3*corr23;
    cov(2,2) = ev3*ev3;
    //
    // total covariance and weight matrix
    //
    cov += covb;
    int ifail;
    wgt = cov.Inverse(ifail);
    if ( ifail ) {
      //edm::LogWarning("FcnBeamSpotFitPV") 
      cout << "Inversion failed" << endl;
      return -1.e30;
    }
    //
    // difference of vertex and beamspot positions
    //
    dv(0) = v1 - vb1;
    dv(1) = v2 - vb2;
    dv(2) = v3 - vb3;
    //
    // exponential term and normalization
    // (neglecting the additional cut)
    //
    sumLL += ROOT::Math::Similarity(dv,wgt);
    double det;
    cov.Det2(det);
    sumLL += log(det);
  }

  return sumLL;
}