Belle II Software  light-2205-abys
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 
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  {
43  config,
44  false
45  );
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 
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_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 
148  {
149  return m_decaychain->posIndex(particle);
150  }
151 
153  {
154  return m_decaychain->momIndex(particle);
155  }
156 
158  {
159  return m_decaychain->tauIndex(particle);
160  }
161 
162  void FitManager::getCovFromPB(const ParticleBase* pb, TMatrixFSym& returncov) const
163  {
164 
165  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> cov = m_fitparams->getCovariance().selfadjointView<Eigen::Lower>();
166  int posindex = pb->posIndex();
167  // hack: for tracks and photons, use the production vertex
168  if (posindex < 0 && pb->mother()) {
169  posindex = pb->mother()->posIndex();
170  }
171  int momindex = pb->momIndex();
172  if (pb->hasEnergy()) {
173  // if particle has energy, get full p4 from fitparams and put them directly in the return type
174  // very important! Belle2 uses p,E,x! Change order here!
175  for (int row = 0; row < 4; ++row) {
176  for (int col = 0; col < 4; ++col) {
177  returncov(row, col) = cov(momindex + row, momindex + col);
178  }
179  }
180 
181  for (int row = 0; row < 3; ++row) {
182  for (int col = 0; col < 3; ++col) {
183  returncov(row + 4, col + 4) = cov(posindex + row, posindex + col);
184  }
185  }
186 
187  } else {
188  Eigen::Matrix<double, 6, 6> cov6 =
189  Eigen::Matrix<double, 6, 6>::Zero(6, 6);
190 
191  for (int row = 0; row < 3; ++row) {
192  for (int col = 0; col < 3; ++col) {
193  cov6(row, col) = cov(momindex + row, momindex + col);
194  cov6(row + 3, col + 3) = cov(posindex + row, posindex + col);
195  }
196  }
197 
198  double mass = pb->pdgMass();
199  Eigen::Matrix<double, 3, 1> momVec =
200  m_fitparams->getStateVector().segment(momindex, 3);
201 
202  double energy2 = momVec.transpose() * momVec;
203  energy2 += mass * mass;
204  double energy = sqrt(energy2);
205 
206  Eigen::Matrix<double, 7, 6> jacobian =
207  Eigen::Matrix<double, 7, 6>::Zero(7, 6);
208 
209  for (int col = 0; col < 3; ++col) {
210  jacobian(col, col) = 1;
211  jacobian(3, col) = m_fitparams->getStateVector()(momindex + col) / energy;
212  jacobian(col + 4, col + 3) = 1;
213  }
214 
215  Eigen::Matrix<double, 7, 7> cov7
216  = jacobian * cov6.selfadjointView<Eigen::Lower>() * jacobian.transpose();
217 
218  for (int row = 0; row < 7; ++row) {
219  for (int col = 0; col < 7; ++col) {
220  returncov(row, col) = cov7(row, col);
221  }
222  }
223  } // else
224  }
225 
226  bool FitManager::updateCand(Belle2::Particle& cand, const bool isTreeHead) const
227  {
228  const ParticleBase* pb = m_decaychain->locate(&cand);
229  if (pb) {
230  updateCand(*pb, cand, isTreeHead);
231  } else {
232  B2ERROR("Can't find candidate " << cand.getName() << " in tree " << m_particle->getName());
233  }
234  return pb != nullptr;
235  }
236 
238  Belle2::Particle& cand, const bool isTreeHead) const
239  {
240  int posindex = pb.posIndex();
241  if (posindex < 0 && pb.mother()) {
242  posindex = pb.mother()->posIndex();
243  }
244 
245  if (m_updateDaugthers || isTreeHead) {
246  if (posindex >= 0) {
247  const ROOT::Math::XYZVector 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  if (pb.mother()) {
259  int motherPosIndex = pb.mother()->posIndex();
260  if (motherPosIndex >= 0) {
261  setExtraInfo(&cand, "prodVertexX", m_fitparams->getStateVector()(motherPosIndex));
262  setExtraInfo(&cand, "prodVertexY", m_fitparams->getStateVector()(motherPosIndex + 1));
263  if (pb.mother()->dim() > 2)
264  setExtraInfo(&cand, "prodVertexZ", m_fitparams->getStateVector()(motherPosIndex + 2));
265  }
266  }
267  }
268 
269  const int momindex = pb.momIndex();
270  ROOT::Math::PxPyPzEVector p;
271  p.SetPx(m_fitparams->getStateVector()(momindex));
272  p.SetPy(m_fitparams->getStateVector()(momindex + 1));
273  p.SetPz(m_fitparams->getStateVector()(momindex + 2));
274  if (pb.hasEnergy()) {
275  p.SetE(m_fitparams->getStateVector()(momindex + 3));
276  cand.set4Vector(p);
277  } else {
278  const double mass = cand.getPDGMass();
279  p.SetE(std::sqrt(p.P2() + mass * mass));
280  cand.set4Vector(p);
281  }
282  TMatrixFSym cov7b2(7);
283  getCovFromPB(&pb, cov7b2);
284  cand.setMomentumVertexErrorMatrix(cov7b2);
285  }
286 
287  if (pb.tauIndex() > 0) {
288  std::tuple<double, double>tau = getDecayLength(cand);
289  std::tuple<double, double>life = getLifeTime(cand);
290  setExtraInfo(&cand, std::string("decayLength"), std::get<0>(tau));
291  setExtraInfo(&cand, std::string("decayLengthErr"), std::get<1>(tau));
292  setExtraInfo(&cand, std::string("lifeTime"), std::get<0>(life));
293  setExtraInfo(&cand, std::string("lifeTimeErr"), std::get<1>(life));
294  }
295  }
296 
297  void FitManager::updateTree(Belle2::Particle& cand, const bool isTreeHead) const
298  {
299  const bool updateableMother = updateCand(cand, isTreeHead);
300 
301  if (updateableMother and not cand.hasExtraInfo("bremsCorrected")) {
302  const int ndaughters = cand.getNDaughters();
303  for (int i = 0; i < ndaughters; i++) {
304  auto* daughter = const_cast<Belle2::Particle*>(cand.getDaughter(i));
305  updateTree(*daughter, false);
306  }
307  }
308  }
309 
310  std::tuple<double, double> FitManager::getLifeTime(Belle2::Particle& cand) const
311  {
312  const ParticleBase* pb = m_decaychain->locate(&cand);
313 
314  if (pb && pb->tauIndex() >= 0 && pb->mother()) {
315  const int momindex = pb->momIndex();
316  const int tauIndex = pb->tauIndex();
317  const Eigen::Matrix<double, 1, 3> mom_vec = m_fitparams->getStateVector().segment(momindex, 3);
318 
319  const Eigen::Matrix<double, 3, 3> mom_cov = m_fitparams->getCovariance().block<3, 3>(momindex, momindex);
320  Eigen::Matrix<double, 4, 4> comb_cov = Eigen::Matrix<double, 4, 4>::Zero(4, 4);
321 
322  std::tuple<double, double> lenTuple = getDecayLength(cand);
323 
324  const double lenErr = std::get<1>(lenTuple);
325  comb_cov(0, 0) = lenErr * lenErr;
326  comb_cov(1, 0) = m_fitparams->getCovariance()(momindex, tauIndex);
327  comb_cov(2, 0) = m_fitparams->getCovariance()(momindex + 1, tauIndex);
328  comb_cov(3, 0) = m_fitparams->getCovariance()(momindex + 2, tauIndex);
329 
330  comb_cov.block<3, 3>(1, 1) = mom_cov;
331 
332  const double mass = pb->pdgMass();
333  const double mBYc = mass / Belle2::Const::speedOfLight;
334  const double mom = mom_vec.norm();
335  const double mom3 = mom * mom * mom;
336 
337  const double len = std::get<0>(lenTuple);
338  const double t = len / mom * mBYc;
339 
340  Eigen::Matrix<double, 1, 4> jac = Eigen::Matrix<double, 1, 4>::Zero();
341  jac(0) = 1. / mom * mBYc;
342  jac(1) = -1. * len * mom_vec(0) / mom3 * mBYc;
343  jac(2) = -1. * len * mom_vec(1) / mom3 * mBYc;
344  jac(3) = -1. * len * mom_vec(2) / mom3 * mBYc;
345 
346  const double tErr2 = jac * comb_cov.selfadjointView<Eigen::Lower>() * jac.transpose();
347  // time in nanosec
348  return std::make_tuple(t, std::sqrt(tErr2));
349  }
350  return std::make_tuple(-999, -999);
351  }
352 
353  std::tuple<double, double> FitManager::getDecayLength(const ParticleBase* pb) const
354  {
355  // returns the decaylength in the lab frame
356  return getDecayLength(pb, *m_fitparams);
357  }
358 
359  std::tuple<double, double> FitManager::getDecayLength(const ParticleBase* pb, const FitParams& fitparams) const
360  {
361  if (pb->tauIndex() >= 0 && pb->mother()) {
362  const int tauindex = pb->tauIndex();
363  const double len = fitparams.getStateVector()(tauindex);
364  const double lenErr2 = fitparams.getCovariance()(tauindex, tauindex);
365  return std::make_tuple(len, std::sqrt(lenErr2));
366  }
367  return std::make_tuple(-999, -999);
368  }
369 
370  std::tuple<double, double> FitManager::getDecayLength(Belle2::Particle& cand) const
371  {
372  std::tuple<double, double> rc = std::make_tuple(-999, -999);
373  const ParticleBase* pb = m_decaychain->locate(&cand) ;
374  if (pb && pb->tauIndex() >= 0 && pb->mother()) {
375  rc = getDecayLength(pb);
376  }
377  return rc;
378  }
379 
380 }//end module namespace
static const double speedOfLight
[cm/ns]
Definition: Const.h:575
Class to store reconstructed particles.
Definition: Particle.h:74
void setExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1306
std::string getName() const override
Return name of this particle.
Definition: Particle.cc:1157
void setVertex(const ROOT::Math::XYZVector &vertex)
Sets position (decay vertex)
Definition: Particle.h:288
bool hasExtraInfo(const std::string &name) const
Return whether the extra info with the given name is set.
Definition: Particle.cc:1255
unsigned getNDaughters(void) const
Returns number of daughter particles.
Definition: Particle.h:656
double getPDGMass(void) const
Returns uncertainty on the invariant mass (requires valid momentum error matrix)
Definition: Particle.cc:636
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1325
void set4Vector(const ROOT::Math::PxPyPzEVector &p4)
Sets Lorentz vector.
Definition: Particle.h:276
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:425
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:331
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition: Particle.cc:654
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:74
const ParticleBase * cand()
get candidate
Definition: DecayChain.h:64
const ParticleBase * locate(Belle2::Particle *bc) const
convert Belle2::particle into particle base(fitter base particle)
Definition: DecayChain.cc:115
int momIndex(Belle2::Particle *bc) const
get momentum index of the particle in the state vector
Definition: DecayChain.cc:152
int dim() const
get dimension
Definition: DecayChain.h:49
ErrCode filterWithReference(FitParams &par, const FitParams &ref)
filter with respect to a previous iteration for better stability
Definition: DecayChain.cc:97
int tauIndex(Belle2::Particle *bc) const
get tau (i.e.
Definition: DecayChain.cc:161
ErrCode filter(FitParams &par)
filter down the chain
Definition: DecayChain.cc:85
int posIndex(Belle2::Particle *bc) const
get the vertex index of the particle in state vector
Definition: DecayChain.cc:143
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
int momIndex(Belle2::Particle *particle) const
getter for the index of the momentum in the state vector
Definition: FitManager.cc:152
const ConstraintConfiguration m_config
config container
Definition: FitManager.h:146
DecayChain * m_decaychain
the decay tree
Definition: FitManager.h:116
FitManager()
constructor
Definition: FitManager.h:29
int posIndex(Belle2::Particle *particle) const
getter for the index of the vertex position in the state vector
Definition: FitManager.cc:147
int tauIndex(Belle2::Particle *particle) const
getter for the index of tau, the lifetime in the statevector
Definition: FitManager.cc:157
void setExtraInfo(Belle2::Particle *part, const std::string &name, const double value) const
add extrainfo to particle
Definition: FitManager.cc:55
ErrCode m_errCode
errorcode
Definition: FitManager.h:131
int nDof() const
getter for degrees of freedom of the fitparameters
Definition: FitManager.cc:142
void getCovFromPB(const ParticleBase *pb, TMatrixFSym &returncov) const
extract cov from particle base
Definition: FitManager.cc:162
void updateTree(Belle2::Particle &particle, const bool isTreeHead) const
update the Belle2::Particles with the fit results
Definition: FitManager.cc:297
FitParams * m_fitparams
parameters to be fitted
Definition: FitManager.h:140
const bool m_updateDaugthers
if this is set all daughters will be updated otherwise only the head of the tree
Definition: FitManager.h:134
~FitManager()
destructor does stuff
Definition: FitManager.cc:49
int m_ndf
number of degrees of freedom for this topology
Definition: FitManager.h:137
double m_prec
precision that is needed for status:converged (delta chi2)
Definition: FitManager.h:128
bool fit()
main fit function that uses the kalman filter
Definition: FitManager.cc:66
int m_status
status of the current iteration
Definition: FitManager.h:119
Belle2::Particle * particle()
getter for the head of the tree
Definition: FitManager.h:109
Belle2::Particle * m_particle
head of the tree
Definition: FitManager.h:113
int m_niter
iteration index
Definition: FitManager.h:125
double m_chiSquare
chi2 of the current iteration
Definition: FitManager.h:122
bool m_useReferencing
use referencing
Definition: FitManager.h:143
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:310
std::tuple< double, double > getDecayLength(const ParticleBase *pb) const
get decay length
Definition: FitManager.cc:353
bool updateCand(Belle2::Particle &particle, const bool isTreeHead) const
update particles parameters with the fit results
Definition: FitManager.cc:226
Class to store and manage fitparams (statevector)
Definition: FitParams.h:20
Eigen::Matrix< double, -1, -1, 0, MAX_MATRIX_SIZE, MAX_MATRIX_SIZE > & getCovariance()
getter for the states covariance
Definition: FitParams.h:53
bool testCovariance() const
test if the covariance makes sense
Definition: FitParams.cc:46
Eigen::Matrix< double, -1, 1, 0, MAX_MATRIX_SIZE, 1 > & getStateVector()
getter for the fit parameters/statevector
Definition: FitParams.h:65
double chiSquare() const
get chi2 of statevector
Definition: FitParams.cc:56
int nDof() const
get numer of degrees of freedom
Definition: FitParams.cc:61
base class for all particles
Definition: ParticleBase.h:25
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:123
virtual int posIndex() const
get vertex index (in statevector!)
Definition: ParticleBase.h:129
virtual int momIndex() const
get momentum index
Definition: ParticleBase.h:135
virtual bool hasEnergy() const
get momentum dimension
Definition: ParticleBase.h:139
virtual int tauIndex() const
get tau index
Definition: ParticleBase.h:132
const ParticleBase * mother() const
getMother() / hasMother()
double pdgMass() const
get pdg mass
Definition: ParticleBase.h:151