17 #include "analysis/OrcaKinFit/JetFitObject.h"
18 #include <framework/logging/Logger.h>
36 namespace OrcaKinFit {
39 JetFitObject::JetFitObject(
double E,
double theta,
double phi,
40 double DE,
double Dtheta,
double Dphi,
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)
47 assert(
int(NPAR) <=
int(BaseDefs::MAXPAR));
58 adjustEThetaPhi(m, E, theta, phi);
60 setParam(1, theta,
true);
61 setParam(2, phi,
true);
70 paramCycl[2] = 2.*M_PI;
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]);
80 JetFitObject::~JetFitObject() =
default;
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)
105 if (
const auto* psource =
dynamic_cast<const JetFitObject*
>(&source)) {
106 if (psource !=
this) {
120 case 1:
return "theta";
121 case 2:
return "phi";
131 int iE = getGlobalParNum(0);
132 int ith = getGlobalParNum(1);
133 int iph = getGlobalParNum(2);
134 assert(iE >= 0 && iE < idim);
135 assert(ith >= 0 && ith < idim);
136 assert(iph >= 0 && iph < idim);
143 B2INFO(
"JetFitObject::updateParams: mirrored E!\n");
149 double massPlusEpsilon =
mass * (1.0000001);
150 if (e < massPlusEpsilon) e = massPlusEpsilon;
152 bool result = ((e - par[0]) * (e - par[0]) > eps2 * cov[0][0]) ||
153 ((th - par[1]) * (th - par[1]) > eps2 * cov[1][1]) ||
154 ((ph - par[2]) * (ph - par[2]) > eps2 * cov[2][2]);
167 assert(ilocal >= 0 && ilocal < NPAR);
168 if (!cachevalid) updateCache();
170 case 0:
return dpxdE;
171 case 1:
return dpxdtheta;
179 assert(ilocal >= 0 && ilocal < NPAR);
180 if (!cachevalid) updateCache();
182 case 0:
return dpydE;
183 case 1:
return dpydtheta;
191 assert(ilocal >= 0 && ilocal < NPAR);
192 if (!cachevalid) updateCache();
194 case 0:
return dpzdE;
203 assert(ilocal >= 0 && ilocal < NPAR);
214 assert(ilocal >= 0 && ilocal < NPAR);
215 B2INFO(
"JetFitObject::getError (ilocal = " << ilocal <<
") = " << std::sqrt(cov[ilocal][ilocal]));
216 return std::sqrt(cov[ilocal][ilocal]);
221 assert(ilocal >= 0 && ilocal < NPAR);
222 assert(jlocal >= 0 && jlocal < NPAR);
223 return cov[ilocal][jlocal];
253 assert(metaSet == 0);
256 return getDE(ilocal);
275 double JetFitObject::getSecondDerivative_Meta_Local(
int iMeta,
int ilocal ,
int jlocal ,
int metaSet)
const
277 assert(metaSet == 0);
278 if (!cachevalid) updateCache();
280 if (jlocal < ilocal) {
286 double d2pdE2 = (
mass != 0) ? -
mass *
mass / (p * p * p) : 0;
287 double d2ptdE2 = d2pdE2 * stheta;
296 if (ilocal == 0 && jlocal == 0)
return d2ptdE2 * cphi;
297 else if (ilocal == 0 && jlocal == 1)
return dpzdE * cphi;
298 else if (ilocal == 0 && jlocal == 2)
return -dpydE;
299 else if (ilocal == 1 && jlocal == 1)
return -px;
300 else if (ilocal == 1 && jlocal == 2)
return -dpydtheta;
301 else if (ilocal == 2 && jlocal == 2)
return -px;
304 if (ilocal == 0 && jlocal == 0)
return d2ptdE2 * sphi;
305 else if (ilocal == 0 && jlocal == 1)
return dpzdE * sphi;
306 else if (ilocal == 0 && jlocal == 2)
return dpxdE;
307 else if (ilocal == 1 && jlocal == 1)
return -py;
308 else if (ilocal == 1 && jlocal == 2)
return dpxdtheta;
309 else if (ilocal == 2 && jlocal == 2)
return -py;
312 if (ilocal == 0 && jlocal == 0)
return d2pdE2 * ctheta;
314 else if (ilocal == 0 && jlocal == 1)
return -dptdE;
315 else if (ilocal == 0 && jlocal == 2)
return 0;
316 else if (ilocal == 1 && jlocal == 1)
return -pz;
317 else if (ilocal == 1 && jlocal == 2)
return 0;
318 else if (ilocal == 2 && jlocal == 2)
return 0;
328 void JetFitObject::updateCache()
const
332 double theta = par[1];
350 fourMomentum.setValues(e, px, py, pz);
353 dptdE = dpdE * stheta;
354 dpxdE = dptdE * cphi;
355 dpydE = dptdE * sphi;
356 dpzdE = dpdE * ctheta;
358 dpxdtheta = pz * cphi;
359 dpydtheta = pz * sphi;
376 B2INFO(
"JetFitObject::adjustEThetaPhi: mirrored E!\n");
378 theta = M_PI - theta;
386 if (theta < -M_PI || theta > M_PI) {
387 while (theta < -M_PI) theta += 2 * M_PI;
388 while (theta > M_PI) theta -= 2 * M_PI;
393 B2INFO(
"JetFitObject::adjustEThetaPhi: mirrored theta!\n");
395 phi = phi > 0 ? phi - M_PI : phi + M_PI;
397 }
else if (theta > M_PI) {
398 B2INFO(
"JetFitObject::adjustEThetaPhi: mirrored theta!\n");
399 theta = 2 * M_PI - theta;
400 phi = phi > 0 ? phi - M_PI : phi + M_PI;
403 if (phi < -M_PI || phi > M_PI) {
404 while (phi < -M_PI) phi += 2 * M_PI;
405 while (phi > M_PI) phi -= 2 * M_PI;