#include <memory>
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "DQMServices/Core/interface/DQMStore.h"
#include "DQMServices/Core/interface/MonitorElement.h"
#include <TF3.h>
#include <iostream>
#include <fstream>
#include <vector>

Classes
struct	VertexType
class	Vx3DHLTAnalyzer
Defines
#define	DIM 3
Functions
void	Gauss3DFunc (int &, double , double &fval, double par, int)
Variables
bool	considerVxCovariance
unsigned int	counterVx
double	maxLongLength
double	maxTransRadius
double	pi
std::vector< VertexType >	Vertices
double	VxErrCorr
double	xPos
double	yPos
double	zPos

Define Documentation

#define DIM 3

Definition at line 48 of file Vx3DHLTAnalyzer.h.

Referenced by Vx3DHLTAnalyzer::analyze(), and Gauss3DFunc().

Function Documentation

void Gauss3DFunc	(	int &	,
		double *	,
		double &	fval,
		double *	par,
		int
	)

Definition at line 176 of file Vx3DHLTAnalyzer.cc.

References considerVxCovariance, counterVx, DIM, i, funct::log(), maxLongLength, maxTransRadius, pi, mathSSE::sqrt(), Vertices, VxErrCorr, ExpressReco_HICollisions_FallBack::x, xPos, ExpressReco_HICollisions_FallBack::y, yPos, 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++;
        }
    }
  
  fval = sumlog;
}