Belle II Software light-2406-ragdoll
FitManager.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#include <TMath.h>
10#include <Math/Vector4D.h>
11
12#include <analysis/dataobjects/Particle.h>
13
14#include <framework/logging/Logger.h>
15#include <framework/gearbox/Const.h>
16
17#include <analysis/VertexFitting/TreeFitter/FitManager.h>
18#include <analysis/VertexFitting/TreeFitter/FitParams.h>
19#include <analysis/VertexFitting/TreeFitter/DecayChain.h>
20#include <analysis/VertexFitting/TreeFitter/ParticleBase.h>
21
22namespace TreeFitter {
23
25 const ConstraintConfiguration& config,
26 double prec,
27 bool updateDaughters,
28 const bool useReferencing
29 ) :
30 m_particle(particle),
31 m_decaychain(nullptr),
32 m_status(VertexStatus::UnFitted),
33 m_chiSquare(-1),
34 m_prec(prec),
35 m_updateDaugthers(updateDaughters),
36 m_ndf(0),
37 m_fitparams(nullptr),
38 m_useReferencing(useReferencing),
39 m_config(config)
40 {
41 m_decaychain = new DecayChain(particle, config, false);
43 }
44
46 {
47 delete m_decaychain;
48 delete m_fitparams;
49 }
50
52 {
53 const int nitermax = 100;
54 const int maxndiverging = 3;
55 const double dChisqConv = m_prec;
56 m_chiSquare = -1;
58
59 if (m_status == VertexStatus::UnFitted) {
61 }
62
63 if (m_errCode.failure()) {
64 m_status = VertexStatus::BadInput;
65 } else {
66 m_status = VertexStatus::UnFitted;
67 int ndiverging = 0;
68 bool finished = false;
69 int niter = 0;
70 for (niter = 0; niter < nitermax && !finished; ++niter) {
71 if (niter == 0) {
73 } else if (m_useReferencing) {
74 auto* tempState = new FitParams(*m_fitparams);
76 delete tempState;
77 }
79 double chisq = m_fitparams->chiSquare();
80 double deltachisq = chisq - m_chiSquare;
81 if (m_errCode.failure()) {
82 finished = true ;
83 m_status = VertexStatus::Failed;
84 m_particle->writeExtraInfo("failed", 1);
85 } else {
86 if (niter > 0) {
87 if ((std::abs(deltachisq) / m_chiSquare < dChisqConv)) {
88 m_chiSquare = chisq;
89 m_status = VertexStatus::Success;
90 finished = true ;
91 m_particle->writeExtraInfo("failed", 0);
92 } else if (deltachisq > 0 && ++ndiverging >= maxndiverging) {
93 m_particle->writeExtraInfo("failed", 2);
94 m_status = VertexStatus::NonConverged;
95 m_errCode = ErrCode(ErrCode::Status::slowdivergingfit);
96 finished = true ;
97 }
98 }
99 if (deltachisq < 0) {
100 ndiverging = 0;
101 }
102 m_chiSquare = chisq;
103 }
104 }
105 if (niter == nitermax && m_status != VertexStatus::Success) {
106 m_particle->writeExtraInfo("failed", 3);
107 m_status = VertexStatus::NonConverged;
108 }
109 if (!(m_fitparams->testCovariance())) {
110 m_particle->writeExtraInfo("failed", 4);
111 m_status = VertexStatus::Failed;
112 }
113 }
114
115 if (m_status == VertexStatus::Success) {
116 // mass constraints comes after kine so we have to
117 // update the mothers with the values set by the mass constraint
118 if (m_config.m_massConstraintListPDG.size() != 0) {
120 }
121 updateTree(*m_particle, true);
122 }
123
124 return (m_status == VertexStatus::Success);
125 }
126
127 void FitManager::getCovFromPB(const ParticleBase* pb, TMatrixFSym& returncov) const
128 {
129
130 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> cov = m_fitparams->getCovariance().selfadjointView<Eigen::Lower>();
131 int posindex = pb->posIndex();
132 // hack: for tracks and photons, use the production vertex
133 if (posindex < 0 && pb->mother()) {
134 posindex = pb->mother()->posIndex();
135 }
136 int momindex = pb->momIndex();
137 if (pb->hasEnergy()) {
138 // if particle has energy, get full p4 from fitparams and put them directly in the return type
139 // very important! Belle2 uses p,E,x! Change order here!
140 for (int row = 0; row < 4; ++row) {
141 for (int col = 0; col < 4; ++col) {
142 returncov(row, col) = cov(momindex + row, momindex + col);
143 }
144 }
145
146 for (int row = 0; row < 3; ++row) {
147 for (int col = 0; col < 3; ++col) {
148 returncov(row + 4, col + 4) = cov(posindex + row, posindex + col);
149 }
150 }
151
152 } else {
153 Eigen::Matrix<double, 6, 6> cov6 =
154 Eigen::Matrix<double, 6, 6>::Zero(6, 6);
155
156 for (int row = 0; row < 3; ++row) {
157 for (int col = 0; col < 3; ++col) {
158 cov6(row, col) = cov(momindex + row, momindex + col);
159 cov6(row + 3, col + 3) = cov(posindex + row, posindex + col);
160 }
161 }
162
163 double mass = pb->particle()->getPDGMass();
164 Eigen::Matrix<double, 3, 1> momVec =
165 m_fitparams->getStateVector().segment(momindex, 3);
166
167 double energy2 = momVec.transpose() * momVec;
168 energy2 += mass * mass;
169 double energy = sqrt(energy2);
170
171 Eigen::Matrix<double, 7, 6> jacobian =
172 Eigen::Matrix<double, 7, 6>::Zero(7, 6);
173
174 for (int col = 0; col < 3; ++col) {
175 jacobian(col, col) = 1;
176 jacobian(3, col) = m_fitparams->getStateVector()(momindex + col) / energy;
177 jacobian(col + 4, col + 3) = 1;
178 }
179
180 Eigen::Matrix<double, 7, 7> cov7
181 = jacobian * cov6.selfadjointView<Eigen::Lower>() * jacobian.transpose();
182
183 for (int row = 0; row < 7; ++row) {
184 for (int col = 0; col < 7; ++col) {
185 returncov(row, col) = cov7(row, col);
186 }
187 }
188 } // else
189 }
190
191 bool FitManager::updateCand(Belle2::Particle& cand, const bool isTreeHead) const
192 {
193 const ParticleBase* pb = m_decaychain->locate(&cand);
194 if (pb) {
195 updateCand(*pb, cand, isTreeHead);
196 } else {
197 B2ERROR("Can't find candidate " << cand.getName() << " in tree " << m_particle->getName());
198 }
199 return pb != nullptr;
200 }
201
203 Belle2::Particle& cand, const bool isTreeHead) const
204 {
205 int posindex = pb.posIndex();
206 if (posindex < 0 && pb.mother()) {
207 posindex = pb.mother()->posIndex();
208 }
209
210 if (m_updateDaugthers || isTreeHead) {
211 if (posindex >= 0) {
212 const ROOT::Math::XYZVector pos(m_fitparams->getStateVector()(posindex),
213 m_fitparams->getStateVector()(posindex + 1),
214 m_fitparams->getStateVector()(posindex + 2));
215 cand.setVertex(pos);
216 if (&pb == m_decaychain->cand()) { // if head
217 const double fitparchi2 = m_fitparams->chiSquare();
218 cand.setPValue(TMath::Prob(fitparchi2, m_ndf));//if m_ndf<1, this is 0.
219 cand.writeExtraInfo("chiSquared", fitparchi2);
220 cand.writeExtraInfo("modifiedPValue", TMath::Prob(fitparchi2, 3));
221 cand.writeExtraInfo("ndf", m_ndf);
222 }
223 if (pb.mother()) {
224 int motherPosIndex = pb.mother()->posIndex();
225 if (motherPosIndex >= 0) {
226 cand.writeExtraInfo("prodVertexX", m_fitparams->getStateVector()(motherPosIndex));
227 cand.writeExtraInfo("prodVertexY", m_fitparams->getStateVector()(motherPosIndex + 1));
228 if (pb.mother()->dim() > 2)
229 cand.writeExtraInfo("prodVertexZ", m_fitparams->getStateVector()(motherPosIndex + 2));
230 }
231 }
232 }
233
234 const int momindex = pb.momIndex();
235 ROOT::Math::PxPyPzEVector p;
236 p.SetPx(m_fitparams->getStateVector()(momindex));
237 p.SetPy(m_fitparams->getStateVector()(momindex + 1));
238 p.SetPz(m_fitparams->getStateVector()(momindex + 2));
239 if (pb.hasEnergy()) {
240 p.SetE(m_fitparams->getStateVector()(momindex + 3));
242 } else {
243 const double mass = cand.getPDGMass();
244 p.SetE(std::sqrt(p.P2() + mass * mass));
246 }
247 TMatrixFSym cov7b2(7);
248 getCovFromPB(&pb, cov7b2);
249 cand.setMomentumVertexErrorMatrix(cov7b2);
250 }
251
252 if (pb.tauIndex() > 0) {
253 const std::tuple<double, double>tau = getDecayLength(cand);
254 const std::tuple<double, double>life = getLifeTime(cand);
255 cand.writeExtraInfo("decayLength", std::get<0>(tau));
256 cand.writeExtraInfo("decayLengthErr", std::get<1>(tau));
257 cand.writeExtraInfo("lifeTime", std::get<0>(life));
258 cand.writeExtraInfo("lifeTimeErr", std::get<1>(life));
259 }
260 }
261
262 void FitManager::updateTree(Belle2::Particle& cand, const bool isTreeHead) const
263 {
264 const bool updateableMother = updateCand(cand, isTreeHead);
265
266 if (updateableMother and not cand.hasExtraInfo("bremsCorrected") and
267 not(cand.hasExtraInfo("treeFitterTreatMeAsInvisible") and cand.getExtraInfo("treeFitterTreatMeAsInvisible") == 1)) {
268 const int ndaughters = cand.getNDaughters();
269 for (int i = 0; i < ndaughters; i++) {
270 auto* daughter = const_cast<Belle2::Particle*>(cand.getDaughter(i));
271 updateTree(*daughter, false);
272 }
273 }
274 }
275
276 std::tuple<double, double> FitManager::getLifeTime(Belle2::Particle& cand) const
277 {
278 const ParticleBase* pb = m_decaychain->locate(&cand);
279
280 if (pb && pb->tauIndex() >= 0 && pb->mother()) {
281 const int momindex = pb->momIndex();
282 const int tauIndex = pb->tauIndex();
283 const Eigen::Matrix<double, 1, 3> mom_vec = m_fitparams->getStateVector().segment(momindex, 3);
284
285 const Eigen::Matrix<double, 3, 3> mom_cov = m_fitparams->getCovariance().block<3, 3>(momindex, momindex);
286 Eigen::Matrix<double, 4, 4> comb_cov = Eigen::Matrix<double, 4, 4>::Zero(4, 4);
287
288 const std::tuple<double, double> lenTuple = getDecayLength(cand);
289
290 const double lenErr = std::get<1>(lenTuple);
291 comb_cov(0, 0) = lenErr * lenErr;
292 comb_cov(1, 0) = m_fitparams->getCovariance()(momindex, tauIndex);
293 comb_cov(2, 0) = m_fitparams->getCovariance()(momindex + 1, tauIndex);
294 comb_cov(3, 0) = m_fitparams->getCovariance()(momindex + 2, tauIndex);
295
296 comb_cov.block<3, 3>(1, 1) = mom_cov;
297
298 const double mass = pb->particle()->getPDGMass();
299 const double mBYc = mass / Belle2::Const::speedOfLight;
300 const double mom = mom_vec.norm();
301 const double mom3 = mom * mom * mom;
302
303 const double len = std::get<0>(lenTuple);
304 const double t = len / mom * mBYc;
305
306 Eigen::Matrix<double, 1, 4> jac = Eigen::Matrix<double, 1, 4>::Zero();
307 jac(0) = 1. / mom * mBYc;
308 jac(1) = -1. * len * mom_vec(0) / mom3 * mBYc;
309 jac(2) = -1. * len * mom_vec(1) / mom3 * mBYc;
310 jac(3) = -1. * len * mom_vec(2) / mom3 * mBYc;
311
312 const double tErr2 = jac * comb_cov.selfadjointView<Eigen::Lower>() * jac.transpose();
313 // time in nanosec
314 return std::make_tuple(t, std::sqrt(tErr2));
315 }
316 return std::make_tuple(-999, -999);
317 }
318
319 std::tuple<double, double> FitManager::getDecayLength(const ParticleBase* pb) const
320 {
321 // returns the decaylength in the lab frame
322 return getDecayLength(pb, *m_fitparams);
323 }
324
325 std::tuple<double, double> FitManager::getDecayLength(const ParticleBase* pb, const FitParams& fitparams) const
326 {
327 if (pb->tauIndex() >= 0 && pb->mother()) {
328 const int tauindex = pb->tauIndex();
329 const double len = fitparams.getStateVector()(tauindex);
330 const double lenErr2 = fitparams.getCovariance()(tauindex, tauindex);
331 return std::make_tuple(len, std::sqrt(lenErr2));
332 }
333 return std::make_tuple(-999, -999);
334 }
335
336 std::tuple<double, double> FitManager::getDecayLength(Belle2::Particle& cand) const
337 {
338 std::tuple<double, double> rc = std::make_tuple(-999, -999);
339 const ParticleBase* pb = m_decaychain->locate(&cand) ;
340 if (pb && pb->tauIndex() >= 0 && pb->mother()) {
341 rc = getDecayLength(pb);
342 }
343 return rc;
344 }
345
346}//end module namespace
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
Class to store reconstructed particles.
Definition: Particle.h:75
std::string getName() const override
Return name of this particle.
Definition: Particle.cc:1165
void writeExtraInfo(const std::string &name, const double value)
Sets the user defined extraInfo.
Definition: Particle.cc:1308
void setVertex(const ROOT::Math::XYZVector &vertex)
Sets position (decay vertex)
Definition: Particle.h:295
void set4VectorDividingByMomentumScaling(const ROOT::Math::PxPyPzEVector &p4)
Sets Lorentz vector dividing by the momentum scaling factor.
Definition: Particle.h:283
bool hasExtraInfo(const std::string &name) const
Return whether the extra info with the given name is set.
Definition: Particle.cc:1266
unsigned getNDaughters(void) const
Returns number of daughter particles.
Definition: Particle.h:727
double getPDGMass(void) const
Returns uncertainty on the invariant mass (requires valid momentum error matrix)
Definition: Particle.cc:604
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:393
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:366
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition: Particle.cc:631
double getExtraInfo(const std::string &name) const
Return given value if set.
Definition: Particle.cc:1289
constraint configuration class
const std::vector< int > m_massConstraintListPDG
list of pdg codes to mass constrain
this class does a lot of stuff: Build decaytree structure allowing to index particles and handle the ...
Definition: DecayChain.h:23
ErrCode initialize(FitParams &par)
initialize the chain
Definition: DecayChain.cc:70
const ParticleBase * locate(Belle2::Particle *bc) const
convert Belle2::particle into particle base(fitter base particle)
Definition: DecayChain.cc:111
int dim() const
get dimension
Definition: DecayChain.h:49
const ParticleBase * cand()
get candidate
Definition: DecayChain.h:64
ErrCode filterWithReference(FitParams &par, const FitParams &ref)
filter with respect to a previous iteration for better stability
Definition: DecayChain.cc:93
ErrCode filter(FitParams &par)
filter down the chain
Definition: DecayChain.cc:81
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
void reset()
reset the errorcode to default (success!)
Definition: ErrCode.h:55
const ConstraintConfiguration m_config
config container
Definition: FitManager.h:113
DecayChain * m_decaychain
the decay tree
Definition: FitManager.h:86
FitManager()
constructor
Definition: FitManager.h:29
Belle2::Particle * particle()
getter for the head of the tree
Definition: FitManager.h:79
ErrCode m_errCode
errorcode
Definition: FitManager.h:98
void getCovFromPB(const ParticleBase *pb, TMatrixFSym &returncov) const
extract cov from particle base
Definition: FitManager.cc:127
void updateTree(Belle2::Particle &particle, const bool isTreeHead) const
update the Belle2::Particles with the fit results
Definition: FitManager.cc:262
FitParams * m_fitparams
parameters to be fitted
Definition: FitManager.h:107
const bool m_updateDaugthers
if this is set all daughters will be updated otherwise only the head of the tree
Definition: FitManager.h:101
~FitManager()
destructor does stuff
Definition: FitManager.cc:45
int m_ndf
number of degrees of freedom for this topology
Definition: FitManager.h:104
double m_prec
precision that is needed for status:converged (delta chi2)
Definition: FitManager.h:95
bool fit()
main fit function that uses the kalman filter
Definition: FitManager.cc:51
int m_status
status of the current iteration
Definition: FitManager.h:89
Belle2::Particle * m_particle
head of the tree
Definition: FitManager.h:83
double m_chiSquare
chi2 of the current iteration
Definition: FitManager.h:92
bool m_useReferencing
use referencing
Definition: FitManager.h:110
VertexStatus
status flag of the fit-itereation (the step in the newton method)
Definition: FitManager.h:26
std::tuple< double, double > getLifeTime(Belle2::Particle &cand) const
get lifetime
Definition: FitManager.cc:276
std::tuple< double, double > getDecayLength(const ParticleBase *pb) const
get decay length
Definition: FitManager.cc:319
bool updateCand(Belle2::Particle &particle, const bool isTreeHead) const
update particles parameters with the fit results
Definition: FitManager.cc:191
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
bool testCovariance() const
test if the covariance makes sense
Definition: FitParams.cc:46
double chiSquare() const
get chi2 of statevector
Definition: FitParams.h:86
int nDof() const
get number of degrees of freedom
Definition: FitParams.cc:56
Eigen::Matrix< double, -1, -1, 0, MAX_MATRIX_SIZE, MAX_MATRIX_SIZE > & getCovariance()
getter for the states covariance
Definition: FitParams.h:53
base class for all particles
Definition: ParticleBase.h:25
Belle2::Particle * particle() const
get basf2 particle
Definition: ParticleBase.h:92
virtual int dim() const =0
get dimension of constraint
virtual void forceP4Sum(FitParams &) const
force p4 sum conservation all along the tree
Definition: ParticleBase.h:116
virtual int posIndex() const
get vertex index (in statevector!)
Definition: ParticleBase.h:122
virtual int momIndex() const
get momentum index
Definition: ParticleBase.h:128
virtual bool hasEnergy() const
get momentum dimension
Definition: ParticleBase.h:132
virtual int tauIndex() const
get tau index
Definition: ParticleBase.h:125
const ParticleBase * mother() const
getMother() / hasMother()
Definition: ParticleBase.h:98