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