Belle II Software  light-2205-abys
InternalParticle.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/dataobjects/Particle.h>
11 
12 #include <analysis/VertexFitting/TreeFitter/InternalParticle.h>
13 #include <analysis/VertexFitting/TreeFitter/FitParams.h>
14 #include <analysis/VertexFitting/TreeFitter/HelixUtils.h>
15 #include <framework/logging/Logger.h>
16 #include <mdst/dataobjects/Track.h>
17 
18 using std::vector;
19 
20 namespace TreeFitter {
21 
22  inline bool sortByType(const ParticleBase* lhs, const ParticleBase* rhs)
23  {
24  int lhstype = lhs->type() ;
25  int rhstype = rhs->type() ;
26  bool rc = false ;
27  if (lhstype == rhstype &&
28  lhstype == ParticleBase::TFParticleType::kRecoTrack) {
29 
30  rc = lhs->particle()->getMomentum().Rho() > rhs->particle()->getMomentum().Rho();
31  } else if (lhs->particle() && rhs->particle() && lhs->particle()->getNDaughters() > 0 &&
32  rhs->particle()->getNDaughters() > 0) {
33  rc = lhs->nFinalChargedCandidates() > rhs->nFinalChargedCandidates();
34  } else {
35  rc = lhstype < rhstype;
36  }
37  return rc;
38  }
39 
41  const ParticleBase* mother,
42  const ConstraintConfiguration& config,
43  bool forceFitAll
44  ) :
45  ParticleBase(particle, mother, &config),// config pointer here to allow final states not to have it
46  m_massconstraint(false),
47  m_lifetimeconstraint(false),
48  m_isconversion(false),
49  m_automatic_vertex_constraining(config.m_automatic_vertex_constraining)
50  {
51  if (particle) {
52  for (auto daughter : particle->getDaughters()) {
53  addDaughter(daughter, config, forceFitAll);
54  }
55  } else {
56  B2ERROR("Trying to create an InternalParticle from NULL. This should never happen.");
57  }
58 
59  m_massconstraint = std::find(config.m_massConstraintListPDG.begin(), config.m_massConstraintListPDG.end(),
60  std::abs(m_particle->getPDGCode())) != config.m_massConstraintListPDG.end();
61 
63  // if this is a hadronically decaying resonance it is useful to constraint the decay vertex to its mothers decay vertex.
64  //
66  config.m_fixedToMotherVertexListPDG.end(),
67  std::abs(m_particle->getPDGCode())) != config.m_fixedToMotherVertexListPDG.end() && this->mother();
68 
69  // use geo constraint if this particle is in the list to constrain
70  m_geo_constraint = std::find(config.m_geoConstraintListPDG.begin(),
71  config.m_geoConstraintListPDG.end(),
72  std::abs(m_particle->getPDGCode())) != config.m_geoConstraintListPDG.end() && this->mother() && !m_shares_vertex_with_mother;
73  } else {
76  }
77  }
78 
80  {
81 
82  return lhs->particle()->getMomentum().Rho() > rhs->particle()->getMomentum().Rho();
83  }
84 
86  {
87  ErrCode status ;
88  for (auto daughter : m_daughters) {
89  status |= daughter->initMotherlessParticle(fitparams);
90  }
91 
92  int posindex = posIndex();
93  // FIXME update this and check if this position is already initialised
94  // lower stream vertices might be better
95  if (hasPosition()) {
96  fitparams.getStateVector().segment(posindex, 3) = Eigen::Matrix<double, 3, 1>::Zero(3);
97 
98  std::vector<ParticleBase*> alldaughters;
99  ParticleBase::collectVertexDaughters(alldaughters, posindex);
100 
101  std::vector<ParticleBase*> vtxdaughters;
102 
103  vector<RecoTrack*> trkdaughters;
104  for (auto daughter : alldaughters) {
105  if (daughter->type() == ParticleBase::TFParticleType::kRecoTrack) {
106  trkdaughters.push_back(static_cast<RecoTrack*>(daughter));
107  } else if (daughter->hasPosition()
108  && fitparams.getStateVector()(daughter->posIndex()) != 0) {
109  vtxdaughters.push_back(daughter);
110  }
111  }
112 
114 
115  if (trkdaughters.size() >= 2) {
116  std::sort(trkdaughters.begin(), trkdaughters.end(), compTrkTransverseMomentum);
117 
118  RecoTrack* dau1 = trkdaughters[0];
119  RecoTrack* dau2 = trkdaughters[1];
120 
122  dau1->particle()->getPDGCode())))->getHelix();
124  dau2->particle()->getPDGCode())))->getHelix();
125 
126  double flt1(0), flt2(0);
127  HelixUtils::helixPoca(helix1, helix2, flt1, flt2, v, m_isconversion);
128 
129  fitparams.getStateVector()(posindex) = v.x();
130  fitparams.getStateVector()(posindex + 1) = v.y();
131  fitparams.getStateVector()(posindex + 2) = v.z();
132 
133  dau1->setFlightLength(flt1);
134  dau2->setFlightLength(flt2);
135 
136  } else if (false && trkdaughters.size() + vtxdaughters.size() >= 2) {
137  // TODO switched off waiting for refactoring of init1 and init2 functions (does not affect performance)
138  } else if (mother() && mother()->posIndex() >= 0) {
139  const int posindexmother = mother()->posIndex();
140  const int dim = m_config->m_originDimension; //TODO access mother
141  fitparams.getStateVector().segment(posindex, dim) = fitparams.getStateVector().segment(posindexmother, dim);
142  } else {
144  fitparams.getStateVector().segment(posindex, 3) = Eigen::Matrix<double, 1, 3>::Zero(3);
145  }
146  }
147 
148  for (auto daughter : m_daughters) {
149  daughter->initParticleWithMother(fitparams);
150  }
151 
152  initMomentum(fitparams);
153  return status;
154  }
155 
157  {
158  int posindex = posIndex();
159  if (hasPosition() &&
160  mother() &&
161  fitparams.getStateVector()(posindex) == 0 &&
162  fitparams.getStateVector()(posindex + 1) == 0 && \
163  fitparams.getStateVector()(posindex + 2) == 0) {
164  const int posindexmom = mother()->posIndex();
165  const int dim = m_config->m_originDimension; //TODO access mother?
166  fitparams.getStateVector().segment(posindex, dim) = fitparams.getStateVector().segment(posindexmom, dim);
167  }
168  return initTau(fitparams);
169  }
170 
172  {
173  int momindex = momIndex();
174  fitparams.getStateVector().segment(momindex, 4) = Eigen::Matrix<double, 4, 1>::Zero(4);
175 
176  for (auto daughter : m_daughters) {
177  int daumomindex = daughter->momIndex();
178  int maxrow = daughter->hasEnergy() ? 4 : 3;
179 
180  double e2 = fitparams.getStateVector().segment(daumomindex, maxrow).squaredNorm();
181  fitparams.getStateVector().segment(momindex, maxrow) += fitparams.getStateVector().segment(daumomindex, maxrow);
182 
183  if (maxrow == 3) {
184  double mass = daughter->pdgMass();
185  fitparams.getStateVector()(momindex + 3) += std::sqrt(e2 + mass * mass);
186  }
187  }
188  return ErrCode(ErrCode::Status::success);
189  }
190 
192  {
193  ErrCode status;
194  ParticleBase::initCovariance(fitparams);
195  for (auto daughter : m_daughters) {
196  status |= daughter->initCovariance(fitparams);
197  }
198  return status;
199  }
200 
201 
203  Projection& p) const
204  {
205  const int momindex = momIndex();
206 
207  // `this` always has an energy row
208  p.getResiduals().segment(0, 4) = fitparams.getStateVector().segment(momindex, 4);
209 
210  for (int imom = 0; imom < 4; ++imom) {
211  p.getH()(imom, momindex + imom) = 1;
212  }
213 
214  for (const auto daughter : m_daughters) {
215  const int daumomindex = daughter->momIndex();
216  const Eigen::Matrix<double, 1, 3> p3_vec = fitparams.getStateVector().segment(daumomindex, 3);
217 
218  // three momentum is easy just subtract the vectors
219  p.getResiduals().segment(0, 3) -= p3_vec;
220 
221  // energy depends on the parametrisation!
222  if (daughter->hasEnergy()) {
223  p.getResiduals()(3) -= fitparams.getStateVector()(daumomindex + 3);
224  p.getH()(3, daumomindex + 3) = -1; // d/dE -E
225  } else {
226  // m^2 + p^2 = E^2
227  // so
228  // E = sqrt(m^2 + p^2)
229  const double mass = daughter->pdgMass();
230  const double p2 = p3_vec.squaredNorm();
231  const double energy = std::sqrt(mass * mass + p2);
232  p.getResiduals()(3) -= energy;
233 
234  for (unsigned i = 0; i < 3; ++i) {
235  // d/dpx_i sqrt(m^2 + p^2)
236  p.getH()(3, daumomindex + i) = -1 * p3_vec(i) / energy;
237  }
238  }
239 
240  // this has to be in any case
241  // d/dp_i p_i
242  for (unsigned i = 0; i < 3; ++i) {
243  p.getH()(i, daumomindex + i) = -1;
244  }
245  }
246  return ErrCode(ErrCode::Status::success);
247  }
248 
249 
251  const FitParams& fitparams,
252  Projection& p) const
253  {
254  ErrCode status;
255  switch (type) {
256  case Constraint::mass:
257  status |= projectMassConstraint(fitparams, p);
258  break;
259  case Constraint::geometric:
260  status |= projectGeoConstraint(fitparams, p);
261  break;
262  case Constraint::kinematic:
263  status |= projectKineConstraint(fitparams, p);
264  break;
265  default:
266  status |= ParticleBase::projectConstraint(type, fitparams, p);
267  }
268 
269  return status;
270  }
271 
273  {
274  // for example B0 and D* can share the same vertex
275  return m_shares_vertex_with_mother ? this->mother()->posIndex() : index();
276  }
277 
279  {
280  // { x, y, z, tau, px, py, pz, E }
281  // the last 4 always exists for composite particles
282  // tau index only exist with a vertex and geo constraint
283  //
284  if (m_shares_vertex_with_mother) { return 4; }
285  else if (!m_geo_constraint) { return 7; }
286  else { return 8; }
287  }
288 
290  {
292  return m_geo_constraint ? index() + 3 : -1;
293  }
294 
296  {
301  if (m_geo_constraint && !m_shares_vertex_with_mother) { return this->index() + 4; }
302 
303  if (m_shares_vertex_with_mother) { return this->index(); }
304 
305  if (!m_geo_constraint) {return index() + 3 ;}
306 
307  // this will crash the initialisation
308  return -1;
309  }
310 
312  {
313  // does this particle have a position index (decayvertex)
314  // true in any case
315  return true;
316  }
317 
319  int depth) const
320  {
321 
322  for (auto daughter : m_daughters) {
323  daughter->addToConstraintList(list, depth - 1);
324  }
325  if (tauIndex() >= 0 && m_lifetimeconstraint) {
326  list.push_back(Constraint(this, Constraint::lifetime, depth, 1));
327  }
328  if (momIndex() >= 0) {
329  list.push_back(Constraint(this, Constraint::kinematic, depth, 4, 3));
330  }
331  if (m_geo_constraint) {
332  assert(m_config);
333  const int dim = m_config->m_originDimension == 2 && std::abs(m_particle->getPDGCode()) == m_config->m_headOfTreePDG ? 2 : 3;
334  list.push_back(Constraint(this, Constraint::geometric, depth, dim, 3));
335  }
336  if (m_massconstraint) {
337  list.push_back(Constraint(this, Constraint::mass, depth, 1, 3));
338  }
339 
340  }
341 
342  std::string InternalParticle::parname(int thisindex) const
343  {
344  int id = thisindex;
345  if (!mother() && id >= 3) {++id;}
346  return ParticleBase::parname(id);
347  }
348 
350  {
351  for (const auto daughter : m_daughters) {
352  daughter->forceP4Sum(fitparams);
353  }
354  const int momindex = momIndex();
355  if (momindex > 0) {
356  const int dim = hasEnergy() ? 4 : 3;
357  Projection p(fitparams.getDimensionOfState(), dim);
358  projectKineConstraint(fitparams, p);
359  fitparams.getStateVector().segment(momindex, dim) -= p.getResiduals().segment(0, dim);
360  }
361  }
362 
363 }
364 
DataType y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:421
DataType z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:423
DataType x() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:419
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:470
This class represents an ideal helix in perigee parameterization.
Definition: Helix.h:59
Class to store reconstructed particles.
Definition: Particle.h:74
const Track * getTrack() const
Returns the pointer to the Track object that was used to create this Particle (ParticleType == c_Trac...
Definition: Particle.cc:837
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:403
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:660
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:497
Helix getHelix() const
Conversion to framework Helix (without covariance).
const TrackFitResult * getTrackFitResultWithClosestMass(const Const::ChargedStable &requestedType) const
Return the track fit for a fit hypothesis with the closest mass.
Definition: Track.cc:68
const int m_originDimension
dimension of the origin constraint and ALL geometric gcosntraints
const std::vector< int > m_geoConstraintListPDG
list of pdg codes to mass constrain
const std::vector< int > m_fixedToMotherVertexListPDG
list of pdg codes to mass constrain
int m_headOfTreePDG
PDG code of the head particle.
const std::vector< int > m_massConstraintListPDG
list of pdg codes to mass constrain
class to manage the order of constraints and their filtering
Definition: Constraint.h:20
Type
type of constraints the order of these constraints is important: it is the order in which they are ap...
Definition: Constraint.h:27
abstract errorocode be aware that the default is success
Definition: ErrCode.h:14
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
int getDimensionOfState() const
get the states dimension
Definition: FitParams.h:92
static double helixPoca(const Belle2::Helix &helix1, const Belle2::Helix &helix2, double &flt1, double &flt2, Belle2::B2Vector3D &vertex, bool parallel=false)
POCA between two tracks.
Definition: HelixUtils.cc:224
virtual void forceP4Sum(FitParams &) const override
enforce conservation of momentum sum
virtual int tauIndex() const override
tau index in fit params only if it has a mother
bool m_isconversion
is conversion
virtual std::string parname(int index) const override
name
ErrCode projectKineConstraint(const FitParams &, Projection &) const
project kinematical constraint
virtual int dim() const override
space reserved in fit params, if has mother then it has tau
bool m_massconstraint
has mass constraint
virtual ErrCode initCovariance(FitParams &) const override
init covariance
virtual int type() const override
type
bool m_automatic_vertex_constraining
automatically figure out if mother and particle vertex should be the same and also add geometric cons...
virtual bool hasEnergy() const override
has energy in fitparams
virtual bool hasPosition() const override
has position index
ErrCode projectConstraint(const Constraint::Type type, const FitParams &fitparams, Projection &p) const override
find out which constraint it is and project
virtual ErrCode initParticleWithMother(FitParams &fitparams) override
init particle in case it has a mother
virtual int momIndex() const override
momentum index in fit params depending on whether it has a mother
ErrCode initMomentum(FitParams &fitparams) const
init momentum of *this and daughters
virtual void addToConstraintList(constraintlist &list, int depth) const override
add to constraint list
virtual ErrCode initMotherlessParticle(FitParams &fitparams) override
init particle in case it has no mother
static bool compTrkTransverseMomentum(const RecoTrack *lhs, const RecoTrack *rhs)
compare transverse track momentum
virtual int posIndex() const override
position index in fit params
InternalParticle(Belle2::Particle *particle, const ParticleBase *mother, const ConstraintConfiguration &config, bool forceFitAll)
constructor
bool m_lifetimeconstraint
has lifetime constraint
bool m_geo_constraint
use a geo metric constraint
bool m_shares_vertex_with_mother
shares vertex with mother, that means decay vertex = productionvertex
base class for all particles
Definition: ParticleBase.h:25
Belle2::Particle * particle() const
get basf2 particle
Definition: ParticleBase.h:96
ErrCode initTau(FitParams &par) const
initialises tau as a length
virtual ParticleBase * addDaughter(Belle2::Particle *, const ConstraintConfiguration &config, bool forceFitAll=false)
add daughter
virtual ErrCode projectMassConstraint(const FitParams &, Projection &) const
project mass constraint abstract
const ConstraintConfiguration * m_config
has all the constraint config
Definition: ParticleBase.h:224
virtual std::string parname(int index) const
get name of parameter i
virtual ErrCode projectConstraint(Constraint::Type, const FitParams &, Projection &) const
project constraint.
bool m_isStronglyDecayingResonance
decay length less than 1 micron
Definition: ParticleBase.h:221
virtual ErrCode initCovariance(FitParams &) const
init covariance matrix
Belle2::Particle * m_particle
pointer to framework type
Definition: ParticleBase.h:212
virtual ErrCode projectGeoConstraint(const FitParams &, Projection &) const
project geometrical constraint
virtual int posIndex() const
get vertex index (in statevector!)
Definition: ParticleBase.h:129
int index() const
get index
Definition: ParticleBase.h:99
std::vector< ParticleBase * > m_daughters
daughter container
Definition: ParticleBase.h:218
void collectVertexDaughters(std::vector< ParticleBase * > &particles, int posindex)
get vertex daughters
const ParticleBase * mother() const
getMother() / hasMother()
std::vector< Constraint > constraintlist
alias
Definition: ParticleBase.h:56
class to store the projected residuals and the corresponding jacobian as well as the covariance matri...
Definition: Projection.h:18
representation of all charged final states FIXME rename since this name is taken in tracking
Definition: RecoTrack.h:18
void setFlightLength(double flt)
setter for the flight length
Definition: RecoTrack.h:64