Belle II Software light-2406-ragdoll
FitManager Class Reference

this class More...

#include <FitManager.h>

Collaboration diagram for FitManager:

Public Types

enum  VertexStatus {
  Success = 0 ,
  NonConverged ,
  BadInput ,
  Failed ,
  UnFitted
}
 status flag of the fit-itereation (the step in the newton method) More...
 

Public Member Functions

 FitManager ()
 constructor

 
 FitManager (Belle2::Particle *particle, const ConstraintConfiguration &config, double prec=0.01, bool updateDaughters=false, const bool useReferencing=false)
 constructor

 
 FitManager (const FitManager &other)=delete
 use default copy constructor
 
FitManageroperator= (const FitManager &other)=delete
 use default assignment op
 
 ~FitManager ()
 destructor does stuff
 
bool fit ()
 main fit function that uses the kalman filter
 
bool updateCand (Belle2::Particle &particle, const bool isTreeHead) const
 update particles parameters with the fit results
 
void updateCand (const ParticleBase &pb, Belle2::Particle &cand, const bool isTreeHead) const
 locate particle base for a belle2 particle and update the particle with the values from particle base
 
void updateTree (Belle2::Particle &particle, const bool isTreeHead) const
 update the Belle2::Particles with the fit results

 
void getCovFromPB (const ParticleBase *pb, TMatrixFSym &returncov) const
 extract cov from particle base
 
std::tuple< double, double > getLifeTime (Belle2::Particle &cand) const
 get lifetime
 
std::tuple< double, double > getDecayLength (const ParticleBase *pb) const
 get decay length
 
std::tuple< double, double > getDecayLength (const ParticleBase *pb, const FitParams &fitparams) const
 get decay length
 
std::tuple< double, double > getDecayLength (Belle2::Particle &cand) const
 get decay length
 
Belle2::Particleparticle ()
 getter for the head of the tree
 

Private Attributes

Belle2::Particlem_particle
 head of the tree

 
DecayChainm_decaychain
 the decay tree
 
int m_status
 status of the current iteration
 
double m_chiSquare
 chi2 of the current iteration
 
double m_prec
 precision that is needed for status:converged (delta chi2)
 
ErrCode m_errCode
 errorcode
 
const bool m_updateDaugthers
 if this is set all daughters will be updated otherwise only the head of the tree
 
int m_ndf
 number of degrees of freedom for this topology
 
FitParamsm_fitparams
 parameters to be fitted
 
bool m_useReferencing
 use referencing
 
const ConstraintConfiguration m_config
 config container
 

Detailed Description

this class

Definition at line 22 of file FitManager.h.

Member Enumeration Documentation

◆ VertexStatus

status flag of the fit-itereation (the step in the newton method)

Definition at line 26 of file FitManager.h.

26{ Success = 0, NonConverged, BadInput, Failed, UnFitted };

Constructor & Destructor Documentation

◆ FitManager() [1/2]

FitManager ( )
inline

constructor

Definition at line 29 of file FitManager.h.

29 : m_particle(0), m_decaychain(0), m_status(VertexStatus::UnFitted),
30 m_chiSquare(-1), m_prec(0.01), m_updateDaugthers(false), m_ndf(0),
32 {}
DecayChain * m_decaychain
the decay tree
Definition: FitManager.h:86
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
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
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

◆ FitManager() [2/2]

FitManager ( Belle2::Particle particle,
const ConstraintConfiguration config,
double  prec = 0.01,
bool  updateDaughters = false,
const bool  useReferencing = false 
)

constructor

Definition at line 24 of file FitManager.cc.

29 :
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);
42 m_fitparams = new FitParams(m_decaychain->dim());
43 }
int dim() const
get dimension
Definition: DecayChain.h:49
const ConstraintConfiguration m_config
config container
Definition: FitManager.h:113
Belle2::Particle * particle()
getter for the head of the tree
Definition: FitManager.h:79

◆ ~FitManager()

~FitManager ( )

destructor does stuff

Definition at line 45 of file FitManager.cc.

46 {
47 delete m_decaychain;
48 delete m_fitparams;
49 }

Member Function Documentation

◆ fit()

bool fit ( )

main fit function that uses the kalman filter

Definition at line 51 of file FitManager.cc.

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 }
void writeExtraInfo(const std::string &name, const double value)
Sets the user defined extraInfo.
Definition: Particle.cc:1308
const std::vector< int > m_massConstraintListPDG
list of pdg codes to mass constrain
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
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
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
ErrCode m_errCode
errorcode
Definition: FitManager.h:98
void updateTree(Belle2::Particle &particle, const bool isTreeHead) const
update the Belle2::Particles with the fit results
Definition: FitManager.cc:262
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
virtual void forceP4Sum(FitParams &) const
force p4 sum conservation all along the tree
Definition: ParticleBase.h:116

