Belle II Software  release-05-01-25
RecoKlong.cc
1 /**************************************************************************
2  *
3  * BASF2 (Belle Analysis Framework 2) *
4  * Copyright(C) 2018 - Belle II Collaboration *
5  * *
6  * Author: The Belle II Collaboration *
7  * Contributor: Wouter Hulsbergen, Francesco Tenchini, Jo-Frederik Krohn *
8  * *
9  * This software is provided "as is" without any warranty. *
10  **************************************************************************/
11 
12 #include <framework/logging/Logger.h>
13 
14 #include <analysis/dataobjects/Particle.h>
15 #include <mdst/dataobjects/KLMCluster.h>
16 
17 #include <analysis/VertexFitting/TreeFitter/RecoKlong.h>
18 #include <analysis/VertexFitting/TreeFitter/FitParams.h>
19 #include <analysis/VertexFitting/TreeFitter/ErrCode.h>
20 
21 namespace TreeFitter {
22 
23  RecoKlong::RecoKlong(Belle2::Particle* particle, const ParticleBase* mother) :
24  RecoParticle(particle, mother),
25  m_dim(3),
26  m_init(false),
27  m_useEnergy(true),
28  m_clusterPars(),
29  m_covariance()
30  {
31  initParams() ;
32  }
33 
35  {
36  const int posindexmother = mother()->posIndex();
37 
38  Eigen::Matrix<double, 1, 3> vertexToCluster = Eigen::Matrix<double, 1, 3>::Zero(1, 3);
39  for (unsigned int i = 0; i < 3; ++i) {
40  vertexToCluster(i) = m_clusterPars(i) - fitparams.getStateVector()(posindexmother + i);
41  }
42 
43  const double distanceToMother = vertexToCluster.norm();
44  const double energy = m_clusterPars(3);
45  const double energy2 = energy * energy;
46  const double pdgMass2 = ParticleBase::pdgMass() * ParticleBase::pdgMass();
47  const double absMom = -1 * std::sqrt(energy2 - pdgMass2);
48 
49  const int momindex = momIndex();
50 
51  for (unsigned int i = 0; i < 3; ++i) {
52  //px = |p| dx/|dx|
53  fitparams.getStateVector()(momindex + i) = absMom * vertexToCluster(i) / distanceToMother;
54  }
55 
56  fitparams.getStateVector()(momindex + 3) = energy ;
57 
58  return ErrCode(ErrCode::Status::success);
59  }
60 
61  ErrCode RecoKlong::initMotherlessParticle([[gnu::unused]] FitParams& fitparams)
62  {
63  return ErrCode(ErrCode::Status::success);
64  }
65 
66  ErrCode RecoKlong::initCovariance(FitParams& fitparams) const
67  {
68  const int momindex = momIndex();
69  const int posindex = mother()->posIndex();
70 
71  const double factorE = 1000 * m_covariance(3, 3);
72  const double factorX = 1000; // ~ 10cm error on initial vertex
73 
74  fitparams.getCovariance().block<4, 4>(momindex, momindex) =
75  Eigen::Matrix<double, 4, 4>::Identity(4, 4) * factorE;
76 
77  fitparams.getCovariance().block<3, 3>(posindex, posindex) =
78  Eigen::Matrix<double, 3, 3>::Identity(3, 3) * factorX;
79 
80  return ErrCode(ErrCode::Status::success);
81  }
82 
84  {
85  const Belle2::KLMCluster* cluster = particle()->getKLMCluster();
86 
87  const TVector3 cluster_pos = cluster->getClusterPosition();
88 
89  m_init = true;
90  m_covariance = Eigen::Matrix<double, 4, 4>::Zero(4, 4);
91 
92  TMatrixDSym cov7 = cluster->getError7x7();
93 
94  for (int row = 0; row < 3; ++row) {
95  for (int col = 0; col < 3; ++col) {
96  m_covariance(row, col) = cov7[row + 4][col + 4] ;
97  }
98  }
99 
103  if (0 == m_covariance(3, 3)) {m_covariance(3, 3) = .214;}
104 
105  m_clusterPars(0) = cluster_pos.X();
106  m_clusterPars(1) = cluster_pos.Y();
107  m_clusterPars(2) = cluster_pos.Z();
108  m_clusterPars(3) = sqrt(particle()->getPDGMass() * particle()->getPDGMass() + cluster->getMomentumMag() *
109  cluster->getMomentumMag());
110 
111  auto p_vec = particle()->getMomentum();
112  // find highest momentum, eliminate dim with highest mom
113  if ((std::abs(p_vec(0)) >= std::abs(p_vec(1))) && (std::abs(p_vec(0)) >= std::abs(p_vec(2)))) {
114  m_i1 = 0; m_i2 = 1; m_i3 = 2;
115  } else if ((std::abs(p_vec(1)) >= std::abs(p_vec(0))) && (std::abs(p_vec(1)) >= std::abs(p_vec(2)))) {
116  m_i1 = 1; m_i2 = 0; m_i3 = 2;
117  } else if ((std::abs(p_vec(2)) >= std::abs(p_vec(1))) && (std::abs(p_vec(2)) >= std::abs(p_vec(0)))) {
118  m_i1 = 2; m_i2 = 1; m_i3 = 0;
119  } else {
120  B2ERROR("Could not estimate highest momentum for Klong constraint. Aborting this fit.\n px: "
121  << p_vec(0) << " py: " << p_vec(1) << " pz: " << p_vec(2) << " calculated from Ec: " << m_clusterPars(3));
122  return ErrCode(ErrCode::Status::photondimerror);
123  }
124 
125  return ErrCode(ErrCode::Status::success);
126  }
127 
128  ErrCode RecoKlong::projectRecoConstraint(const FitParams& fitparams, Projection& p) const
129  {
130  ErrCode status ;
131  const int momindex = momIndex() ;
132  const int posindex = mother()->posIndex();
133 
134  const Eigen::Matrix<double, 1, 3> x_vertex = fitparams.getStateVector().segment(posindex, 3);
135  const Eigen::Matrix<double, 1, 3> p_vec = fitparams.getStateVector().segment(momindex, 3);
136 
137  if (0 == p_vec[m_i1]) { return ErrCode(ErrCode::Status::klongdimerror); }
138 
139  // p_vec[m_i1] must not be 0
140  const double elim = (m_clusterPars[m_i1] - x_vertex[m_i1]) / p_vec[m_i1];
141  const double mom = p_vec.norm();
142  const double energy = fitparams.getStateVector()(momindex + 3);
143 
144  // r'
145  Eigen::Matrix<double, 3, 1> residual3 = Eigen::Matrix<double, 3, 1>::Zero(3, 1);
146  residual3(0) = m_clusterPars[m_i2] - x_vertex[m_i2] - p_vec[m_i2] * elim;
147  residual3(1) = m_clusterPars[m_i3] - x_vertex[m_i3] - p_vec[m_i3] * elim;
148  residual3(2) = m_clusterPars[3] - energy;
149 
150  Eigen::Matrix<double, 3, 4> P = Eigen::Matrix<double, 3, 4>::Zero(3, 4);
151 
152  // dr'/dm | m:={xc,yc,zc,Ec} the measured quantities
153  P(0, m_i2) = 1;
154  P(0, m_i1) = - p_vec[m_i2] / p_vec[m_i1];
155 
156  P(1, m_i3) = 1;
157  P(1, m_i1) = - p_vec[m_i3] / p_vec[m_i1];
158  P(2, 3) = 1; // dE/dEc
159 
160  p.getResiduals().segment(0, 3) = residual3;
161 
162  p.getV() = P * m_covariance.selfadjointView<Eigen::Lower>() * P.transpose();
163 
164  // dr'/dm | m:={x,y,z,px,py,pz,E}
165  // x := x_vertex (decay vertex of mother)
166  p.getH()(0, posindex + m_i1) = p_vec[m_i2] / p_vec[m_i1];
167  p.getH()(0, posindex + m_i2) = -1.0;
168  p.getH()(0, posindex + m_i3) = 0;
169 
170  p.getH()(1, posindex + m_i1) = p_vec[m_i3] / p_vec[m_i1];
171  p.getH()(1, posindex + m_i2) = 0;
172  p.getH()(1, posindex + m_i3) = -1.0;
173 
174  p.getH()(0, momindex + m_i1) = p_vec[m_i2] * elim / p_vec[m_i1];
175  p.getH()(0, momindex + m_i2) = -1. * elim;
176  p.getH()(0, momindex + m_i3) = 0;
177 
178  p.getH()(1, momindex + m_i1) = p_vec[m_i3] * elim / p_vec[m_i1];;
179  p.getH()(1, momindex + m_i2) = 0;
180  p.getH()(1, momindex + m_i3) = -1. * elim;
181 
182  p.getH()(2, momindex + m_i1) = -1. * p_vec[m_i1] / mom;
183  p.getH()(2, momindex + m_i2) = -1. * p_vec[m_i2] / mom;
184  p.getH()(2, momindex + m_i3) = -1. * p_vec[m_i3] / mom;
185  p.getH()(2, momindex + 3) = -1;
186 
187  return ErrCode(ErrCode::Status::success);
188  }
189 
190 }
191 
192 
TreeFitter::RecoKlong::initParams
ErrCode initParams()
update or init params
Definition: RecoKlong.cc:92
TreeFitter::RecoKlong::m_clusterPars
Eigen::Matrix< double, 1, 4 > m_clusterPars
constrains measured params (x_c, y_c, z_c, E_c)
Definition: RecoKlong.h:90
TreeFitter::ErrCode
abstract errorocode be aware that the default is succes
Definition: ErrCode.h:23
TreeFitter::RecoKlong::initMotherlessParticle
virtual ErrCode initMotherlessParticle(FitParams &fitparams) override
init particle without mother
Definition: RecoKlong.cc:70
TreeFitter::ParticleBase::pdgMass
double pdgMass() const
get pdg mass
Definition: ParticleBase.h:164
TreeFitter::RecoKlong::projectRecoConstraint
ErrCode projectRecoConstraint(const FitParams &fitparams, Projection &p) const override
project klong constraint
Definition: RecoKlong.cc:137
TreeFitter::FitParams
Class to store and manage fitparams (statevector)
Definition: FitParams.h:29
TreeFitter::RecoKlong::RecoKlong
RecoKlong(Belle2::Particle *bc, const ParticleBase *mother)
constructor
Definition: RecoKlong.cc:32
TreeFitter::RecoKlong::initParticleWithMother
virtual ErrCode initParticleWithMother(FitParams &fitparams) override
init particle with mother
Definition: RecoKlong.cc:43
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::RecoKlong::m_i2
int m_i2
random index
Definition: RecoKlong.h:98
Belle2::Particle::getKLMCluster
const KLMCluster * getKLMCluster() const
Returns the pointer to the KLMCluster object that was used to create this Particle (ParticleType == c...
Definition: Particle.cc:850
TreeFitter::RecoKlong::initCovariance
ErrCode initCovariance(FitParams &fitparams) const override
init covariance
Definition: RecoKlong.cc:75
TreeFitter::RecoKlong::m_i3
int m_i3
another random index
Definition: RecoKlong.h:100
TreeFitter::ParticleBase::posIndex
virtual int posIndex() const
get vertex index (in statevector!)
Definition: ParticleBase.h:142
Belle2::KLMCluster
KLM cluster data.
Definition: KLMCluster.h:38
Belle2::Particle::getMomentum
TVector3 getMomentum() const
Returns momentum vector.
Definition: Particle.h:475
TreeFitter::RecoParticle::momIndex
virtual int momIndex() const override
get momentum index
Definition: RecoParticle.h:59
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
TreeFitter::RecoKlong::m_i1
int m_i1
index with the highest momentum.
Definition: RecoKlong.h:96
TreeFitter::ParticleBase::mother
const ParticleBase * mother() const
getMother() / hasMother()
Definition: ParticleBase.cc:295
TreeFitter::RecoKlong::m_covariance
Eigen::Matrix< double, 4, 4 > m_covariance
covariance (x_c,y_c,z_c,E_c) of measured pars
Definition: RecoKlong.h:93
TreeFitter::RecoKlong::m_init
bool m_init
was initialized*
Definition: RecoKlong.h:84
TreeFitter::ParticleBase::particle
Belle2::Particle * particle() const
get basf2 particle
Definition: ParticleBase.h:109