Belle II Software light-2406-ragdoll
ParticleFitObject.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/ParticleFitObject.h"
18#include <framework/logging/Logger.h>
19
20#include <iostream>
21#undef NDEBUG
22#include <cassert>
23#include <cmath>
24using std::isfinite;
25using std::endl;
26
27// #include <TMatrixDSym.h>
28#include <gsl/gsl_matrix.h>
29#include <gsl/gsl_linalg.h>
30
31
32namespace Belle2 {
37 namespace OrcaKinFit {
38
40 : mass(0), fourMomentum(FourVector(0, 0, 0, 0)), paramCycl{}
41 {
42 for (double& i : paramCycl)
43 i = -1;
44 }
45
47
49 : BaseFitObject(rhs), mass(0), fourMomentum(FourVector(0, 0, 0, 0)), paramCycl{}
50 {
51 //B2INFO( "copying ParticleFitObject with name" << rhs.name );
53 }
54
56 {
57 if (this != &rhs) {
58 assign(rhs); // calls virtual function assign of derived class
59 }
60 return *this;
61 }
62
63 bool ParticleFitObject::setMass(double mass_)
64 {
65 if (!isfinite(mass_)) return false;
66 if (mass == mass_) return true;
68 mass = std::abs(mass_);
69 return true;
70 }
71
72
74 {
75 if (const auto* psource = dynamic_cast<const ParticleFitObject*>(&source)) {
76 if (psource != this) {
78 mass = psource->mass;
79 for (int i = 0; i < BaseDefs::MAXPAR; ++i)
80 paramCycl[i] = psource->paramCycl[i];
81 }
82 } else {
83 assert(0);
84 }
85 return *this;
86 }
87
88
90 {
91 return mass;
92 }
93
94 std::ostream& ParticleFitObject::print4Vector(std::ostream& os) const
95 {
96 os << "[" << getE() << ", " << getPx() << ", "
97 << getPy() << ", " << getPz() << "]";
98 return os;
99 }
100
101 FourVector ParticleFitObject::getFourMomentum() const
102 {
103 if (!cachevalid) updateCache();
104 return fourMomentum;
105 }
107 {
108 return getFourMomentum().getE();
109 }
111 {
112 return getFourMomentum().getPx();
113 }
115 {
116 return getFourMomentum().getPy();
117 }
119 {
120 return getFourMomentum().getPz();
121 }
123 {
124 return getFourMomentum().getP();
125 }
127 {
128 return getFourMomentum().getP2();
129 }
131 {
132 return getFourMomentum().getPt();
133 }
135 {
136 return getFourMomentum().getPt2();
137 }
138
139
140 std::ostream& ParticleFitObject::print(std::ostream& os) const
141 {
142
143 if (!cachevalid) updateCache();
144
145 printParams(os);
146 os << " => ";
147 print4Vector(os);
148 os << std::endl;
149 return os;
150 }
151
152 void ParticleFitObject::addToGlobalChi2DerVectorNum(double* y, int idim, double eps)
153 {
154 for (int ilocal = 0; ilocal < getNPar(); ++ilocal) {
155 int iglobal = getGlobalParNum(ilocal);
156 assert(iglobal >= 0 && iglobal < idim);
157 y[iglobal] += num1stDerivative(ilocal, eps);
158 }
159 }
160
161
162 void ParticleFitObject::addToGlobalChi2DerMatrixNum(double* M, int idim, double eps)
163 {
164 for (int ilocal1 = 0; ilocal1 < getNPar(); ++ilocal1) {
165 int iglobal1 = getGlobalParNum(ilocal1);
166 for (int ilocal2 = ilocal1; ilocal2 < getNPar(); ++ilocal2) {
167 int iglobal2 = getGlobalParNum(ilocal2);
168 M[idim * iglobal1 + iglobal2] += num2ndDerivative(ilocal1, eps, ilocal2, eps);
169 }
170 }
171 }
172
173 void ParticleFitObject::getDerivatives(double der[], int idim) const
174 {
175 for (int ilocal = 0; ilocal < getNPar(); ++ilocal) {
176 assert(ilocal < idim);
177 der [4 * ilocal] = getDE(ilocal);
178 der [4 * ilocal + 1] = getDPx(ilocal);
179 der [4 * ilocal + 2] = getDPy(ilocal);
180 der [4 * ilocal + 3] = getDPz(ilocal);
181 }
182 }
183
184 void ParticleFitObject::test1stDerivatives()
185 {
186 B2INFO("ParticleFitObject::test1stDerivatives, object " << getName() << "\n");
187 double ycalc[100], ynum[100];
188 for (int i = 0; i < 100; ++i) ycalc[i] = ynum[i] = 0;
189 addToGlobalChi2DerVector(ycalc, 100);
190 double eps = 0.00001;
191 addToGlobalChi2DerVectorNum(ynum, 100, eps);
192 for (int ilocal = 0; ilocal < getNPar(); ++ilocal) {
193 int iglobal = getGlobalParNum(ilocal);
194 double calc = ycalc[iglobal];
195 double num = ynum[iglobal];
196 B2INFO("fo: " << getName() << " par " << ilocal << "/"
197 << iglobal << " (" << getParamName(ilocal)
198 << ") calc: " << calc << " - num: " << num << " = " << calc - num);
199 }
200 }
201
202 void ParticleFitObject::test2ndDerivatives()
203 {
204 B2INFO("ParticleFitObject::test2ndDerivatives, object " << getName() << "\n");
205 const int idim = 100;
206 auto* Mnum = new double[idim * idim];
207 auto* Mcalc = new double[idim * idim];
208 for (int i = 0; i < idim * idim; ++i) Mnum[i] = Mcalc[i] = 0;
209 addToGlobalChi2DerMatrix(Mcalc, idim);
210 double eps = 0.0001;
211 B2INFO("eps=" << eps);
212 addToGlobalChi2DerMatrixNum(Mnum, idim, eps);
213 for (int ilocal1 = 0; ilocal1 < getNPar(); ++ilocal1) {
214 int iglobal1 = getGlobalParNum(ilocal1);
215 for (int ilocal2 = ilocal1; ilocal2 < getNPar(); ++ilocal2) {
216 int iglobal2 = getGlobalParNum(ilocal2);
217 double calc = Mcalc[idim * iglobal1 + iglobal2];
218 double num = Mnum[idim * iglobal1 + iglobal2];
219 B2INFO("fo: " << getName() << " par " << ilocal1 << "/"
220 << iglobal1 << " (" << getParamName(ilocal1)
221 << "), par " << ilocal2 << "/"
222 << iglobal2 << " (" << getParamName(ilocal2)
223 << ") calc: " << calc << " - num: " << num << " = " << calc - num);
224 }
225 }
226 delete[] Mnum;
227 delete[] Mcalc;
228 }
229
230 double ParticleFitObject::num1stDerivative(int ilocal, double eps)
231 {
232 double save = getParam(ilocal);
233 setParam(ilocal, save + eps);
234 double v1 = getChi2();
235 setParam(ilocal, save - eps);
236 double v2 = getChi2();
237 double result = (v1 - v2) / (2 * eps);
238 setParam(ilocal, save);
239 return result;
240 }
241
242 double ParticleFitObject::num2ndDerivative(int ilocal1, double eeps1,
243 int ilocal2, double eeps2)
244 {
245 double result;
246
247 if (ilocal1 == ilocal2) {
248 double save = getParam(ilocal1);
249 double v0 = getChi2();
250 setParam(ilocal1, save + eeps1);
251 double v1 = getChi2();
252 setParam(ilocal1, save - eeps1);
253 double v2 = getChi2();
254 result = (v1 + v2 - 2 * v0) / (eeps1 * eeps1);
255 setParam(ilocal1, save);
256 } else {
257 double save1 = getParam(ilocal1);
258 double save2 = getParam(ilocal2);
259 setParam(ilocal1, save1 + eeps1);
260 setParam(ilocal2, save2 + eeps2);
261 double v11 = getChi2();
262 setParam(ilocal2, save2 - eeps2);
263 double v12 = getChi2();
264 setParam(ilocal1, save1 - eeps1);
265 double v22 = getChi2();
266 setParam(ilocal2, save2 + eeps2);
267 double v21 = getChi2();
268 result = (v11 + v22 - v12 - v21) / (4 * eeps1 * eeps2);
269 setParam(ilocal1, save1);
270 setParam(ilocal2, save2);
271 }
272 return result;
273 }
274
275
277 {
278 // reimplemented here to take account of cyclical variables e.g azimuthal angle phi - DJeans
279
280 if (not covinvvalid and not calculateCovInv()) return 0;
281
282 double resid[BaseDefs::MAXPAR] = {0};
283 for (int i = 0; i < getNPar(); i++) {
284 if (isParamMeasured(i) && !isParamFixed(i)) {
285 resid[i] = par[i] - mpar[i];
286 if (paramCycl[i] > 0) {
287 resid[i] = fmod(resid[i], paramCycl[i]);
288 if (resid[i] > paramCycl[i] / 2) resid[i] -= paramCycl[i];
289 if (resid[i] < -paramCycl[i] / 2) resid[i] += paramCycl[i];
290 }
291 }
292 }
293
294
295 double chi2 = 0;
296 for (int i = 0; i < getNPar(); i++) {
297 if (isParamMeasured(i) && !isParamFixed(i)) {
298 for (int j = 0; j < getNPar(); j++) {
299 if (isParamMeasured(j) && !isParamFixed(j)) {
300 chi2 += resid[i] * covinv[i][j] * resid[j];
301 // B2INFO( getName () << " === " << i << " " << j << " : " <<
302 // resid[i] << "*" << covinv[i][j] << "*" << resid[j] << " = " << resid[i]*covinv[i][j]*resid[j] << " , sum " << chi2);
303 }
304 }
305 }
306 }
307
308 // B2INFO("ParticleFitObject::getChi2 () chi2 = " << chi2 );
309 return chi2;
310 }
311
312 }// end OrcaKinFit namespace
314} // end Belle2 namespace
bool covinvvalid
flag for valid inverse covariance matrix
virtual int getNPar() const =0
Get total number of parameters of this FitObject.
virtual void addToGlobalChi2DerMatrix(double *M, int idim) const
Add 2nd derivatives of chi squared to global derivative matrix.
virtual const char * getParamName(int ilocal) const =0
Get name of parameter ilocal.
virtual BaseFitObject & assign(const BaseFitObject &source)
Assign from anther object, if of same type.
virtual bool setParam(int ilocal, double par_, bool measured_, bool fixed_=false)
Set value and measured flag of parameter i; return: significant change.
double covinv[BaseDefs::MAXPAR][BaseDefs::MAXPAR]
inverse pf local covariance matrix
bool cachevalid
flag for valid cache
virtual int getGlobalParNum(int ilocal) const
Get global parameter number of parameter ilocal.
virtual const char * getName() const
Get object's name.
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 mpar[BaseDefs::MAXPAR]
measured parameters
virtual double getParam(int ilocal) const
Get current value of parameter ilocal.
virtual int addToGlobalChi2DerVector(double *y, int idim) const
Add derivatives of chi squared to global derivative vector.
virtual std::ostream & printParams(std::ostream &os) const
print the parameters and errors
virtual bool calculateCovInv() const
Calculate the inverse of the covariance matrix.
virtual bool isParamMeasured(int ilocal) const
Get measured flag for parameter ilocal.
Yet another four vector class, with metric +—.
Definition: FourVector.h:43
double getPx() const
Returns the x momentum / 1 component.
Definition: FourVector.h:128
double getPz() const
Returns the z momentum / 3 component.
Definition: FourVector.h:130
double getPt2() const
Returns the transverse momentum squared / magnitude of the 1 and 2 component vector squared.
Definition: FourVector.h:132
double getPy() const
Returns the y momentum / 2 component.
Definition: FourVector.h:129
double getP2() const
Returns the momentum squared / magnitude of the three vector squared.
Definition: FourVector.h:135
double getPt() const
Returns the transverse momentum / magnitude of the 1 and 2 component vector.
Definition: FourVector.h:133
double getP() const
Returns the momentum / magnitude of the three vector.
Definition: FourVector.h:136
double getE() const
Returns the energy / 0 component.
Definition: FourVector.h:127
virtual double getPx() const
Return px.
virtual ~ParticleFitObject()
Virtual destructor.
virtual void addToGlobalChi2DerMatrixNum(double *M, int idim, double eps)
Add numerically determined derivatives of chi squared to global covariance matrix.
ParticleFitObject & operator=(const ParticleFitObject &rhs)
Assignment.
virtual double getPz() const
Return pz.
virtual bool setMass(double mass_)
Set mass of particle; return=success.
double num2ndDerivative(int ilocal1, double eps1, int ilocal2, double eps2)
Evaluates numerically the 2nd derivative of chi2 w.r.t. 2 parameters.
virtual void addToGlobalChi2DerVectorNum(double *y, int idim, double eps)
Add numerically determined derivatives of chi squared to global derivative vector.
virtual std::ostream & print(std::ostream &os) const override
print object to ostream
virtual double getPt2() const
Return pt (transverse momentum) squared.
virtual double getDE(int ilocal) const =0
Return d E / d par_ilocal (derivative of E w.r.t. local parameter ilocal)
virtual double getPy() const
Return py.
virtual double getDPx(int ilocal) const =0
Return d p_x / d par_ilocal (derivative of px w.r.t. local parameter ilocal)
virtual double getDPz(int ilocal) const =0
Return d p_z / d par_ilocal (derivative of pz w.r.t. local parameter ilocal)
virtual double getP2() const
Return p (momentum) squared.
virtual double getPt() const
Return pt (transverse momentum)
virtual double getP() const
Return p (momentum)
virtual ParticleFitObject & assign(const BaseFitObject &source) override
Assign from anther object, if of same type.
virtual double getE() const
Return E.
double num1stDerivative(int ilocal, double eps)
Evaluates numerically the 1st derivative of chi2 w.r.t. a parameter.
virtual double getChi2() const override
Get chi squared from measured and fitted parameters.
virtual double getDPy(int ilocal) const =0
Return d p_y / d par_ilocal (derivative of py w.r.t. local parameter ilocal)
virtual double getMass() const
Get mass of particle.
virtual std::ostream & print4Vector(std::ostream &os) const
print the four-momentum (E, px, py, pz)
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24