Belle II Software light-2406-ragdoll
JetFitObject.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * Forked from https://github.com/iLCSoft/MarlinKinfit *
6 * *
7 * Further information about the fit engine and the user interface *
8 * provided in MarlinKinfit can be found at *
9 * https://www.desy.de/~blist/kinfit/doc/html/ *
10 * and in the LCNotes LC-TOOL-2009-001 and LC-TOOL-2009-004 available *
11 * from http://www-flc.desy.de/lcnotes/ *
12 * *
13 * See git log for contributors and copyright holders. *
14 * This file is licensed under LGPL-3.0, see LICENSE.md. *
15 **************************************************************************/
16
17#include "analysis/OrcaKinFit/JetFitObject.h"
18#include <framework/logging/Logger.h>
19#include <cmath>
20
21#undef NDEBUG
22#include <cassert>
23
24#include <iostream>
25
26using std::sqrt;
27using std::sin;
28using std::cos;
29using std::endl;
30
31namespace Belle2 {
36 namespace OrcaKinFit {
37
38// constructor
39 JetFitObject::JetFitObject(double E, double theta, double phi,
40 double DE, double Dtheta, double Dphi,
41 double m)
42 : ctheta(0), stheta(0), cphi(0), sphi(0),
43 p2(0), p(0), pt(0), px(0), py(0), pz(0), dpdE(0), dptdE(0),
44 dpxdE(0), dpydE(0), dpzdE(0), dpxdtheta(0), dpydtheta(0), chi2(0)
45 {
46
47 assert(int(NPAR) <= int(BaseDefs::MAXPAR));
48
49 initCov();
50// assert( !isinf(E) ); assert( !isnan(E) );
51// assert( !isinf(theta) ); assert( !isnan(theta) );
52// assert( !isinf(phi) ); assert( !isnan(phi) );
53// assert( !isinf(DE) ); assert( !isnan(DE) );
54// assert( !isinf(Dtheta) ); assert( !isnan(Dtheta) );
55// assert( !isinf(Dphi) ); assert( !isnan(Dphi) );
56// assert( !isinf(m) ); assert( !isnan(m) );
57 setMass(m);
58 adjustEThetaPhi(m, E, theta, phi);
59 setParam(0, E, true);
60 setParam(1, theta, true);
61 setParam(2, phi, true);
62 setMParam(0, E);
63 setMParam(1, theta);
64 setMParam(2, phi);
65 setError(0, DE);
66 setError(1, Dtheta);
67 setError(2, Dphi);
68
69 // parameter 2 repeats every 2*pi
70 paramCycl[2] = 2.*M_PI;
71
72 invalidateCache();
73 B2DEBUG(15, "JetFitObject::JetFitObject: E = " << E);
74 B2DEBUG(15, "JetFitObject::JetFitObject: getParam(0) = " << getParam(0));
75 B2DEBUG(15, "JetFitObject::JetFitObject: " << *this);
76 B2DEBUG(15, "mpar= " << mpar[0] << ", " << mpar[1] << ", " << mpar[2]);
77 }
78
79// destructor
80 JetFitObject::~JetFitObject() = default;
81
82 JetFitObject::JetFitObject(const JetFitObject& rhs)
83 : ParticleFitObject(rhs), ctheta(0), stheta(0), cphi(0), sphi(0),
84 p2(0), p(0), pt(0), px(0), py(0), pz(0), dpdE(0), dptdE(0),
85 dpxdE(0), dpydE(0), dpzdE(0), dpxdtheta(0), dpydtheta(0), chi2(0)
86 {
88 }
89
91 {
92 if (this != &rhs) {
93 assign(rhs); // calls virtual function assign of derived class
94 }
95 return *this;
96 }
97
99 {
100 return new JetFitObject(*this);
101 }
102
104 {
105 if (const auto* psource = dynamic_cast<const JetFitObject*>(&source)) {
106 if (psource != this) {
108 // only mutable data members, need not to be copied, if cache is invalid
109 }
110 } else {
111 assert(0);
112 }
113 return *this;
114 }
115
116 const char* JetFitObject::getParamName(int ilocal) const
117 {
118 switch (ilocal) {
119 case 0: return "E";
120 case 1: return "theta";
121 case 2: return "phi";
122 }
123 return "undefined";
124 }
125
126
127 bool JetFitObject::updateParams(double pp[], int idim)
128 {
130
131 bool result = false;
132 bool mirrored_e = false;
133
134 if (!isParamFixed(0)) {
135 int iE = getGlobalParNum(0);
136 assert(iE >= 0 && iE < idim);
137 double e = pp[iE];
138
139 if (e < 0) {
140 B2INFO("JetFitObject::updateParams: mirrored E!\n");
141 e = -e;
142 mirrored_e = true;
143 }
144
145 double massPlusEpsilon = mass * (1.0000001);
146 if (e < massPlusEpsilon) e = massPlusEpsilon;
147 result = result || ((e - par[0]) * (e - par[0]) > eps2 * cov[0][0]);
148 par[0] = e;
149 pp[iE] = par[0];
150 }
151
152 if (!isParamFixed(1)) {
153 int ith = getGlobalParNum(1);
154 assert(ith >= 0 && ith < idim);
155 double th = pp[ith];
156
157 if (mirrored_e) {
158 th = M_PI - th;
159 }
160
161 result = result || ((th - par[1]) * (th - par[1]) > eps2 * cov[1][1]);
162 par[1] = th;
163 pp[ith] = par[1];
164 }
165
166 if (!isParamFixed(2)) {
167 int iph = getGlobalParNum(2);
168 assert(iph >= 0 && iph < idim);
169 double ph = pp[iph];
170
171 if (mirrored_e) {
172 ph = M_PI + ph;
173 }
174
175 result = result || ((ph - par[2]) * (ph - par[2]) > eps2 * cov[2][2]);
176 par[2] = ph;
177 pp[iph] = par[2];
178 }
179
180
181 return result;
182 }
183
184 double JetFitObject::getDPx(int ilocal) const
185 {
186 assert(ilocal >= 0 && ilocal < NPAR);
187 if (!cachevalid) updateCache();
188 switch (ilocal) {
189 case 0: return dpxdE;
190 case 1: return dpxdtheta;
191 case 2: return -py;
192 }
193 return 0;
194 }
195
196 double JetFitObject::getDPy(int ilocal) const
197 {
198 assert(ilocal >= 0 && ilocal < NPAR);
199 if (!cachevalid) updateCache();
200 switch (ilocal) {
201 case 0: return dpydE;
202 case 1: return dpydtheta;
203 case 2: return px;
204 }
205 return 0;
206 }
207
208 double JetFitObject::getDPz(int ilocal) const
209 {
210 assert(ilocal >= 0 && ilocal < NPAR);
211 if (!cachevalid) updateCache();
212 switch (ilocal) {
213 case 0: return dpzdE;
214 case 1: return -pt;
215 case 2: return 0;
216 }
217 return 0;
218 }
219
220 double JetFitObject::getDE(int ilocal) const
221 {
222 assert(ilocal >= 0 && ilocal < NPAR);
223 switch (ilocal) {
224 case 0: return 1;
225 case 1: return 0;
226 case 2: return 0;
227 }
228 return 0;
229 }
230
231 double JetFitObject::getError(int ilocal) const
232 {
233 assert(ilocal >= 0 && ilocal < NPAR);
234 B2INFO("JetFitObject::getError (ilocal = " << ilocal << ") = " << std::sqrt(cov[ilocal][ilocal]));
235 return std::sqrt(cov[ilocal][ilocal]);
236 }
237
238 double JetFitObject::getCov(int ilocal, int jlocal) const
239 {
240 assert(ilocal >= 0 && ilocal < NPAR);
241 assert(jlocal >= 0 && jlocal < NPAR);
242 return cov[ilocal][jlocal];
243 }
244
245//void JetFitObject::addToDerivatives (double der[], int idim,
246// double efact, double pxfact,
247// double pyfact, double pzfact) const {
248// int i_E = globalParNum[0];
249// int i_theta = globalParNum[1];
250// int i_phi = globalParNum[2];
251// assert (i_E >= 0 && i_E < idim);
252// assert (i_theta >= 0 && i_theta < idim);
253// assert (i_phi >= 0 && i_phi < idim);
254//
255// if (!cachevalid) updateCache();
256// // for numerical accuracy, add up derivatives first,
257// // then add them to global vector
258// double der_E = efact;
259// double der_theta = 0;
260// double der_phi = 0;
261//
262// if (pxfact != 0) {
263// der_E += pxfact*dpxdE;
264// der_theta += pxfact*dpxdtheta;
265// der_phi -= pxfact*py;
266
267
268
269 double JetFitObject::getFirstDerivative_Meta_Local(int iMeta, int ilocal, int metaSet) const
270 {
271
272 assert(metaSet == 0);
273 switch (iMeta) {
274 case 0:
275 return getDE(ilocal);
276 break;
277 case 1:
278 return getDPx(ilocal);
279 break;
280 case 2:
281 return getDPy(ilocal);
282 break;
283 case 3:
284 return getDPz(ilocal);
285 break;
286 default:
287 assert(0);
288
289 }
290 return -999;
291 }
292
293
294 double JetFitObject::getSecondDerivative_Meta_Local(int iMeta, int ilocal, int jlocal, int metaSet) const
295 {
296 assert(metaSet == 0);
297 if (!cachevalid) updateCache();
298
299 if (jlocal < ilocal) {
300 int temp = jlocal;
301 jlocal = ilocal;
302 ilocal = temp;
303 }
304
305 double d2pdE2 = (mass != 0) ? -mass * mass / (p * p * p) : 0;
306 double d2ptdE2 = d2pdE2 * stheta;
307
308 // daniel hasn't checked these, copied from orig code
309 switch (iMeta) {
310
311 case 0:
312 return 0;
313 break;
314 case 1:
315 if (ilocal == 0 && jlocal == 0) return d2ptdE2 * cphi;
316 else if (ilocal == 0 && jlocal == 1) return dpzdE * cphi;
317 else if (ilocal == 0 && jlocal == 2) return -dpydE;
318 else if (ilocal == 1 && jlocal == 1) return -px;
319 else if (ilocal == 1 && jlocal == 2) return -dpydtheta;
320 else if (ilocal == 2 && jlocal == 2) return -px;
321 break;
322 case 2:
323 if (ilocal == 0 && jlocal == 0) return d2ptdE2 * sphi;
324 else if (ilocal == 0 && jlocal == 1) return dpzdE * sphi;
325 else if (ilocal == 0 && jlocal == 2) return dpxdE;
326 else if (ilocal == 1 && jlocal == 1) return -py;
327 else if (ilocal == 1 && jlocal == 2) return dpxdtheta;
328 else if (ilocal == 2 && jlocal == 2) return -py;
329 break;
330 case 3:
331 if (ilocal == 0 && jlocal == 0) return d2pdE2 * ctheta;
332 // else if ( ilocal==0 && jlocal==1 ) return dptdE;
333 else if (ilocal == 0 && jlocal == 1) return -dptdE; // this is "-" in the orig JetFitObject, DJ fixed 2015may27
334 else if (ilocal == 0 && jlocal == 2) return 0;
335 else if (ilocal == 1 && jlocal == 1) return -pz;
336 else if (ilocal == 1 && jlocal == 2) return 0;
337 else if (ilocal == 2 && jlocal == 2) return 0;
338 break;
339 default:
340 assert(0);
341 }
342 return -999;
343
344 }
345
346
347 void JetFitObject::updateCache() const
348 {
349
350 double e = par[0];
351 double theta = par[1];
352 double phi = par[2];
353
354 ctheta = cos(theta);
355 stheta = sin(theta);
356 cphi = cos(phi);
357 sphi = sin(phi);
358
359 p2 = std::abs(e * e - mass * mass);
360 p = std::sqrt(p2);
361 assert(p != 0);
362
363 pt = p * stheta;
364
365 px = pt * cphi;
366 py = pt * sphi;
367 pz = p * ctheta;
368
369 fourMomentum.setValues(e, px, py, pz);
370
371 dpdE = e / p;
372 dptdE = dpdE * stheta;
373 dpxdE = dptdE * cphi;
374 dpydE = dptdE * sphi;
375 dpzdE = dpdE * ctheta;
376
377 dpxdtheta = pz * cphi;
378 dpydtheta = pz * sphi;
379
380 cachevalid = true;
381
382 }
383
384//double JetFitObject::getChi2 () const {
385// if (!cachevalid) updateCache();
386// return chi2;
387//}
388//
389
390 bool JetFitObject::adjustEThetaPhi(const double& m, double& E, double& theta, double& phi)
391 {
392 bool result = false;
393
394 if (E < 0) {
395 B2INFO("JetFitObject::adjustEThetaPhi: mirrored E!\n");
396 E = -E;
397 theta = M_PI - theta;
398 phi = M_PI + phi;
399 result = true;
400 }
401 if (E < m) {
402 E = m;
403 result = true;
404 }
405 if (theta < -M_PI || theta > M_PI) {
406 while (theta < -M_PI) theta += 2 * M_PI;
407 while (theta > M_PI) theta -= 2 * M_PI;
408 result = true;
409 }
410
411 if (theta < 0) {
412 B2INFO("JetFitObject::adjustEThetaPhi: mirrored theta!\n");
413 theta = -theta;
414 phi = phi > 0 ? phi - M_PI : phi + M_PI;
415 result = true;
416 } else if (theta > M_PI) {
417 B2INFO("JetFitObject::adjustEThetaPhi: mirrored theta!\n");
418 theta = 2 * M_PI - theta;
419 phi = phi > 0 ? phi - M_PI : phi + M_PI;
420 result = true;
421 }
422 if (phi < -M_PI || phi > M_PI) {
423 while (phi < -M_PI) phi += 2 * M_PI;
424 while (phi > M_PI) phi -= 2 * M_PI;
425 result = true;
426 }
427
428 return result;
429 }
430 }// end OrcaKinFit namespace
432} // end Belle2 namespace
433
434
bool cachevalid
flag for valid cache
virtual int getGlobalParNum(int ilocal) const
Get global parameter number of parameter ilocal.
virtual bool isParamFixed(int ilocal) const
Returns whether parameter is fixed.
double par[BaseDefs::MAXPAR]
fit parameters
void invalidateCache() const
invalidate any cached quantities
double cov[BaseDefs::MAXPAR][BaseDefs::MAXPAR]
local covariance matrix
Class for jets with (E, eta, phi) in kinematic fits.
Definition: JetFitObject.h:43
virtual double getDPx(int ilocal) const override
Return d p_x / d par_ilocal (derivative of px w.r.t. local parameter ilocal)
virtual JetFitObject & assign(const BaseFitObject &source) override
Assign from anther object, if of same type.
static bool adjustEThetaPhi(const double &m, double &E, double &theta, double &phi)
Adjust E, theta and phi such that E>=m, 0<=theta<=pi, -pi <= phi < pi; returns true if anything was c...
virtual double getDPy(int ilocal) const override
Return d p_y / d par_ilocal (derivative of py w.r.t. local parameter ilocal)
virtual double getDE(int ilocal) const override
Return d E / d par_ilocal (derivative of E w.r.t. local parameter ilocal)
JetFitObject & operator=(const JetFitObject &rhs)
Assignment.
Definition: JetFitObject.cc:90
virtual JetFitObject * copy() const override
Return a new copy of itself.
Definition: JetFitObject.cc:98
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 double getFirstDerivative_Meta_Local(int iMeta, int ilocal, int metaSet) const override
add derivatives to vector der of size idim pxfact*dpx/dx_i + pyfact*dpy/dx_i + pzfact*dpz/dx_i + efac...
virtual bool updateParams(double p[], int idim) override
Read values from global vector, readjust vector; return: significant change.
virtual double getCov(int ilocal, int jlocal) const override
Get covariance between parameters ilocal and jlocal.
virtual double getError(int ilocal) const override
Get error of parameter ilocal.
virtual ParticleFitObject & assign(const BaseFitObject &source) override
Assign from anther object, if of same type.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24