Belle II Software light-2406-ragdoll
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/FitParams.h>
11#include <analysis/VertexFitting/TreeFitter/ParticleBase.h>
12#include <analysis/VertexFitting/TreeFitter/Constraint.h>
13#include <analysis/VertexFitting/TreeFitter/KalmanCalculator.h>
14
15namespace TreeFitter {
16
17 bool Constraint::operator<(const Constraint& rhs) const
18 {
19 return m_depth < rhs.m_depth ||
20 (m_depth == rhs.m_depth && m_type < rhs.m_type);
21 }
22
24 {
25 return m_node->projectConstraint(m_type, fitpar, p);
26 }
27
29 {
34 ErrCode status;
37
38 double chisq(0);
39 int iter(0);
40 bool finished(false) ;
41
42 double accumulated_chi2 = 0;
43 while (!finished && !status.failure()) {
44
45 p.resetProjection();
46 status |= project(fitpar, p);
47
48 if (!status.failure()) {
49
50 status |= kalman.calculateGainMatrix(
51 p.getResiduals(),
52 p.getH(),
53 fitpar,
54 &p.getV(),
55 1
56 );
57
58 if (!status.failure()) {
59 kalman.updateState(fitpar);
60
61 // r R^-1 r
62 double newchisq = kalman.getChiSquare();
63
64 double dchisqconverged = 0.001 ;
65
66 double dchisq = newchisq - chisq;
67 bool diverging = iter > 0 && dchisq > 0;
68 bool converged = std::abs(dchisq) < dchisqconverged;
69 finished = ++iter >= m_maxNIter || diverging || converged;
70 chisq = newchisq;
71 accumulated_chi2 += newchisq;
72 }
73 }
74 }
75
76 const unsigned int number_of_constraints = kalman.getConstraintDim();
77 fitpar.addChiSquare(accumulated_chi2, number_of_constraints);
78
79 kalman.updateCovariance(fitpar);
80 return status;
81 }
82
84 {
91 ErrCode status;
94
95 p.resetProjection();
96 status |= project(oldState, p);
97
102 p.getResiduals() += p.getH() * (fitpar.getStateVector() - oldState.getStateVector());
103 if (!status.failure()) {
104 status |= kalman.calculateGainMatrix(
105 p.getResiduals(),
106 p.getH(),
107 fitpar,
108 &p.getV(),
109 1
110 );
111
112 if (!status.failure()) {
113 kalman.updateState(fitpar);
114 }
115 }
116
117 const unsigned int number_of_constraints = kalman.getConstraintDim();
118 fitpar.addChiSquare(kalman.getChiSquare(), number_of_constraints);
119
120 kalman.updateCovariance(fitpar);
121 return status;
122 }
123
124 std::string Constraint::name() const
125 {
126 std::string rc = "unknown constraint!";
127 switch (m_type) {
128 case beamspot: rc = "beamspot"; break;
129 case beamenergy: rc = "beamenergy"; break;
130 case beam: rc = "beam"; break;
131 case origin: rc = "origin"; break;
132 case composite: rc = "composite"; break;
133 case resonance: rc = "resonance"; break;
134 case track: rc = "track"; break;
135 case photon: rc = "photon"; break;
136 case klong: rc = "klong"; 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}
class to manage the order of constraints and their filtering
Definition: Constraint.h:20
Type m_type
type of constraint
Definition: Constraint.h:140
std::string name() const
get name of constraint
Definition: Constraint.cc:124
int m_depth
chi2 coming from the constraint
Definition: Constraint.h:137
unsigned int m_dim
dimension of constraint
Definition: Constraint.h:143
int m_maxNIter
maximum number of iterations for non-linear constraints
Definition: Constraint.h:149
virtual ErrCode filterWithReference(FitParams &fitpar, const FitParams &oldState)
filter this constraint
Definition: Constraint.cc:83
bool operator<(const Constraint &rhs) const
operator used to sort the constraints
Definition: Constraint.cc:17
const ParticleBase * m_node
particle behind the constraint
Definition: Constraint.h:131
virtual ErrCode filter(FitParams &fitpar)
filter this constraint
Definition: Constraint.cc:28
virtual ErrCode project(const FitParams &fitpar, Projection &p) const
call the constraints projection function
Definition: Constraint.cc:23
abstract errorocode be aware that the default is success
Definition: ErrCode.h:14
bool failure() const
returns true if errorcode is error
Definition: ErrCode.h:58
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
virtual ErrCode projectConstraint(Constraint::Type, const FitParams &, Projection &) const
project constraint.
class to store the projected residuals and the corresponding jacobian as well as the covariance matri...
Definition: Projection.h:18