BSpdfsFcn.cc

Go to the documentation of this file.

#include "RecoVertex/BeamSpotProducer/interface/BSpdfsFcn.h"
#include "TMath.h"

#include <cmath>
#include <vector>

//______________________________________________________________________
double BSpdfsFcn::PDFGauss_d(double z, double d, double sigmad, double phi, const std::vector<double>& parms) const {
  //---------------------------------------------------------------------------
  //  PDF for d0 distribution. This PDF is a simple gaussian in the
  //  beam reference frame.
  //---------------------------------------------------------------------------
  double fsqrt2pi = sqrt(2. * TMath::Pi());

  double sig = sqrt(parms[fPar_SigmaBeam] * parms[fPar_SigmaBeam] + sigmad * sigmad);

  double dprime =
      d - ((parms[fPar_X0] + z * parms[fPar_dxdz]) * sin(phi) - (parms[fPar_Y0] + z * parms[fPar_dydz]) * cos(phi));

  double result = (exp(-(dprime * dprime) / (2.0 * sig * sig))) / (sig * fsqrt2pi);

  return result;
}

//______________________________________________________________________
double BSpdfsFcn::PDFGauss_d_resolution(
    double z, double d, double phi, double pt, const std::vector<double>& parms) const {
  //---------------------------------------------------------------------------
  //  PDF for d0 distribution. This PDF is a simple gaussian in the
  //  beam reference frame. The IP resolution is parametrize by a linear
  //  function as a function of 1/pt.
  //---------------------------------------------------------------------------
  double fsqrt2pi = sqrt(2. * TMath::Pi());

  double sigmad = parms[fPar_c0] + parms[fPar_c1] / pt;

  double sig = sqrt(parms[fPar_SigmaBeam] * parms[fPar_SigmaBeam] + sigmad * sigmad);

  double dprime =
      d - ((parms[fPar_X0] + z * parms[fPar_dxdz]) * sin(phi) - (parms[fPar_Y0] + z * parms[fPar_dydz]) * cos(phi));

  double result = (exp(-(dprime * dprime) / (2.0 * sig * sig))) / (sig * fsqrt2pi);

  return result;
}

//______________________________________________________________________
double BSpdfsFcn::PDFGauss_z(double z, double sigmaz, const std::vector<double>& parms) const {
  //---------------------------------------------------------------------------
  //  PDF for z-vertex distribution. This distribution
  // is parametrized by a simple normalized gaussian distribution.
  //---------------------------------------------------------------------------
  double fsqrt2pi = sqrt(2. * TMath::Pi());

  double sig = sqrt(sigmaz * sigmaz + parms[fPar_SigmaZ] * parms[fPar_SigmaZ]);
  //double sig = sqrt(sigmaz*sigmaz+parms[1]*parms[1]);
  double result = (exp(-((z - parms[fPar_Z0]) * (z - parms[fPar_Z0])) / (2.0 * sig * sig))) / (sig * fsqrt2pi);
  //double result = (exp(-((z-parms[0])*(z-parms[0]))/(2.0*sig*sig)))/(sig*fsqrt2pi);

  return result;
}

//______________________________________________________________________
double BSpdfsFcn::operator()(const std::vector<double>& params) const {
  double f = 0.0;

  //std::cout << "fusepdfs=" << fusepdfs << " params.size="<<params.size() << std::endl;

  std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin();

  double pdf = 0;

  for (iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {
    if (fusepdfs == "PDFGauss_z") {
      pdf = PDFGauss_z(iparam->z0(), iparam->sigz0(), params);
    } else if (fusepdfs == "PDFGauss_d") {
      pdf = PDFGauss_d(iparam->z0(), iparam->d0(), iparam->sigd0(), iparam->phi0(), params);
    } else if (fusepdfs == "PDFGauss_d_resolution") {
      pdf = PDFGauss_d_resolution(iparam->z0(), iparam->d0(), iparam->phi0(), iparam->pt(), params);
    } else if (fusepdfs == "PDFGauss_d*PDFGauss_z") {
      //std::cout << "pdf= " << pdf << std::endl;
      pdf = PDFGauss_d(iparam->z0(), iparam->d0(), iparam->sigd0(), iparam->phi0(), params) *
            PDFGauss_z(iparam->z0(), iparam->sigz0(), params);
    } else if (fusepdfs == "PDFGauss_d_resolution*PDFGauss_z") {
      pdf = PDFGauss_d_resolution(iparam->z0(), iparam->d0(), iparam->phi0(), iparam->pt(), params) *
            PDFGauss_z(iparam->z0(), iparam->sigz0(), params);
    }

    f = log(pdf) + f;
  }

  f = -2.0 * f;
  return f;
}