10#include <analysis/VertexFitting/TreeFitter/InternalParticle.h>
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>
18#include <framework/logging/Logger.h>
19#include <mdst/dataobjects/V0.h>
27 int lhstype = lhs->type() ;
28 int rhstype = rhs->type() ;
30 if (lhstype == rhstype &&
31 lhstype == ParticleBase::TFParticleType::kRecoTrack) {
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();
38 rc = lhstype < rhstype;
56 for (
auto daughter :
particle->getDaughters()) {
60 B2ERROR(
"Trying to create an InternalParticle from NULL. This should never happen.");
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");
73 [pdg = std::abs(
m_particle->getPDGCode())](
int val) { return std::abs(val) == pdg; }) != config.m_fixedToMotherVertexListPDG.end()
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()
96 status |= daughter->initMotherlessParticle(fitparams);
103 fitparams.
getStateVector().segment(posindex, 3) = Eigen::Matrix<double, 3, 1>::Zero(3);
105 std::vector<ParticleBase*> alldaughters;
108 std::vector<ParticleBase*> vtxdaughters;
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()
116 vtxdaughters.push_back(daughter);
120 if (trkdaughters.size() >= 2) {
123 auto dummy_vertex = ROOT::Math::XYZVector(0, 0, 0);
125 bool initWithV0 =
false;
126 if (v0 && v0->getFittedVertexPosition() != dummy_vertex) {
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(); });
135 if (recotrack_dau1 == trkdaughters.end() || recotrack_dau2 == trkdaughters.end()) {
136 B2WARNING(
"V0 daughter particles do not match with RecoTracks.");
138 double X_V0(v0->getFittedVertexX()), Y_V0(v0->getFittedVertexY()), Z_V0(v0->getFittedVertexZ());
143 Belle2::Helix helix1 = v0->getTrackFitResults().first->getHelix();
144 Belle2::Helix helix2 = v0->getTrackFitResults().second->getHelix();
146 (*recotrack_dau1)->setFlightLength(helix1.getArcLength2DAtXY(X_V0, Y_V0));
147 (*recotrack_dau2)->setFlightLength(helix2.getArcLength2DAtXY(X_V0, Y_V0));
162 double flt1(0), flt2(0);
172 }
else if (
false && trkdaughters.size() + vtxdaughters.size() >= 2) {
180 fitparams.
getStateVector().segment(posindex, 3) = Eigen::Matrix<double, 1, 3>::Zero(3);
185 daughter->initParticleWithMother(fitparams);
208 fitparams.
getStateVector().segment(momindex, 4) = Eigen::Matrix<double, 4, 1>::Zero(4);
211 int daumomindex = daughter->momIndex();
212 int maxrow = daughter->hasEnergy() ? 4 : 3;
214 double e2 = fitparams.
getStateVector().segment(daumomindex, maxrow).squaredNorm();
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);
225 return ErrCode(ErrCode::Status::success);
233 status |= daughter->initCovariance(fitparams);
245 p.getResiduals().segment(0, 4) = fitparams.
getStateVector().segment(momindex, 4);
247 for (
int imom = 0; imom < 4; ++imom) {
248 p.getH()(imom, momindex + imom) = 1;
252 const int daumomindex = daughter->momIndex();
253 const Eigen::Matrix<double, 1, 3> p3_vec = fitparams.
getStateVector().segment(daumomindex, 3);
256 p.getResiduals().segment(0, 3) -= p3_vec;
259 if (daughter->hasEnergy()) {
260 p.getResiduals()(3) -= fitparams.
getStateVector()(daumomindex + 3);
261 p.getH()(3, daumomindex + 3) = -1;
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;
274 for (
unsigned i = 0; i < 3; ++i) {
276 p.getH()(3, daumomindex + i) = -1 * p3_vec(i) / energy;
282 for (
unsigned i = 0; i < 3; ++i) {
283 p.getH()(i, daumomindex + i) = -1;
286 return ErrCode(ErrCode::Status::success);
296 const Eigen::Matrix<double, 4, 1> fitMomE = fitparams.
getStateVector().segment(momindex, 4);
298 p.getResiduals() =
m_config->m_beamMomE - fitMomE;
300 for (
int row = 0; row < 4; ++row) {
301 p.getH()(row, momindex + row) = -1;
304 p.getV() =
m_config->m_beamCovariance;
306 return ErrCode(ErrCode::Status::success) ;
316 case Constraint::mass:
319 case Constraint::geometric:
322 case Constraint::kinematic:
325 case Constraint::beam:
386 daughter->addToConstraintList(list, depth - 1);
389 list.push_back(
Constraint(
this, Constraint::lifetime, depth, 1));
392 list.push_back(
Constraint(
this, Constraint::kinematic, depth, 4, 3));
397 list.push_back(
Constraint(
this, Constraint::geometric, depth,
dim, 3));
400 list.push_back(
Constraint(
this, Constraint::mass, depth, 1, 3));
404 list.push_back(
Constraint(
this, Constraint::beam, depth, 4, 3));
412 if (!
mother() &&
id >= 3) {++id;}
419 daughter->forceP4Sum(fitparams);
Class to store reconstructed particles.
const V0 * getV0() const
Returns the pointer to the V0 object that was used to create this Particle (if ParticleType == c_V0).
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
const TrackFitResult * getTrackFitResult() const
Returns the pointer to the TrackFitResult that was used to create this Particle (ParticleType == c_Tr...
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Helix getHelix() const
Conversion to framework Helix (without covariance).
constraint configuration class
class to manage the order of constraints and their filtering
Type
type of constraints the order of these constraints is important: it is the order in which they are ap...
abstract errorocode be aware that the default is success
Class to store and manage fitparams (statevector)
Eigen::Matrix< double, -1, 1, 0, MAX_MATRIX_SIZE, 1 > & getStateVector()
getter for the fit parameters/statevector
int getDimensionOfState() const
get the states dimension
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
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
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...
representation of all charged final states FIXME rename since this name is taken in tracking
void setFlightLength(double flt)
setter for the flight length