Belle II Software light-2406-ragdoll
NeutrinoFitObject.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/NeutrinoFitObject.h"
18#include <framework/logging/Logger.h>
19#include <cmath>
20
21#undef NDEBUG
22#include <cassert>
23
24#include <algorithm>
25
26using std::sqrt;
27using std::sin;
28using std::cos;
29
30namespace Belle2 {
35 namespace OrcaKinFit {
36
37// constructor
38 NeutrinoFitObject::NeutrinoFitObject(double E, double theta, double phi,
39 double DE, double Dtheta, double Dphi)
40 : ctheta(0), stheta(0), cphi(0), sphi(0), pt(0), px(0), py(0), pz(0), dptdE(0),
41 dpxdE(0), dpydE(0), dpxdtheta(0), dpydtheta(0), chi2(0)
42 {
43
44 assert(int(NPAR) <= int(BaseDefs::MAXPAR));
45
46 setMass(0);
47 setParam(0, E, false);
48 setParam(1, theta, false);
49 setParam(2, phi, false);
50 setError(0, DE);
51 setError(1, Dtheta);
52 setError(2, Dphi);
53 invalidateCache();
54 }
55
56// destructor
57 NeutrinoFitObject::~NeutrinoFitObject() = default;
58
59 NeutrinoFitObject::NeutrinoFitObject(const NeutrinoFitObject& rhs)
60 : ParticleFitObject(rhs), ctheta(0), stheta(0), cphi(0), sphi(0), pt(0), px(0), py(0), pz(0), dptdE(0),
61 dpxdE(0), dpydE(0), dpxdtheta(0), dpydtheta(0), chi2(0)
62 {
64 }
65
67 {
68 if (this != &rhs) {
69 assign(rhs);
70 }
71 return *this;
72 }
73
75 {
76 return new NeutrinoFitObject(*this);
77 }
78
80 {
81 if (const auto* psource = dynamic_cast<const NeutrinoFitObject*>(&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* NeutrinoFitObject::getParamName(int ilocal) const
93 {
94 switch (ilocal) {
95 case 0: return "E";
96 case 1: return "theta";
97 case 2: return "phi";
98 }
99 return "undefined";
100 }
101
102
103 bool NeutrinoFitObject::updateParams(double pp[], int idim)
104 {
105
107
108 int iE = getGlobalParNum(0);
109 int ith = getGlobalParNum(1);
110 int iph = getGlobalParNum(2);
111 assert(iE >= 0 && iE < idim);
112 assert(ith >= 0 && ith < idim);
113 assert(iph >= 0 && iph < idim);
114
115 double e = pp[iE];
116 double th = pp[ith];
117 double ph = pp[iph];
118 if (e < 0) {
119 e = -e;
120 th = M_PI - th;
121 ph = M_PI + ph;
122 }
123 if (th < 0 || th > M_PI) {
124 th = M_PI - th;
125 ph = M_PI + ph;
126 }
127
128 bool result = (e - par[0]) * (e - par[0]) > eps2 * cov[0][0] ||
129 (th - par[1]) * (th - par[1]) > eps2 * cov[1][1] ||
130 (ph - par[2]) * (ph - par[2]) > eps2 * cov[2][2];
131 par[0] = e;
132 par[1] = (th >= 0 && th < M_PI) ?
133 th : std::acos(std::cos(th));
134 if (std::abs(ph) > M_PI) ph = atan2(sin(ph), cos(ph));
135 par[2] = ph;
136 pp[iE] = par[0];
137 pp[ith] = par[1];
138 pp[iph] = par[2];
139 return result;
140 }
141
142// these depend on actual parametrisation!
143 double NeutrinoFitObject::getDPx(int ilocal) const
144 {
145 assert(ilocal >= 0 && ilocal < NPAR);
146 if (!cachevalid) updateCache();
147 switch (ilocal) {
148 case 0: return dpxdE;
149 case 1: return dpxdtheta;
150 case 2: return -py;
151 }
152 return 0;
153 }
154
155 double NeutrinoFitObject::getDPy(int ilocal) const
156 {
157 assert(ilocal >= 0 && ilocal < NPAR);
158 if (!cachevalid) updateCache();
159 switch (ilocal) {
160 case 0: return dpydE;
161 case 1: return dpydtheta;
162 case 2: return px;
163 }
164 return 0;
165 }
166
167 double NeutrinoFitObject::getDPz(int ilocal) const
168 {
169 assert(ilocal >= 0 && ilocal < NPAR);
170 if (!cachevalid) updateCache();
171 switch (ilocal) {
172 case 0: return ctheta;
173 case 1: return -pt;
174 case 2: return 0;
175 }
176 return 0;
177 }
178
179 double NeutrinoFitObject::getDE(int ilocal) const
180 {
181 assert(ilocal >= 0 && ilocal < NPAR);
182 switch (ilocal) {
183 case 0: return 1;
184 case 1: return 0;
185 case 2: return 0;
186 }
187 return 0;
188 }
189
190 double NeutrinoFitObject::getFirstDerivative_Meta_Local(int iMeta, int ilocal, int metaSet) const
191 {
192 // iMeta = intermediate variable (i.e. E,px,py,pz)
193 // ilocal = local variable (E, theta, phi)
194 // metaSet = which set of intermediate varlables
195
196 assert(metaSet == 0); // only defined for E,px,py,pz
197
198 switch (iMeta) {
199 case 0: // E
200 return getDE(ilocal);
201 break;
202 case 1: // Px
203 return getDPx(ilocal);
204 break;
205 case 2: // Py
206 return getDPy(ilocal);
207 break;
208 case 3: // Pz
209 return getDPz(ilocal);
210 break;
211 default:
212 assert(0);
213 }
214 return -999;
215 }
216
217
218 double NeutrinoFitObject::getSecondDerivative_Meta_Local(int iMeta, int ilocal, int jlocal, int metaSet) const
219 {
220 assert(metaSet == 0);
221 if (!cachevalid) updateCache();
222
223 if (jlocal < ilocal) {
224 int temp = jlocal;
225 jlocal = ilocal;
226 ilocal = temp;
227 }
228
229 // daniel hasn't checked these, copied from orig code
230 switch (iMeta) {
231
232 case 0:
233 return 0;
234 break;
235 case 1:
236 if (ilocal == 0 && jlocal == 1) return ctheta * cphi;
237 else if (ilocal == 0 && jlocal == 2) return -dpydE;
238 else if (ilocal == 1 && jlocal == 1) return -px;
239 else if (ilocal == 1 && jlocal == 2) return -dpydtheta;
240 else if (ilocal == 2 && jlocal == 2) return -px;
241 else return 0;
242 break;
243 case 2:
244 if (ilocal == 0 && jlocal == 1) return ctheta * sphi;
245 else if (ilocal == 0 && jlocal == 2) return dpxdE;
246 else if (ilocal == 1 && jlocal == 1) return -py;
247 else if (ilocal == 1 && jlocal == 2) return dpxdtheta;
248 else if (ilocal == 2 && jlocal == 2) return -py;
249 else return 0;
250 break;
251 case 3:
252 if (ilocal == 0 && jlocal == 1) return -stheta;
253 else if (ilocal == 1 && jlocal == 1) return -pz;
254 else return 0;
255 break;
256 default:
257 assert(0);
258 }
259 return -999;
260 }
261
262
263 void NeutrinoFitObject::updateCache() const
264 {
265 double e = par[0];
266 double theta = par[1];
267 double phi = par[2];
268
269 ctheta = cos(theta);
270 stheta = sin(theta);
271 cphi = cos(phi);
272 sphi = sin(phi);
273
274 pt = e * stheta;
275
276 px = pt * cphi;
277 py = pt * sphi;
278 pz = e * ctheta;
279
280 fourMomentum.setValues(e, px, py, pz);
281
282 dpxdE = stheta * cphi;
283 dpydE = stheta * sphi;
284 dpxdtheta = pz * cphi;
285 dpydtheta = pz * sphi;
286
287 cachevalid = true;
288 }
289
290 }// end OrcaKinFit namespace
292} // end Belle2 namespace
293
bool cachevalid
flag for valid cache
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 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)
NeutrinoFitObject & operator=(const NeutrinoFitObject &rhs)
Assignment.
virtual double getDE(int ilocal) const override
Return d E / d par_ilocal (derivative of E w.r.t. local parameter ilocal)
virtual NeutrinoFitObject & assign(const BaseFitObject &source) override
Assign from anther object, if of same type.
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.
virtual NeutrinoFitObject * copy() const override
Return a new copy of itself.
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