Vx3DHLTAnalyzer.cc File Reference

#include "DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Framework/interface/LuminosityBlock.h"
#include "FWCore/Utilities/interface/isFinite.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include <Math/Minimizer.h>
#include <Math/Factory.h>
#include <Math/Functor.h>

Go to the source code of this file.

Functions
	DEFINE_FWK_MODULE (Vx3DHLTAnalyzer)
double	Gauss3DFunc (const double *par)

Function Documentation

DEFINE_FWK_MODULE

(

Vx3DHLTAnalyzer

)

double Gauss3DFunc

(

const double *

par

)

Definition at line 176 of file Vx3DHLTAnalyzer.cc.

References considerVxCovariance, counterVx, DIM, i, fff_deleter::log, maxLongLength, maxTransRadius, pi, mathSSE::sqrt(), Vertices, VxErrCorr, x, xPos, detailsBasic3DVector::y, yPos, detailsBasic3DVector::z, and zPos.

Referenced by Vx3DHLTAnalyzer::MyFit().

{
  double K[DIM][DIM]; // Covariance Matrix
  double M[DIM][DIM]; // K^-1
  double det;
  double sumlog = 0.;

//  par[0] = K(0,0) --> Var[X]
//  par[1] = K(1,1) --> Var[Y]
//  par[2] = K(2,2) --> Var[Z]
//  par[3] = K(0,1) = K(1,0) --> Cov[X,Y]
//  par[4] = K(1,2) = K(2,1) --> Cov[Y,Z] --> dy/dz
//  par[5] = K(0,2) = K(2,0) --> Cov[X,Z] --> dx/dz
//  par[6] = mean x
//  par[7] = mean y
//  par[8] = mean z

  counterVx = 0;
  for (unsigned int i = 0; i < Vertices.size(); i++)
    {
      if ((std::sqrt((Vertices[i].x-xPos)*(Vertices[i].x-xPos) + (Vertices[i].y-yPos)*(Vertices[i].y-yPos)) <= maxTransRadius) &&
      (std::fabs(Vertices[i].z-zPos) <= maxLongLength))
    {
      if (considerVxCovariance == true)
        {
          K[0][0] = std::fabs(par[0]) + VxErrCorr*VxErrCorr * std::fabs(Vertices[i].Covariance[0][0]);
          K[1][1] = std::fabs(par[1]) + VxErrCorr*VxErrCorr * std::fabs(Vertices[i].Covariance[1][1]);
          K[2][2] = std::fabs(par[2]) + VxErrCorr*VxErrCorr * std::fabs(Vertices[i].Covariance[2][2]);
          K[0][1] = K[1][0] = par[3] + VxErrCorr*VxErrCorr * Vertices[i].Covariance[0][1];
          K[1][2] = K[2][1] = par[4]*(std::fabs(par[2])-std::fabs(par[1])) - par[5]*par[3] + VxErrCorr*VxErrCorr * Vertices[i].Covariance[1][2];
          K[0][2] = K[2][0] = par[5]*(std::fabs(par[2])-std::fabs(par[0])) - par[4]*par[3] + VxErrCorr*VxErrCorr * Vertices[i].Covariance[0][2];
        }
      else
        {
          K[0][0] = std::fabs(par[0]);
          K[1][1] = std::fabs(par[1]);
          K[2][2] = std::fabs(par[2]);
          K[0][1] = K[1][0] = par[3];
          K[1][2] = K[2][1] = par[4]*(std::fabs(par[2])-std::fabs(par[1])) - par[5]*par[3];
          K[0][2] = K[2][0] = par[5]*(std::fabs(par[2])-std::fabs(par[0])) - par[4]*par[3];
        }

      det = K[0][0]*(K[1][1]*K[2][2] - K[1][2]*K[1][2]) -
        K[0][1]*(K[0][1]*K[2][2] - K[0][2]*K[1][2]) +
        K[0][2]*(K[0][1]*K[1][2] - K[0][2]*K[1][1]);

      M[0][0] = (K[1][1]*K[2][2] - K[1][2]*K[1][2]) / det;
      M[1][1] = (K[0][0]*K[2][2] - K[0][2]*K[0][2]) / det;
      M[2][2] = (K[0][0]*K[1][1] - K[0][1]*K[0][1]) / det;
      M[0][1] = M[1][0] = (K[0][2]*K[1][2] - K[0][1]*K[2][2]) / det;
      M[1][2] = M[2][1] = (K[0][2]*K[0][1] - K[1][2]*K[0][0]) / det;
      M[0][2] = M[2][0] = (K[0][1]*K[1][2] - K[0][2]*K[1][1]) / det;
      
      sumlog += double(DIM)*std::log(2.*pi) + std::log(std::fabs(det)) +
        (M[0][0]*(Vertices[i].x-par[6])*(Vertices[i].x-par[6]) +
         M[1][1]*(Vertices[i].y-par[7])*(Vertices[i].y-par[7]) +
         M[2][2]*(Vertices[i].z-par[8])*(Vertices[i].z-par[8]) +
         2.*M[0][1]*(Vertices[i].x-par[6])*(Vertices[i].y-par[7]) +
         2.*M[1][2]*(Vertices[i].y-par[7])*(Vertices[i].z-par[8]) +
         2.*M[0][2]*(Vertices[i].x-par[6])*(Vertices[i].z-par[8]));
      
      counterVx++;
    }
    }

  return sumlog;
}