BeamSpotProducer/src/FcnBeamSpotFitPV.cc

0001 #include "RecoVertex/BeamSpotProducer/interface/FcnBeamSpotFitPV.h"
0002 #include "Math/SVector.h"
0003 #include "Math/SMatrix.h"
0004 // #include "Math/SMatrixDfwd.h"
0005 //#include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 using namespace std;
0008
0009 //
0010 // constructor from vertex data
0011 //
0012 FcnBeamSpotFitPV::FcnBeamSpotFitPV(const vector<BeamSpotFitPVData>& data) : data_(data), errorDef_(1.) {
0013   //
0014   // default: no additional cuts
0015   //
0016   lowerLimits_[0] = lowerLimits_[1] = lowerLimits_[2] = -1.e30;
0017   upperLimits_[0] = upperLimits_[1] = upperLimits_[2] = 1.e30;
0018 }
0019
0020 void FcnBeamSpotFitPV::setLimits(float xmin, float xmax, float ymin, float ymax, float zmin, float zmax) {
0021   lowerLimits_[0] = xmin;
0022   lowerLimits_[1] = ymin;
0023   lowerLimits_[2] = zmin;
0024   upperLimits_[0] = xmax;
0025   upperLimits_[1] = ymax;
0026   upperLimits_[2] = zmax;
0027 }
0028
0029 unsigned int FcnBeamSpotFitPV::nrOfVerticesUsed() const {
0030   //
0031   // count vertices imposing the current limits
0032   //
0033   unsigned int nVtx = 0;
0034   double v1(0);
0035   double v2(0);
0036   double v3(0);
0037   for (vector<BeamSpotFitPVData>::const_iterator ipv = data_.begin(); ipv != data_.end(); ++ipv) {
0038     v1 = (*ipv).position[0];
0039     if (v1 < lowerLimits_[0] || v1 > upperLimits_[0])
0040       continue;
0041     v2 = (*ipv).position[1];
0042     if (v2 < lowerLimits_[1] || v2 > upperLimits_[1])
0043       continue;
0044     v3 = (*ipv).position[2];
0045     if (v3 < lowerLimits_[2] || v3 > upperLimits_[2])
0046       continue;
0047
0048     ++nVtx;
0049   }
0050
0051   return nVtx;
0052 }
0053
0054 double FcnBeamSpotFitPV::operator()(const std::vector<double>& pars) const {
0055   //
0056   // fit parameters
0057   //
0058   double vb1 = pars[0];
0059   double vb2 = pars[1];
0060   double vb3 = pars[2];
0061   double sigb1 = pars[3];
0062   double corrb12 = pars[4];
0063   double sigb2 = pars[5];
0064   double dxdz = pars[6];
0065   double dydz = pars[7];
0066   double sigb3 = pars[8];
0067   double escale = pars[9];
0068   //
0069   // covariance matrix of the beamspot distribution
0070   //
0071   typedef ROOT::Math::SVector<double, 3> Vector3D;
0072   typedef ROOT::Math::SMatrix<double, 3, 3, ROOT::Math::MatRepSym<double, 3> > Matrix3D;
0073   Matrix3D covb;
0074   double varb1 = sigb1 * sigb1;
0075   double varb2 = sigb2 * sigb2;
0076   double varb3 = sigb3 * sigb3;
0077   // parametrisation: rotation (dx/dz, dy/dz); covxy
0078   covb(0, 0) = varb1;
0079   covb(1, 0) = covb(0, 1) = corrb12 * sigb1 * sigb2;
0080   covb(1, 1) = varb2;
0081   covb(2, 0) = covb(0, 2) = dxdz * (varb3 - varb1) - dydz * covb(1, 0);
0082   covb(2, 1) = covb(1, 2) = dydz * (varb3 - varb2) - dxdz * covb(1, 0);
0083   covb(2, 2) = varb3;
0084
0085   //
0086   // calculation of the likelihood function
0087   //
0088   double sumLL(0);
0089   double v1(0);
0090   double v2(0);
0091   double v3(0);
0092   double ev1(0);
0093   double corr12(0);
0094   double ev2(0);
0095   double corr13(0);
0096   double corr23(0);
0097   double ev3(0);
0098   //
0099   // vertex - beamspot difference and total covariance matrix
0100   //
0101   Vector3D dv;
0102   Matrix3D cov;
0103   Matrix3D wgt;
0104   //
0105   // iteration over vertices
0106   //
0107   for (vector<BeamSpotFitPVData>::const_iterator ipv = data_.begin(); ipv != data_.end(); ++ipv) {
0108     //
0109     // additional selection
0110     //
0111     v1 = (*ipv).position[0];
0112     if (v1 < lowerLimits_[0] || v1 > upperLimits_[0])
0113       continue;
0114     v2 = (*ipv).position[1];
0115     if (v2 < lowerLimits_[1] || v2 > upperLimits_[1])
0116       continue;
0117     v3 = (*ipv).position[2];
0118     if (v3 < lowerLimits_[2] || v3 > upperLimits_[2])
0119       continue;
0120     //
0121     // vertex errors (after scaling) and correlations
0122     //
0123     ev1 = (*ipv).posError[0];
0124     corr12 = (*ipv).posCorr[0];
0125     ev2 = (*ipv).posError[1];
0126     corr13 = (*ipv).posCorr[1];
0127     corr23 = (*ipv).posCorr[2];
0128     ev3 = (*ipv).posError[2];
0129     ev1 *= escale;
0130     ev2 *= escale;
0131     ev3 *= escale;
0132     //
0133     // vertex covariance matrix
0134     //
0135     cov(0, 0) = ev1 * ev1;
0136     cov(1, 0) = cov(0, 1) = ev1 * ev2 * corr12;
0137     cov(1, 1) = ev2 * ev2;
0138     cov(2, 0) = cov(0, 2) = ev1 * ev3 * corr13;
0139     cov(2, 1) = cov(1, 2) = ev2 * ev3 * corr23;
0140     cov(2, 2) = ev3 * ev3;
0141     //
0142     // total covariance and weight matrix
0143     //
0144     cov += covb;
0145     int ifail;
0146     wgt = cov.Inverse(ifail);
0147     if (ifail) {
0148       //edm::LogWarning("FcnBeamSpotFitPV")
0149       cout << "Inversion failed" << endl;
0150       return -1.e30;
0151     }
0152     //
0153     // difference of vertex and beamspot positions
0154     //
0155     dv(0) = v1 - vb1;
0156     dv(1) = v2 - vb2;
0157     dv(2) = v3 - vb3;
0158     //
0159     // exponential term and normalization
0160     // (neglecting the additional cut)
0161     //
0162     sumLL += ROOT::Math::Similarity(dv, wgt);
0163     double det;
0164     cov.Det2(det);
0165     sumLL += log(det);
0166   }
0167
0168   return sumLL;
0169 }