Go to the documentation of this file.00001
00014 #include "RecoVertex/BeamSpotProducer/interface/BSpdfsFcn.h"
00015 #include "TMath.h"
00016
00017 #include <cmath>
00018 #include <vector>
00019
00020
00021 double BSpdfsFcn::PDFGauss_d(double z, double d, double sigmad,
00022 double phi, const std::vector<double>& parms) const {
00023
00024
00025
00026
00027
00028 double fsqrt2pi = sqrt(2.* TMath::Pi());
00029
00030 double sig = sqrt(parms[fPar_SigmaBeam]*parms[fPar_SigmaBeam] +
00031 sigmad*sigmad);
00032
00033 double dprime = d - ( ( parms[fPar_X0] + z*parms[fPar_dxdz] )*sin(phi)
00034 - ( parms[fPar_Y0] + z*parms[fPar_dydz] )*cos(phi) );
00035
00036 double result = (exp(-(dprime*dprime)/(2.0*sig*sig)))/(sig*fsqrt2pi);
00037
00038 return result;
00039
00040 }
00041
00042
00043 double BSpdfsFcn::PDFGauss_d_resolution(double z, double d, double phi, double pt, const std::vector<double>& parms) const {
00044
00045
00046
00047
00048
00049
00050 double fsqrt2pi = sqrt(2.* TMath::Pi());
00051
00052 double sigmad = parms[fPar_c0] + parms[fPar_c1]/pt;
00053
00054 double sig = sqrt(parms[fPar_SigmaBeam]*parms[fPar_SigmaBeam] +
00055 sigmad*sigmad);
00056
00057 double dprime = d - ( ( parms[fPar_X0] + z*parms[fPar_dxdz] )*sin(phi)
00058 - ( parms[fPar_Y0] + z*parms[fPar_dydz] )*cos(phi) );
00059
00060 double result = (exp(-(dprime*dprime)/(2.0*sig*sig)))/(sig*fsqrt2pi);
00061
00062 return result;
00063
00064 }
00065
00066
00067 double BSpdfsFcn::PDFGauss_z(double z, double sigmaz, const std::vector<double>& parms) const {
00068
00069
00070
00071
00072
00073 double fsqrt2pi = sqrt(2.* TMath::Pi());
00074
00075 double sig = sqrt(sigmaz*sigmaz+parms[fPar_SigmaZ]*parms[fPar_SigmaZ]);
00076
00077 double result = (exp(-((z-parms[fPar_Z0])*(z-parms[fPar_Z0]))/(2.0*sig*sig)))/(sig*fsqrt2pi);
00078
00079
00080 return result;
00081
00082 }
00083
00084
00085
00086
00087
00088 double BSpdfsFcn::operator() (const std::vector<double>& params) const {
00089
00090 double f = 0.0;
00091
00092
00093
00094 std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin();
00095
00096 double pdf = 0;
00097
00098 for( iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {
00099
00100
00101 if (fusepdfs == "PDFGauss_z") {
00102 pdf = PDFGauss_z( iparam->z0(), iparam->sigz0(),params);
00103 }
00104 else if (fusepdfs == "PDFGauss_d") {
00105 pdf = PDFGauss_d( iparam->z0(), iparam->d0(),
00106 iparam->sigd0(), iparam->phi0(),params);
00107 }
00108 else if (fusepdfs == "PDFGauss_d_resolution") {
00109 pdf = PDFGauss_d_resolution( iparam->z0(), iparam->d0(),
00110 iparam->phi0(), iparam->pt(),params);
00111 }
00112 else if (fusepdfs == "PDFGauss_d*PDFGauss_z") {
00113
00114 pdf = PDFGauss_d( iparam->z0(), iparam->d0(),
00115 iparam->sigd0(), iparam->phi0(),params)*
00116 PDFGauss_z( iparam->z0(), iparam->sigz0(),params);
00117 }
00118 else if (fusepdfs == "PDFGauss_d_resolution*PDFGauss_z") {
00119 pdf = PDFGauss_d_resolution( iparam->z0(), iparam->d0(),
00120 iparam->phi0(), iparam->pt(),params)*
00121 PDFGauss_z( iparam->z0(), iparam->sigz0(),params);
00122 }
00123
00124 f = log(pdf) + f;
00125 }
00126
00127 f= -2.0*f;
00128 return f;
00129
00130 }