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