◆ getCovFromPB()

void getCovFromPB ( const ParticleBase pb,
TMatrixFSym &  returncov 
) const

extract cov from particle base

Definition at line 127 of file FitManager.cc.

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 }
Eigen::Matrix< double, -1, 1, 0, MAX_MATRIX_SIZE, 1 > & getStateVector()
getter for the fit parameters/statevector
Definition: FitParams.h:65
Eigen::Matrix< double, -1, -1, 0, MAX_MATRIX_SIZE, MAX_MATRIX_SIZE > & getCovariance()
getter for the states covariance
Definition: FitParams.h:53

◆ getDecayLength() [1/3]

std::tuple< double, double > getDecayLength ( Belle2::Particle cand) const

get decay length

Definition at line 336 of file FitManager.cc.

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 }
std::tuple< double, double > getDecayLength(const ParticleBase *pb) const
get decay length
Definition: FitManager.cc:319

◆ getDecayLength() [2/3]

std::tuple< double, double > getDecayLength ( const ParticleBase pb) const

get decay length

Definition at line 319 of file FitManager.cc.

320 {
321 // returns the decaylength in the lab frame
322 return getDecayLength(pb, *m_fitparams);
323 }

◆ getDecayLength() [3/3]

std::tuple< double, double > getDecayLength ( const ParticleBase pb,
const FitParams fitparams 
) const

get decay length

Definition at line 325 of file FitManager.cc.

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 }

◆ getLifeTime()

std::tuple< double, double > getLifeTime ( Belle2::Particle cand) const

get lifetime

Definition at line 276 of file FitManager.cc.

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 }
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
virtual int momIndex() const
get momentum index
Definition: ParticleBase.h:128

◆ particle()

Belle2::Particle * particle ( )
inline

getter for the head of the tree

Definition at line 79 of file FitManager.h.

79{ return m_particle; }

◆ updateCand() [1/2]

bool updateCand ( Belle2::Particle particle,
const bool  isTreeHead 
) const

update particles parameters with the fit results

Definition at line 191 of file FitManager.cc.

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 }
std::string getName() const override
Return name of this particle.
Definition: Particle.cc:1165
bool updateCand(Belle2::Particle &particle, const bool isTreeHead) const
update particles parameters with the fit results
Definition: FitManager.cc:191

◆ updateCand() [2/2]

void updateCand ( const ParticleBase pb,
Belle2::Particle cand,
const bool  isTreeHead 
) const

locate particle base for a belle2 particle and update the particle with the values from particle base

Definition at line 202 of file FitManager.cc.

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 }
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
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 ParticleBase * cand()
get candidate
Definition: DecayChain.h:64
void getCovFromPB(const ParticleBase *pb, TMatrixFSym &returncov) const
extract cov from particle base
Definition: FitManager.cc:127
std::tuple< double, double > getLifeTime(Belle2::Particle &cand) const
get lifetime
Definition: FitManager.cc:276

◆ updateTree()

void updateTree ( Belle2::Particle particle,
const bool  isTreeHead 
) const

update the Belle2::Particles with the fit results

Definition at line 262 of file FitManager.cc.

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 }
Class to store reconstructed particles.
Definition: Particle.h:75

Member Data Documentation

◆ m_chiSquare

double m_chiSquare
private

chi2 of the current iteration

Definition at line 92 of file FitManager.h.

◆ m_config

const ConstraintConfiguration m_config
private

config container

Definition at line 113 of file FitManager.h.

◆ m_decaychain

DecayChain* m_decaychain
private

the decay tree

Definition at line 86 of file FitManager.h.

◆ m_errCode

ErrCode m_errCode
private

errorcode

Definition at line 98 of file FitManager.h.

◆ m_fitparams

FitParams* m_fitparams
private

parameters to be fitted

Definition at line 107 of file FitManager.h.

◆ m_ndf

int m_ndf
private

number of degrees of freedom for this topology

Definition at line 104 of file FitManager.h.

◆ m_particle

Belle2::Particle* m_particle
private

head of the tree

Definition at line 83 of file FitManager.h.

◆ m_prec

double m_prec
private

precision that is needed for status:converged (delta chi2)

Definition at line 95 of file FitManager.h.

◆ m_status

int m_status
private

status of the current iteration

Definition at line 89 of file FitManager.h.

◆ m_updateDaugthers

const bool m_updateDaugthers
private

if this is set all daughters will be updated otherwise only the head of the tree

Definition at line 101 of file FitManager.h.

◆ m_useReferencing

bool m_useReferencing
private

use referencing

Definition at line 110 of file FitManager.h.


The documentation for this class was generated from the following files: