Belle II Software light-2406-ragdoll
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#include <mdst/dataobjects/V0.h>
18
19using std::vector;
20
21namespace 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().Rho() > rhs->particle()->getMomentum().Rho();
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_beamconstraint(false),
49 m_lifetimeconstraint(false),
50 m_isconversion(false),
51 m_automatic_vertex_constraining(config.m_automatic_vertex_constraining)
52 {
53 if (particle) {
54 for (auto daughter : particle->getDaughters()) {
55 addDaughter(daughter, config, forceFitAll);
56 }
57 } else {
58 B2ERROR("Trying to create an InternalParticle from NULL. This should never happen.");
59 }
60
61 m_massconstraint = std::find(config.m_massConstraintListPDG.begin(), config.m_massConstraintListPDG.end(),
62 std::abs(m_particle->getPDGCode())) != config.m_massConstraintListPDG.end();
63
65
67 // if this is a hadronically decaying resonance it is useful to constraint the decay vertex to its mothers decay vertex.
68 //
71 std::abs(m_particle->getPDGCode())) != config.m_fixedToMotherVertexListPDG.end() && this->mother();
72
73 // use geo constraint if this particle is in the list to constrain
74 m_geo_constraint = std::find(config.m_geoConstraintListPDG.begin(),
75 config.m_geoConstraintListPDG.end(),
76 std::abs(m_particle->getPDGCode())) != config.m_geoConstraintListPDG.end() && this->mother() && !m_shares_vertex_with_mother;
77 } else {
80 }
81 }
82
84 {
85
86 return lhs->particle()->getMomentum().Rho() > rhs->particle()->getMomentum().Rho();
87 }
88
90 {
91 ErrCode status ;
92 for (auto daughter : m_daughters) {
93 status |= daughter->initMotherlessParticle(fitparams);
94 }
95
96 int posindex = posIndex();
97 // FIXME update this and check if this position is already initialised
98 // lower stream vertices might be better
99 if (hasPosition()) {
100 fitparams.getStateVector().segment(posindex, 3) = Eigen::Matrix<double, 3, 1>::Zero(3);
101
102 std::vector<ParticleBase*> alldaughters;
103 ParticleBase::collectVertexDaughters(alldaughters, posindex);
104
105 std::vector<ParticleBase*> vtxdaughters;
106
107 vector<RecoTrack*> trkdaughters;
108 for (auto daughter : alldaughters) {
109 if (daughter->type() == ParticleBase::TFParticleType::kRecoTrack) {
110 trkdaughters.push_back(static_cast<RecoTrack*>(daughter));
111 } else if (daughter->hasPosition()
112 && fitparams.getStateVector()(daughter->posIndex()) != 0) {
113 vtxdaughters.push_back(daughter);
114 }
115 }
116
117 if (trkdaughters.size() >= 2) {
118
119 auto v0 = particle()->getV0();
120 auto dummy_vertex = ROOT::Math::XYZVector(0, 0, 0);
121
122 bool initWithV0 = false;
123 if (v0 && v0->getFittedVertexPosition() != dummy_vertex) {
124 auto part_dau1 = particle()->getDaughter(0);
125 auto part_dau2 = particle()->getDaughter(1);
126
127 auto recotrack_dau1 = std::find_if(trkdaughters.begin(), trkdaughters.end(),
128 [&part_dau1](RecoTrack * rt) { return rt->particle()->getMdstArrayIndex() == part_dau1->getMdstArrayIndex(); });
129 auto recotrack_dau2 = std::find_if(trkdaughters.begin(), trkdaughters.end(),
130 [&part_dau2](RecoTrack * rt) { return rt->particle()->getMdstArrayIndex() == part_dau2->getMdstArrayIndex(); });
131
132 if (recotrack_dau1 == trkdaughters.end() || recotrack_dau2 == trkdaughters.end()) {
133 B2WARNING("V0 daughter particles do not match with RecoTracks.");
134 } else {
135 double X_V0(v0->getFittedVertexX()), Y_V0(v0->getFittedVertexY()), Z_V0(v0->getFittedVertexZ());
136 fitparams.getStateVector()(posindex) = X_V0;
137 fitparams.getStateVector()(posindex + 1) = Y_V0;
138 fitparams.getStateVector()(posindex + 2) = Z_V0;
139
140 Belle2::Helix helix1 = v0->getTrackFitResults().first->getHelix();
141 Belle2::Helix helix2 = v0->getTrackFitResults().second->getHelix();
142
143 (*recotrack_dau1)->setFlightLength(helix1.getArcLength2DAtXY(X_V0, Y_V0));
144 (*recotrack_dau2)->setFlightLength(helix2.getArcLength2DAtXY(X_V0, Y_V0));
145
146 initWithV0 = true;
147 }
148 }
149
150 if (!initWithV0) {
151 std::sort(trkdaughters.begin(), trkdaughters.end(), compTrkTransverseMomentum);
152
153 RecoTrack* dau1 = trkdaughters[0];
154 RecoTrack* dau2 = trkdaughters[1];
155
156 Belle2::Helix helix1 = dau1->particle()->getTrackFitResult()->getHelix();
157 Belle2::Helix helix2 = dau2->particle()->getTrackFitResult()->getHelix();
158
159 double flt1(0), flt2(0);
161 HelixUtils::helixPoca(helix1, helix2, flt1, flt2, v, m_isconversion);
162
163 fitparams.getStateVector()(posindex) = v.X();
164 fitparams.getStateVector()(posindex + 1) = v.Y();
165 fitparams.getStateVector()(posindex + 2) = v.Z();
166
167 dau1->setFlightLength(flt1);
168 dau2->setFlightLength(flt2);
169 }
170
171 } else if (false && trkdaughters.size() + vtxdaughters.size() >= 2) {
172 // TODO switched off waiting for refactoring of init1 and init2 functions (does not affect performance)
173 } else if (mother() && mother()->posIndex() >= 0) {
174 const int posindexmother = mother()->posIndex();
175 const int dim = m_config->m_originDimension; //TODO access mother
176 fitparams.getStateVector().segment(posindex, dim) = fitparams.getStateVector().segment(posindexmother, dim);
177 } else {
179 fitparams.getStateVector().segment(posindex, 3) = Eigen::Matrix<double, 1, 3>::Zero(3);
180 }
181 }
182
183 for (auto daughter : m_daughters) {
184 daughter->initParticleWithMother(fitparams);
185 }
186
187 initMomentum(fitparams);
188 return status;
189 }
190
192 {
193 int posindex = posIndex();
194 if (hasPosition() &&
195 mother() &&
196 fitparams.getStateVector()(posindex) == 0 &&
197 fitparams.getStateVector()(posindex + 1) == 0 && \
198 fitparams.getStateVector()(posindex + 2) == 0) {
199 const int posindexmom = mother()->posIndex();
200 const int dim = m_config->m_originDimension; //TODO access mother?
201 fitparams.getStateVector().segment(posindex, dim) = fitparams.getStateVector().segment(posindexmom, dim);
202 }
203 return initTau(fitparams);
204 }
205
207 {
208 int momindex = momIndex();
209 fitparams.getStateVector().segment(momindex, 4) = Eigen::Matrix<double, 4, 1>::Zero(4);
210
211 for (auto daughter : m_daughters) {
212 int daumomindex = daughter->momIndex();
213 int maxrow = daughter->hasEnergy() ? 4 : 3;
214
215 double e2 = fitparams.getStateVector().segment(daumomindex, maxrow).squaredNorm();
216 fitparams.getStateVector().segment(momindex, maxrow) += fitparams.getStateVector().segment(daumomindex, maxrow);
217
218 if (maxrow == 3) {
219 double mass = daughter->particle()->getPDGMass();
220 fitparams.getStateVector()(momindex + 3) += std::sqrt(e2 + mass * mass);
221 }
222 }
223 return ErrCode(ErrCode::Status::success);
224 }
225
227 {
228 ErrCode status;
230 for (auto daughter : m_daughters) {
231 status |= daughter->initCovariance(fitparams);
232 }
233 return status;
234 }
235
236
238 Projection& p) const
239 {
240 const int momindex = momIndex();
241
242 // `this` always has an energy row
243 p.getResiduals().segment(0, 4) = fitparams.getStateVector().segment(momindex, 4);
244
245 for (int imom = 0; imom < 4; ++imom) {
246 p.getH()(imom, momindex + imom) = 1;
247 }
248
249 for (const auto daughter : m_daughters) {
250 const int daumomindex = daughter->momIndex();
251 const Eigen::Matrix<double, 1, 3> p3_vec = fitparams.getStateVector().segment(daumomindex, 3);
252
253 // three momentum is easy just subtract the vectors
254 p.getResiduals().segment(0, 3) -= p3_vec;
255
256 // energy depends on the parametrisation!
257 if (daughter->hasEnergy()) {
258 p.getResiduals()(3) -= fitparams.getStateVector()(daumomindex + 3);
259 p.getH()(3, daumomindex + 3) = -1; // d/dE -E
260 } else {
261 // m^2 + p^2 = E^2
262 // so
263 // E = sqrt(m^2 + p^2)
264 const double mass = daughter->particle()->getPDGMass();
265 const double p2 = p3_vec.squaredNorm();
266 const double energy = std::sqrt(mass * mass + p2);
267 p.getResiduals()(3) -= energy;
268
269 for (unsigned i = 0; i < 3; ++i) {
270 // d/dpx_i sqrt(m^2 + p^2)
271 p.getH()(3, daumomindex + i) = -1 * p3_vec(i) / energy;
272 }
273 }
274
275 // this has to be in any case
276 // d/dp_i p_i
277 for (unsigned i = 0; i < 3; ++i) {
278 p.getH()(i, daumomindex + i) = -1;
279 }
280 }
281 return ErrCode(ErrCode::Status::success);
282 }
283
284
286 Projection& p) const
287 {
288
289 const int momindex = momIndex() ;
290
291 const Eigen::Matrix<double, 4, 1> fitMomE = fitparams.getStateVector().segment(momindex, 4);
292
293 p.getResiduals() = m_config->m_beamMomE - fitMomE;
294
295 for (int row = 0; row < 4; ++row) {
296 p.getH()(row, momindex + row) = -1;
297 }
298
299 p.getV() = m_config->m_beamCovariance;
300
301 return ErrCode(ErrCode::Status::success) ;
302 }
303
304
306 const FitParams& fitparams,
307 Projection& p) const
308 {
309 ErrCode status;
310 switch (type) {
311 case Constraint::mass:
312 status |= projectMassConstraint(fitparams, p);
313 break;
314 case Constraint::geometric:
315 status |= projectGeoConstraint(fitparams, p);
316 break;
317 case Constraint::kinematic:
318 status |= projectKineConstraint(fitparams, p);
319 break;
320 case Constraint::beam:
321 status |= projectBeamConstraint(fitparams, p);
322 break;
323 default:
324 status |= ParticleBase::projectConstraint(type, fitparams, p);
325 }
326
327 return status;
328 }
329
331 {
332 // for example B0 and D* can share the same vertex
333 return m_shares_vertex_with_mother ? this->mother()->posIndex() : index();
334 }
335
337 {
338 // { x, y, z, tau, px, py, pz, E }
339 // the last 4 always exists for composite particles
340 // tau index only exist with a vertex and geo constraint
341 //
342 if (m_shares_vertex_with_mother) { return 4; }
343 else if (!m_geo_constraint) { return 7; }
344 else { return 8; }
345 }
346
348 {
350 return m_geo_constraint ? index() + 3 : -1;
351 }
352
354 {
359 if (m_geo_constraint && !m_shares_vertex_with_mother) { return this->index() + 4; }
360
361 if (m_shares_vertex_with_mother) { return this->index(); }
362
363 if (!m_geo_constraint) {return index() + 3 ;}
364
365 // this will crash the initialisation
366 return -1;
367 }
368
370 {
371 // does this particle have a position index (decayvertex)
372 // true in any case
373 return true;
374 }
375
377 int depth) const
378 {
379
380 for (auto daughter : m_daughters) {
381 daughter->addToConstraintList(list, depth - 1);
382 }
383 if (tauIndex() >= 0 && m_lifetimeconstraint) {
384 list.push_back(Constraint(this, Constraint::lifetime, depth, 1));
385 }
386 if (momIndex() >= 0) {
387 list.push_back(Constraint(this, Constraint::kinematic, depth, 4, 3));
388 }
389 if (m_geo_constraint) {
390 assert(m_config);
391 const int dim = m_config->m_originDimension == 2 && std::abs(m_particle->getPDGCode()) == m_config->m_headOfTreePDG ? 2 : 3;
392 list.push_back(Constraint(this, Constraint::geometric, depth, dim, 3));
393 }
394 if (m_massconstraint) {
395 list.push_back(Constraint(this, Constraint::mass, depth, 1, 3));
396 }
397 if (m_beamconstraint) {
398 assert(m_config);
399 list.push_back(Constraint(this, Constraint::beam, depth, 4, 3));
400 }
401
402 }
403
404 std::string InternalParticle::parname(int thisindex) const
405 {
406 int id = thisindex;
407 if (!mother() && id >= 3) {++id;}
408 return ParticleBase::parname(id);
409 }
410
412 {
413 for (const auto daughter : m_daughters) {
414 daughter->forceP4Sum(fitparams);
415 }
416 const int momindex = momIndex();
417 if (momindex > 0) {
418 const int dim = hasEnergy() ? 4 : 3;
419 Projection p(fitparams.getDimensionOfState(), dim);
420 projectKineConstraint(fitparams, p);
421 fitparams.getStateVector().segment(momindex, dim) -= p.getResiduals().segment(0, dim);
422 }
423 }
424
425}
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:435
DataType X() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:431
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:433
This class represents an ideal helix in perigee parameterization.
Definition: Helix.h:59
double getArcLength2DAtXY(const double &x, const double &y) const
Calculates the two dimensional arc length at which the circle in the xy projection is closest to the ...
Definition: Helix.cc:141
Class to store reconstructed particles.
Definition: Particle.h:75
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:880
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:454
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:637
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:560
const TrackFitResult * getTrackFitResult() const
Returns the pointer to the TrackFitResult that was used to create this Particle (ParticleType == c_Tr...
Definition: Particle.cc:854
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition: Particle.cc:631
Helix getHelix() const
Conversion to framework Helix (without covariance).
constraint configuration class
const int m_originDimension
dimension of the origin constraint and ALL geometric gcosntraints
Eigen::Matrix< double, 4, 1 > m_beamMomE
Beam four-momentum.
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
const int m_beamConstraintPDG
PDG code to beam constraint.
Eigen::Matrix< double, 4, 4 > m_beamCovariance
Beam 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, Belle2::B2Vector3D &vertex, bool parallel=false)
POCA between two tracks.
Definition: HelixUtils.cc:214
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_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
Definition: ParticleBase.h:25
Belle2::Particle * particle() const
get basf2 particle
Definition: ParticleBase.h:92
ErrCode initTau(FitParams &par) const
initialises tau as a length
virtual ParticleBase * addDaughter(Belle2::Particle *, const ConstraintConfiguration &config, bool forceFitAll=false)
add daughter
Definition: ParticleBase.cc:65
virtual ErrCode projectMassConstraint(const FitParams &, Projection &) const
project mass constraint abstract
const ConstraintConfiguration * m_config
has all the constraint config
Definition: ParticleBase.h:204
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:201
virtual ErrCode initCovariance(FitParams &) const
init covariance matrix
Belle2::Particle * m_particle
pointer to framework type
Definition: ParticleBase.h:192
virtual ErrCode projectGeoConstraint(const FitParams &, Projection &) const
project geometrical constraint
virtual int posIndex() const
get vertex index (in statevector!)
Definition: ParticleBase.h:122
int index() const
get index
Definition: ParticleBase.h:95
std::vector< ParticleBase * > m_daughters
daughter container
Definition: ParticleBase.h:198
void collectVertexDaughters(std::vector< ParticleBase * > &particles, int posindex)
get vertex daughters
const ParticleBase * mother() const
getMother() / hasMother()
Definition: ParticleBase.h:98
std::vector< Constraint > constraintlist
alias
Definition: ParticleBase.h:52
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