Belle II Software light-2406-ragdoll
PxPyPzMFitObject.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include "analysis/OrcaKinFit/PxPyPzMFitObject.h"
10#include <framework/logging/Logger.h>
11#include <cmath>
12#include <cassert>
13#include <iostream>
14
15using std::sqrt;
16using std::sin;
17using std::cos;
18
19namespace Belle2 {
24 namespace OrcaKinFit {
25
26// constructor
27 PxPyPzMFitObject::PxPyPzMFitObject(CLHEP::HepLorentzVector& particle, const CLHEP::HepSymMatrix& covmatrix)
28 : cachevalid(false), chi2(0), dEdpx(0), dEdpy(0), dEdpz(0),
29 dE2dpxdpx(0), dE2dpxdpy(0), dE2dpxdpz(0), dE2dpydpy(0), dE2dpydpz(0), dE2dpzdpz(0)
30 {
31
32 assert(int(NPAR) <= int(BaseDefs::MAXPAR));
33
34 initCov();
35 setMass(particle.m());
36
37 setParam(0, particle.px(), true);
38 setParam(1, particle.py(), true);
39 setParam(2, particle.pz(), true);
40
41 setMParam(0, particle.px());
42 setMParam(1, particle.py());
43 setMParam(2, particle.pz());
44
45 //set covariance matrix
46 for (int i = 0; i < int(NPAR); i++) {
47 for (int j = 0; j < int(NPAR); j++) {
48 setCov(i, j, covmatrix[i][j]);
49 }
50 }
51
52 invalidateCache();
53 }
54
55// destructor
56 PxPyPzMFitObject::~PxPyPzMFitObject() = default;
57
58 PxPyPzMFitObject::PxPyPzMFitObject(const PxPyPzMFitObject& rhs)
59 : ParticleFitObject(rhs), cachevalid(false), chi2(0), dEdpx(0), dEdpy(0), dEdpz(0),
60 dE2dpxdpx(0), dE2dpxdpy(0), dE2dpxdpz(0), dE2dpydpy(0), dE2dpydpz(0), dE2dpzdpz(0)
61 {
62
64 }
65
67 {
68 if (this != &rhs) {
69 assign(rhs); // calls virtual function assign of derived class
70 }
71 return *this;
72 }
73
75 {
76 return new PxPyPzMFitObject(*this);
77 }
78
80 {
81 if (const auto* psource = dynamic_cast<const PxPyPzMFitObject*>(&source)) {
82 if (psource != this) {
84 // only mutable data members, need not to be copied, if cache is invalid
85 }
86 } else {
87 assert(0);
88 }
89 return *this;
90 }
91
92 const char* PxPyPzMFitObject::getParamName(int ilocal) const
93 {
94 switch (ilocal) {
95 case 0: return "px";
96 case 1: return "py";
97 case 2: return "pz";
98 }
99 return "undefined";
100 }
101
102
103 bool PxPyPzMFitObject::updateParams(double p[], int idim)
104 {
106
107 int ipx = getGlobalParNum(0);
108 int ipy = getGlobalParNum(1);
109 int ipz = getGlobalParNum(2);
110 assert(ipx >= 0 && ipx < idim);
111 assert(ipy >= 0 && ipy < idim);
112 assert(ipz >= 0 && ipz < idim);
113
114 double px = p[ipx];
115 double py = p[ipy];
116 double pz = p[ipz];
117
118 bool result = ((px - par[0]) * (px - par[0]) > eps2 * cov[0][0]) ||
119 ((py - par[1]) * (py - par[1]) > eps2 * cov[1][1]) ||
120 ((pz - par[2]) * (pz - par[2]) > eps2 * cov[2][2]);
121
122 par[0] = px;
123 par[1] = py;
124 par[2] = pz;
125 p[ipx] = par[0];
126 p[ipy] = par[1];
127 p[ipz] = par[2];
128 return result;
129 }
130
131// these depend on actual parametrisation!
132
133 double PxPyPzMFitObject::getDPx(int ilocal) const
134 {
135 assert(ilocal >= 0 && ilocal < NPAR);
136 if (!cachevalid) updateCache();
137 switch (ilocal) {
138 case 0: return 1.0;
139 case 1: return 0.0;
140 case 2: return 0.0;
141 }
142 return 0;
143 }
144
145 double PxPyPzMFitObject::getDPy(int ilocal) const
146 {
147 assert(ilocal >= 0 && ilocal < NPAR);
148 if (!cachevalid) updateCache();
149 switch (ilocal) {
150 case 0: return 0.0;
151 // case 1: return dpydtheta;
152 case 1: return 1.0;
153 case 2: return 0.0;
154 }
155 return 0;
156 }
157
158 double PxPyPzMFitObject::getDPz(int ilocal) const
159 {
160 assert(ilocal >= 0 && ilocal < NPAR);
161 if (!cachevalid) updateCache();
162 switch (ilocal) {
163 case 0: return 0.0;
164 case 1: return 0.0;
165 case 2: return 1.0;
166 }
167 return 0;
168 }
169
170 double PxPyPzMFitObject::getDE(int ilocal) const
171 {
172 assert(ilocal >= 0 && ilocal < NPAR);
173 if (!cachevalid) updateCache();
174 switch (ilocal) {
175 case 0: return dEdpx;
176 case 1: return dEdpy;
177 case 2: return dEdpz;
178 }
179 return 0;
180 }
181
182 double PxPyPzMFitObject::getFirstDerivative_Meta_Local(int iMeta, int ilocal, int metaSet) const
183 {
184 // iMeta = intermediate variable (i.e. E,px,py,pz)
185 // ilocal = local variable (ptinv, theta, phi)
186 // metaSet = which set of intermediate varlables
187
188 assert(metaSet == 0); // only defined for E,px,py,pz
189
190 switch (iMeta) {
191 case 0: // E
192 return getDE(ilocal);
193 break;
194 case 1: // Px
195 return getDPx(ilocal);
196 break;
197 case 2: // Py
198 return getDPy(ilocal);
199 break;
200 case 3: // Pz
201 return getDPz(ilocal);
202 break;
203 default:
204 assert(0);
205 }
206 }
207
208 double PxPyPzMFitObject::getSecondDerivative_Meta_Local(int iMeta, int ilocal, int jlocal, int metaSet) const
209 {
210 assert(metaSet == 0);
211 if (!cachevalid) updateCache();
212 if (jlocal < ilocal) {
213 int temp = jlocal;
214 jlocal = ilocal;
215 ilocal = temp;
216 }
217
218 switch (iMeta) {
219 case 0: // E
220 if (ilocal == 0 && jlocal == 0) return dE2dpxdpx;
221 else if (ilocal == 0 && jlocal == 1) return dE2dpxdpy;
222 else if (ilocal == 0 && jlocal == 2) return dE2dpxdpz;
223 else if (ilocal == 1 && jlocal == 1) return dE2dpydpy;
224 else if (ilocal == 1 && jlocal == 2) return dE2dpydpz;
225 else if (ilocal == 2 && jlocal == 2) return dE2dpzdpz;
226 else return 0;
227 break;
228 case 1: // px
229 return 0;
230 break;
231 case 2: // py
232 return 0;
233 break;
234 case 3: // pz
235 return 0;
236 break;
237 default:
238 assert(0);
239 }
240 }
241
242 void PxPyPzMFitObject::updateCache() const
243 {
244
245 // get the basic quantities
246 double px = par[0];
247 double py = par[1];
248 double pz = par[2];
249 double px2 = px * px;
250 double py2 = py * py;
251 double pz2 = pz * pz;
252 double p = std::sqrt(px2 + py2 + pz2);
253 assert(p != 0);
254
255 double mass2 = mass * mass;
256 double e2 = px2 + py2 + pz2 + mass2;
257 double e = std::sqrt(e2);
258
259 // get the first derivatives
260 dEdpx = px / e;
261 dEdpy = py / e;
262 dEdpz = pz / e;
263
264 // get second derivatives
265 double e32 = std::pow(e, 1.5);
266 dE2dpxdpx = (py2 + pz2 + mass2) / (e32);
267 dE2dpxdpy = -(px * py) / (e32);
268 dE2dpxdpz = -(px * pz) / (e32);
269 dE2dpydpy = (px2 + pz2 + mass2) / (e32);
270 dE2dpydpz = -(py * pz) / (e32);
271 dE2dpzdpz = (px2 + py2 + mass2) / (e32);
272
273 fourMomentum.setValues(e, px, py, pz);
274
275 cachevalid = true;
276 }
277 }// end OrcaKinFit namespace
279} // end Belle2 namespace
280
virtual int getGlobalParNum(int ilocal) const
Get global parameter number of parameter ilocal.
double par[BaseDefs::MAXPAR]
fit parameters
void invalidateCache() const
invalidate any cached quantities
double cov[BaseDefs::MAXPAR][BaseDefs::MAXPAR]
local covariance matrix
virtual ParticleFitObject & assign(const BaseFitObject &source) override
Assign from anther object, if of same type.
virtual double getDPx(int ilocal) const override
Return d p_x / d par_ilocal (derivative of px w.r.t. local parameter ilocal)
virtual double getDPy(int ilocal) const override
Return d p_y / d par_ilocal (derivative of py w.r.t. local parameter ilocal)
virtual PxPyPzMFitObject & assign(const BaseFitObject &source) override
Assign from anther object, if of same type.
virtual double getDE(int ilocal) const override
Return d E / d par_ilocal (derivative of E w.r.t. local parameter ilocal)
virtual PxPyPzMFitObject * copy() const override
Return a new copy of itself.
virtual const char * getParamName(int ilocal) const override
Get name of parameter ilocal.
virtual double getDPz(int ilocal) const override
Return d p_z / d par_ilocal (derivative of pz w.r.t. local parameter ilocal)
virtual bool updateParams(double p[], int idim) override
Read values from global vector, readjust vector; return: significant change.
PxPyPzMFitObject & operator=(const PxPyPzMFitObject &rhs)
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24