Belle II Software development
Constraint.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * External Contributor: Wouter Hulsbergen *
5 * *
6 * See git log for contributors and copyright holders. *
7 * This file is licensed under LGPL-3.0, see LICENSE.md. *
8 **************************************************************************/
9
10#include <analysis/VertexFitting/TreeFitter/Constraint.h>
11#include <analysis/VertexFitting/TreeFitter/FitParams.h>
12#include <analysis/VertexFitting/TreeFitter/KalmanCalculator.h>
13#include <analysis/VertexFitting/TreeFitter/ParticleBase.h>
14#include <analysis/VertexFitting/TreeFitter/Projection.h>
15
16namespace TreeFitter {
17
18 bool Constraint::operator<(const Constraint& rhs) const
19 {
20 return m_depth < rhs.m_depth ||
21 (m_depth == rhs.m_depth && m_type < rhs.m_type);
22 }
23
25 {
26 return m_node->projectConstraint(m_type, fitpar, p);
27 }
28
30 {
35 ErrCode status;
38
39 double chisq(0);
40 int iter(0);
41 bool finished(false) ;
42
43 double accumulated_chi2 = 0;
44 while (!finished && !status.failure()) {
45
46 p.resetProjection();
47 status |= project(fitpar, p);
48
49 if (!status.failure()) {
50
51 status |= kalman.calculateGainMatrix(
52 p.getResiduals(),
53 p.getH(),
54 fitpar,
55 &p.getV(),
56 1
57 );
58
59 if (!status.failure()) {
60 kalman.updateState(fitpar);
61
62 // r R^-1 r
63 double newchisq = kalman.getChiSquare();
64
65 double dchisqconverged = 0.001 ;
66
67 double dchisq = newchisq - chisq;
68 bool diverging = iter > 0 && dchisq > 0;
69 bool converged = std::abs(dchisq) < dchisqconverged;
70 finished = ++iter >= m_maxNIter || diverging || converged;
71 chisq = newchisq;
72 accumulated_chi2 += newchisq;
73 }
74 }
75 }
76
77 const unsigned int number_of_constraints = kalman.getConstraintDim();
78 fitpar.addChiSquare(accumulated_chi2, number_of_constraints);
79
80 kalman.updateCovariance(fitpar);
81 return status;
82 }
83
85 {
92 ErrCode status;
95
96 p.resetProjection();
97 status |= project(oldState, p);
98
103 p.getResiduals() += p.getH() * (fitpar.getStateVector() - oldState.getStateVector());
104 if (!status.failure()) {
105 status |= kalman.calculateGainMatrix(
106 p.getResiduals(),
107 p.getH(),
108 fitpar,
109 &p.getV(),
110 1
111 );
112
113 if (!status.failure()) {
114 kalman.updateState(fitpar);
115 }
116 }
117
118 const unsigned int number_of_constraints = kalman.getConstraintDim();
119 fitpar.addChiSquare(kalman.getChiSquare(), number_of_constraints);
120
121 kalman.updateCovariance(fitpar);
122 return status;
123 }
124
125 std::string Constraint::name() const
126 {
127 std::string rc = "unknown constraint!";
128 switch (m_type) {
129 case beamspot: rc = "beamspot"; break;
130 case beamenergy: rc = "beamenergy"; break;
131 case beam: rc = "beam"; break;
132 case origin: rc = "origin"; break;
133 case composite: rc = "composite"; break;
134 case resonance: rc = "resonance"; break;
135 case track: rc = "track"; break;
136 case neutralHadron: rc = "neutralHadron"; break;
137 case kinematic: rc = "kinematic"; break;
138 case geometric: rc = "geometric"; break;
139 case mass: rc = "mass"; break;
140 case massEnergy: rc = "massEnergy"; break;
141 case lifetime: rc = "lifetime"; break;
142 case merged: rc = "merged"; break;
143 case conversion: rc = "conversion"; break;
144 case helix: rc = "helix"; break;
145 case ntypes:
146 case unknown:
147 break;
148 }
149 return rc;
150 }
151}
Type m_type
type of constraint
Definition Constraint.h:139
std::string name() const
get name of constraint
Constraint()
constructor
Definition Constraint.h:69
int m_depth
chi2 coming from the constraint
Definition Constraint.h:136
unsigned int m_dim
dimension of constraint
Definition Constraint.h:142
int m_maxNIter
maximum number of iterations for non-linear constraints
Definition Constraint.h:148
virtual ErrCode filterWithReference(FitParams &fitpar, const FitParams &oldState)
filter this constraint
Definition Constraint.cc:84
bool operator<(const Constraint &rhs) const
operator used to sort the constraints
Definition Constraint.cc:18
const ParticleBase * m_node
particle behind the constraint
Definition Constraint.h:130
virtual ErrCode filter(FitParams &fitpar)
filter this constraint
Definition Constraint.cc:29
virtual ErrCode project(const FitParams &fitpar, Projection &p) const
call the constraints projection function
Definition Constraint.cc:24
abstract errorocode be aware that the default is success
Definition ErrCode.h:14
Class to store and manage fitparams (statevector)
Definition FitParams.h:20
Eigen::Matrix< double, -1, 1, 0, MAX_MATRIX_SIZE, 1 > & getStateVector()
getter for the fit parameters/statevector
Definition FitParams.h:65
void addChiSquare(double chisq, int nconstraints)
increment global chi2
Definition FitParams.h:92
int getDimensionOfState() const
get the states dimension
Definition FitParams.h:80
does the calculation of the gain matrix, updates the cov and fitpars
void updateState(FitParams &fitparams)
update statevector
ErrCode calculateGainMatrix(const Eigen::Matrix< double, -1, 1, 0, 7, 1 > &residuals, const Eigen::Matrix< double, -1, -1, 0, 7, MAX_MATRIX_SIZE > &G, const FitParams &fitparams, const Eigen::Matrix< double, -1, -1, 0, 7, 7 > *V=0, double weight=1)
init the kalman machienery
double getConstraintDim() const
get dimension of the constraint
void updateCovariance(FitParams &fitparams)
update the statevectors covariance
double getChiSquare() const
get chi2 of this iteration
class to store the projected residuals and the corresponding jacobian as well as the covariance matri...
Definition Projection.h:18