Belle II Software development
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/VertexFitting/TreeFitter/InternalParticle.h>
11
12#include <analysis/VertexFitting/TreeFitter/ConstraintConfiguration.h>
13#include <analysis/VertexFitting/TreeFitter/FitParams.h>
14#include <analysis/VertexFitting/TreeFitter/HelixUtils.h>
15#include <analysis/VertexFitting/TreeFitter/Projection.h>
16#include <analysis/VertexFitting/TreeFitter/RecoTrack.h>
17
18#include <framework/logging/Logger.h>
19#include <mdst/dataobjects/V0.h>
20
21using std::vector;
22
23namespace TreeFitter {
24
25 inline bool sortByType(const ParticleBase* lhs, const ParticleBase* rhs)
26 {
27 int lhstype = lhs->type() ;
28 int rhstype = rhs->type() ;
29 bool rc = false ;
30 if (lhstype == rhstype &&
31 lhstype == ParticleBase::TFParticleType::kRecoTrack) {
32
33 rc = lhs->particle()->getMomentum().Rho() > rhs->particle()->getMomentum().Rho();
34 } else if (lhs->particle() && rhs->particle() && lhs->particle()->getNDaughters() > 0 &&
35 rhs->particle()->getNDaughters() > 0) {
36 rc = lhs->nFinalChargedCandidates() > rhs->nFinalChargedCandidates();
37 } else {
38 rc = lhstype < rhstype;
39 }
40 return rc;
41 }
42
44 const ParticleBase* mother,
45 const ConstraintConfiguration& config,
46 bool forceFitAll
47 ) :
48 ParticleBase(particle, mother, &config),// config pointer here to allow final states not to have it
49 m_massconstraint(false),
50 m_beamconstraint(false),
52 m_isconversion(false),
54 {
55 if (particle) {
56 for (auto daughter : particle->getDaughters()) {
57 addDaughter(daughter, config, forceFitAll);
58 }
59 } else {
60 B2ERROR("Trying to create an InternalParticle from NULL. This should never happen.");
61 }
62
63 m_massconstraint = std::find_if(config.m_massConstraintListPDG.begin(), config.m_massConstraintListPDG.end(),
64 [pdg = std::abs(m_particle->getPDGCode())](int val) { return std::abs(val) == pdg; }) != config.m_massConstraintListPDG.end()
65 or m_particle->hasExtraInfo("treeFitterMassConstraint");
66
67 m_beamconstraint = (std::abs(m_particle->getPDGCode()) == std::abs(config.m_beamConstraintPDG));
68
70 // if this is a hadronically decaying resonance it is useful to constrain the decay vertex to its mother's decay vertex.
71 //
72 m_shares_vertex_with_mother = std::find_if(config.m_fixedToMotherVertexListPDG.begin(), config.m_fixedToMotherVertexListPDG.end(),
73 [pdg = std::abs(m_particle->getPDGCode())](int val) { return std::abs(val) == pdg; }) != config.m_fixedToMotherVertexListPDG.end()
74 and this->mother();
75
76 // use geo constraint if this particle is in the list to constrain
77 m_geo_constraint = std::find_if(config.m_geoConstraintListPDG.begin(), config.m_geoConstraintListPDG.end(),
78 [pdg = std::abs(m_particle->getPDGCode())](int val) { return std::abs(val) == pdg; }) != config.m_geoConstraintListPDG.end()
79 and this->mother() and !m_shares_vertex_with_mother;
80 } else {
83 }
84 }
85
87 {
88
89 return lhs->particle()->getMomentum().Rho() > rhs->particle()->getMomentum().Rho();
90 }
91
93 {
94 ErrCode status ;
95 for (auto daughter : m_daughters) {
96 status |= daughter->initMotherlessParticle(fitparams);
97 }
98
99 int posindex = posIndex();
100 // FIXME update this and check if this position is already initialised
101 // lower stream vertices might be better
102 if (hasPosition()) {
103 fitparams.getStateVector().segment(posindex, 3) = Eigen::Matrix<double, 3, 1>::Zero(3);
104
105 std::vector<ParticleBase*> alldaughters;
106 ParticleBase::collectVertexDaughters(alldaughters, posindex);
107
108 std::vector<ParticleBase*> vtxdaughters;
109
110 vector<RecoTrack*> trkdaughters;
111 for (auto daughter : alldaughters) {
112 if (daughter->type() == ParticleBase::TFParticleType::kRecoTrack) {
113 trkdaughters.push_back(static_cast<RecoTrack*>(daughter));
114 } else if (daughter->hasPosition()
115 && fitparams.getStateVector()(daughter->posIndex()) != 0) {
116 vtxdaughters.push_back(daughter);
117 }
118 }
119
120 if (trkdaughters.size() >= 2) {
121
122 auto v0 = particle()->getV0();
123 auto dummy_vertex = ROOT::Math::XYZVector(0, 0, 0);
124
125 bool initWithV0 = false;
126 if (v0 && v0->getFittedVertexPosition() != dummy_vertex) {
127 auto part_dau1 = particle()->getDaughter(0);
128 auto part_dau2 = particle()->getDaughter(1);
129
130 auto recotrack_dau1 = std::find_if(trkdaughters.begin(), trkdaughters.end(),
131 [&part_dau1](RecoTrack * rt) { return rt->particle()->getMdstSource() == part_dau1->getMdstSource(); });
132 auto recotrack_dau2 = std::find_if(trkdaughters.begin(), trkdaughters.end(),
133 [&part_dau2](RecoTrack * rt) { return rt->particle()->getMdstSource() == part_dau2->getMdstSource(); });
134
135 if (recotrack_dau1 == trkdaughters.end() || recotrack_dau2 == trkdaughters.end()) {
136 B2WARNING("V0 daughter particles do not match with RecoTracks.");
137 } else {
138 double X_V0(v0->getFittedVertexX()), Y_V0(v0->getFittedVertexY()), Z_V0(v0->getFittedVertexZ());
139 fitparams.getStateVector()(posindex) = X_V0;
140 fitparams.getStateVector()(posindex + 1) = Y_V0;
141 fitparams.getStateVector()(posindex + 2) = Z_V0;
142
143 Belle2::Helix helix1 = v0->getTrackFitResults().first->getHelix();
144 Belle2::Helix helix2 = v0->getTrackFitResults().second->getHelix();
145
146 (*recotrack_dau1)->setFlightLength(helix1.getArcLength2DAtXY(X_V0, Y_V0));
147 (*recotrack_dau2)->setFlightLength(helix2.getArcLength2DAtXY(X_V0, Y_V0));
148
149 initWithV0 = true;
150 }
151 }
152
153 if (!initWithV0) {
154 std::sort(trkdaughters.begin(), trkdaughters.end(), compTrkTransverseMomentum);
155
156 RecoTrack* dau1 = trkdaughters[0];
157 RecoTrack* dau2 = trkdaughters[1];
158
159 Belle2::Helix helix1 = dau1->particle()->getTrackFitResult()->getHelix();
160 Belle2::Helix helix2 = dau2->particle()->getTrackFitResult()->getHelix();
161
162 double flt1(0), flt2(0);
163 Eigen::Vector3d v;
164 HelixUtils::helixPoca(helix1, helix2, flt1, flt2, v, m_isconversion);
165
166 fitparams.getStateVector().segment<3>(posindex) = v;
167
168 dau1->setFlightLength(flt1);
169 dau2->setFlightLength(flt2);
170 }
171
172 } else if (false && trkdaughters.size() + vtxdaughters.size() >= 2) {
173 // TODO switched off waiting for refactoring of init1 and init2 functions (does not affect performance)
174 } else if (mother() && mother()->posIndex() >= 0) {
175 const int posindexmother = mother()->posIndex();
176 const int dim = m_config->m_originDimension; //TODO access mother
177 fitparams.getStateVector().segment(posindex, dim) = fitparams.getStateVector().segment(posindexmother, dim);
178 } else {
180 fitparams.getStateVector().segment(posindex, 3) = Eigen::Matrix<double, 1, 3>::Zero(3);
181 }
182 }
183
184 for (auto daughter : m_daughters) {
185 daughter->initParticleWithMother(fitparams);
186 }
187
188 initMomentum(fitparams);
189 return status;
190 }
191
193 {
194 int posindex = posIndex();
195 if (hasPosition() &&
196 mother() &&
197 fitparams.getStateVector().segment<3>(posindex).isZero()) {
198 const int posindexmom = mother()->posIndex();
199 const int dim = m_config->m_originDimension; //TODO access mother?
200 fitparams.getStateVector().segment(posindex, dim) = fitparams.getStateVector().segment(posindexmom, dim);
201 }
202 return initTau(fitparams);
203 }
204
206 {
207 int momindex = momIndex();
208 fitparams.getStateVector().segment(momindex, 4) = Eigen::Matrix<double, 4, 1>::Zero(4);
209
210 for (auto daughter : m_daughters) {
211 int daumomindex = daughter->momIndex();
212 int maxrow = daughter->hasEnergy() ? 4 : 3;
213
214 double e2 = fitparams.getStateVector().segment(daumomindex, maxrow).squaredNorm();
215 fitparams.getStateVector().segment(momindex, maxrow) += fitparams.getStateVector().segment(daumomindex, maxrow);
216
217 if (maxrow == 3) {
218 double mass = 0;
219 if (daughter->particle()->hasExtraInfo("treeFitterMassConstraintValue")) {
220 mass = daughter->particle()->getExtraInfo("treeFitterMassConstraintValue");
221 } else mass = daughter->particle()->getPDGMass();
222 fitparams.getStateVector()(momindex + 3) += std::sqrt(e2 + mass * mass);
223 }
224 }
225 return ErrCode(ErrCode::Status::success);
226 }
227
229 {
230 ErrCode status;
232 for (auto daughter : m_daughters) {
233 status |= daughter->initCovariance(fitparams);
234 }
235 return status;
236 }
237
238
240 Projection& p) const
241 {
242 const int momindex = momIndex();
243
244 // `this` always has an energy row
245 p.getResiduals().segment(0, 4) = fitparams.getStateVector().segment(momindex, 4);
246
247 for (int imom = 0; imom < 4; ++imom) {
248 p.getH()(imom, momindex + imom) = 1;
249 }
250
251 for (const auto daughter : m_daughters) {
252 const int daumomindex = daughter->momIndex();
253 const Eigen::Matrix<double, 1, 3> p3_vec = fitparams.getStateVector().segment(daumomindex, 3);
254
255 // three momentum is easy just subtract the vectors
256 p.getResiduals().segment(0, 3) -= p3_vec;
257
258 // energy depends on the parametrisation!
259 if (daughter->hasEnergy()) {
260 p.getResiduals()(3) -= fitparams.getStateVector()(daumomindex + 3);
261 p.getH()(3, daumomindex + 3) = -1; // d/dE -E
262 } else {
263 // m^2 + p^2 = E^2
264 // so
265 // E = sqrt(m^2 + p^2)
266 double mass = 0;
267 if (daughter->particle()->hasExtraInfo("treeFitterMassConstraintValue")) {
268 mass = daughter->particle()->getExtraInfo("treeFitterMassConstraintValue");
269 } else mass = daughter->particle()->getPDGMass();
270 const double p2 = p3_vec.squaredNorm();
271 const double energy = std::sqrt(mass * mass + p2);
272 p.getResiduals()(3) -= energy;
273
274 for (unsigned i = 0; i < 3; ++i) {
275 // d/dpx_i sqrt(m^2 + p^2)
276 p.getH()(3, daumomindex + i) = -1 * p3_vec(i) / energy;
277 }
278 }
279
280 // this has to be in any case
281 // d/dp_i p_i
282 for (unsigned i = 0; i < 3; ++i) {
283 p.getH()(i, daumomindex + i) = -1;
284 }
285 }
286 return ErrCode(ErrCode::Status::success);
287 }
288
289
291 Projection& p) const
292 {
293
294 const int momindex = momIndex() ;
295
296 const Eigen::Matrix<double, 4, 1> fitMomE = fitparams.getStateVector().segment(momindex, 4);
297
298 p.getResiduals() = m_config->m_beamMomE - fitMomE;
299
300 for (int row = 0; row < 4; ++row) {
301 p.getH()(row, momindex + row) = -1;
302 }
303
304 p.getV() = m_config->m_beamCovariance;
305
306 return ErrCode(ErrCode::Status::success) ;
307 }
308
309
311 const FitParams& fitparams,
312 Projection& p) const
313 {
314 ErrCode status;
315 switch (type) {
316 case Constraint::mass:
317 status |= projectMassConstraint(fitparams, p);
318 break;
319 case Constraint::geometric:
320 status |= projectGeoConstraint(fitparams, p);
321 break;
322 case Constraint::kinematic:
323 status |= projectKineConstraint(fitparams, p);
324 break;
325 case Constraint::beam:
326 status |= projectBeamConstraint(fitparams, p);
327 break;
328 default:
329 status |= ParticleBase::projectConstraint(type, fitparams, p);
330 }
331
332 return status;
333 }
334
336 {
337 // for example B0 and D* can share the same vertex
338 return m_shares_vertex_with_mother ? this->mother()->posIndex() : index();
339 }
340
342 {
343 // { x, y, z, tau, px, py, pz, E }
344 // the last 4 always exists for composite particles
345 // tau index only exist with a vertex and geo constraint
346 //
347 if (m_shares_vertex_with_mother) { return 4; }
348 else if (!m_geo_constraint) { return 7; }
349 else { return 8; }
350 }
351
353 {
355 return m_geo_constraint ? index() + 3 : -1;
356 }
357
359 {
363
364 if (m_geo_constraint && !m_shares_vertex_with_mother) { return this->index() + 4; }
365
366 if (m_shares_vertex_with_mother) { return this->index(); }
367
368 if (!m_geo_constraint) {return index() + 3 ;}
369
370 // this will crash the initialisation
371 return -1;
372 }
373
375 {
376 // does this particle have a position index (decayvertex)
377 // true in any case
378 return true;
379 }
380
382 int depth) const
383 {
384
385 for (auto daughter : m_daughters) {
386 daughter->addToConstraintList(list, depth - 1);
387 }
388 if (tauIndex() >= 0 && m_lifetimeconstraint) {
389 list.push_back(Constraint(this, Constraint::lifetime, depth, 1));
390 }
391 if (momIndex() >= 0) {
392 list.push_back(Constraint(this, Constraint::kinematic, depth, 4, 3));
393 }
394 if (m_geo_constraint) {
395 assert(m_config);
396 const int dim = m_config->m_originDimension == 2 && std::abs(m_particle->getPDGCode()) == m_config->m_headOfTreePDG ? 2 : 3;
397 list.push_back(Constraint(this, Constraint::geometric, depth, dim, 3));
398 }
399 if (m_massconstraint) {
400 list.push_back(Constraint(this, Constraint::mass, depth, 1, 3));
401 }
402 if (m_beamconstraint) {
403 assert(m_config);
404 list.push_back(Constraint(this, Constraint::beam, depth, 4, 3));
405 }
406
407 }
408
409 std::string InternalParticle::parname(int thisindex) const
410 {
411 int id = thisindex;
412 if (!mother() && id >= 3) {++id;}
413 return ParticleBase::parname(id);
414 }
415
417 {
418 for (const auto daughter : m_daughters) {
419 daughter->forceP4Sum(fitparams);
420 }
421 const int momindex = momIndex();
422 if (momindex > 0) {
423 const int dim = hasEnergy() ? 4 : 3;
424 Projection p(fitparams.getDimensionOfState(), dim);
425 projectKineConstraint(fitparams, p);
426 fitparams.getStateVector().segment(momindex, dim) -= p.getResiduals().segment(0, dim);
427 }
428 }
429
430}
Helix parameter class.
Definition Helix.h:48
Class to store reconstructed particles.
Definition Particle.h:76
const V0 * getV0() const
Returns the pointer to the V0 object that was used to create this Particle (if ParticleType == c_V0).
Definition Particle.cc:956
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition Particle.h:580
const TrackFitResult * getTrackFitResult() const
Returns the pointer to the TrackFitResult that was used to create this Particle (ParticleType == c_Tr...
Definition Particle.cc:925
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition Particle.cc:662
Helix getHelix() const
Conversion to framework Helix (without covariance).
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:80
static double helixPoca(const Belle2::Helix &helix1, const Belle2::Helix &helix2, double &flt1, double &flt2, Eigen::Vector3d &vertex, bool parallel=false)
POCA between two tracks.
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
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_beamconstraint
has beam constraint
bool m_massconstraint
has mass constraint
ErrCode projectBeamConstraint(const FitParams &, Projection &) const
project beam four momentum 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
Belle2::Particle * particle() const
get basf2 particle
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
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
virtual ErrCode initCovariance(FitParams &) const
init covariance matrix
Belle2::Particle * m_particle
pointer to framework type
ParticleBase(Belle2::Particle *particle, const ParticleBase *mother, const ConstraintConfiguration *config=nullptr)
default constructor
virtual ErrCode projectGeoConstraint(const FitParams &, Projection &) const
project geometrical constraint
virtual int posIndex() const
get vertex index (in statevector!)
int index() const
get index
std::vector< ParticleBase * > m_daughters
daughter container
void collectVertexDaughters(std::vector< ParticleBase * > &particles, int posindex)
get vertex daughters
const ParticleBase * mother() const
getMother() / hasMother()
std::vector< Constraint > constraintlist
alias
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