File indexing completed on 2024-10-10 23:05:30
0001
0002
0003
0004
0005
0006
0007
0008
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
0023
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
0042
0043
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
0063
0064
0065 double fsqrt2pi = sqrt(2. * TMath::Pi());
0066
0067 double sig = sqrt(sigmaz * sigmaz + parms[fPar_SigmaZ] * parms[fPar_SigmaZ]);
0068
0069 double result = (exp(-((z - parms[fPar_Z0]) * (z - parms[fPar_Z0])) / (2.0 * sig * sig))) / (sig * fsqrt2pi);
0070
0071
0072 return result;
0073 }
0074
0075
0076 double BSpdfsFcn::operator()(const std::vector<double>& params) const {
0077 double f = 0.0;
0078
0079
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
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 }