-
Notifications
You must be signed in to change notification settings - Fork 4.4k
/
Copy pathBSpdfsFcn.cc
106 lines (77 loc) · 4.29 KB
/
BSpdfsFcn.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
/**_________________________________________________________________
class: BSpdfsFcn.cc
package: RecoVertex/BeamSpotProducer
author: Francisco Yumiceva, Fermilab ([email protected])
________________________________________________________________**/
#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;
}