BeamSpotProducer/src/BSpdfsFcn.cc

0001
0002 /**_________________________________________________________________
0003    class:   BSpdfsFcn.cc
0004    package: RecoVertex/BeamSpotProducer
0005
0006
0007
0008  author: Francisco Yumiceva, Fermilab (yumiceva@fnal.gov)
0009
0010
0011 ________________________________________________________________**/
0012
0013 #include "RecoVertex/BeamSpotProducer/interface/BSpdfsFcn.h"
0014 #include "TMath.h"
0015
0016 #include <cmath>
0017 #include <vector>
0018
0019 //______________________________________________________________________
0020 double BSpdfsFcn::PDFGauss_d(double z, double d, double sigmad, double phi, const std::vector<double>& parms) const {
0021   //---------------------------------------------------------------------------
0022   //  PDF for d0 distribution. This PDF is a simple gaussian in the
0023   //  beam reference frame.
0024   //---------------------------------------------------------------------------
0025   double fsqrt2pi = sqrt(2. * TMath::Pi());
0026
0027   double sig = sqrt(parms[fPar_SigmaBeam] * parms[fPar_SigmaBeam] + sigmad * sigmad);
0028
0029   double dprime =
0030       d - ((parms[fPar_X0] + z * parms[fPar_dxdz]) * sin(phi) - (parms[fPar_Y0] + z * parms[fPar_dydz]) * cos(phi));
0031
0032   double result = (exp(-(dprime * dprime) / (2.0 * sig * sig))) / (sig * fsqrt2pi);
0033
0034   return result;
0035 }
0036
0037 //______________________________________________________________________
0038 double BSpdfsFcn::PDFGauss_d_resolution(
0039     double z, double d, double phi, double pt, const std::vector<double>& parms) const {
0040   //---------------------------------------------------------------------------
0041   //  PDF for d0 distribution. This PDF is a simple gaussian in the
0042   //  beam reference frame. The IP resolution is parametrize by a linear
0043   //  function as a function of 1/pt.
0044   //---------------------------------------------------------------------------
0045   double fsqrt2pi = sqrt(2. * TMath::Pi());
0046
0047   double sigmad = parms[fPar_c0] + parms[fPar_c1] / pt;
0048
0049   double sig = sqrt(parms[fPar_SigmaBeam] * parms[fPar_SigmaBeam] + sigmad * sigmad);
0050
0051   double dprime =
0052       d - ((parms[fPar_X0] + z * parms[fPar_dxdz]) * sin(phi) - (parms[fPar_Y0] + z * parms[fPar_dydz]) * cos(phi));
0053
0054   double result = (exp(-(dprime * dprime) / (2.0 * sig * sig))) / (sig * fsqrt2pi);
0055
0056   return result;
0057 }
0058
0059 //______________________________________________________________________
0060 double BSpdfsFcn::PDFGauss_z(double z, double sigmaz, const std::vector<double>& parms) const {
0061   //---------------------------------------------------------------------------
0062   //  PDF for z-vertex distribution. This distribution
0063   // is parametrized by a simple normalized gaussian distribution.
0064   //---------------------------------------------------------------------------
0065   double fsqrt2pi = sqrt(2. * TMath::Pi());
0066
0067   double sig = sqrt(sigmaz * sigmaz + parms[fPar_SigmaZ] * parms[fPar_SigmaZ]);
0068   //double sig = sqrt(sigmaz*sigmaz+parms[1]*parms[1]);
0069   double result = (exp(-((z - parms[fPar_Z0]) * (z - parms[fPar_Z0])) / (2.0 * sig * sig))) / (sig * fsqrt2pi);
0070   //double result = (exp(-((z-parms[0])*(z-parms[0]))/(2.0*sig*sig)))/(sig*fsqrt2pi);
0071
0072   return result;
0073 }
0074
0075 //______________________________________________________________________
0076 double BSpdfsFcn::operator()(const std::vector<double>& params) const {
0077   double f = 0.0;
0078
0079   //std::cout << "fusepdfs=" << fusepdfs << " params.size="<<params.size() << std::endl;
0080
0081   std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin();
0082
0083   double pdf = 0;
0084
0085   for (iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {
0086     if (fusepdfs == "PDFGauss_z") {
0087       pdf = PDFGauss_z(iparam->z0(), iparam->sigz0(), params);
0088     } else if (fusepdfs == "PDFGauss_d") {
0089       pdf = PDFGauss_d(iparam->z0(), iparam->d0(), iparam->sigd0(), iparam->phi0(), params);
0090     } else if (fusepdfs == "PDFGauss_d_resolution") {
0091       pdf = PDFGauss_d_resolution(iparam->z0(), iparam->d0(), iparam->phi0(), iparam->pt(), params);
0092     } else if (fusepdfs == "PDFGauss_d*PDFGauss_z") {
0093       //std::cout << "pdf= " << pdf << std::endl;
0094       pdf = PDFGauss_d(iparam->z0(), iparam->d0(), iparam->sigd0(), iparam->phi0(), params) *
0095             PDFGauss_z(iparam->z0(), iparam->sigz0(), params);
0096     } else if (fusepdfs == "PDFGauss_d_resolution*PDFGauss_z") {
0097       pdf = PDFGauss_d_resolution(iparam->z0(), iparam->d0(), iparam->phi0(), iparam->pt(), params) *
0098             PDFGauss_z(iparam->z0(), iparam->sigz0(), params);
0099     }
0100
0101     f = log(pdf) + f;
0102   }
0103
0104   f = -2.0 * f;
0105   return f;
0106 }