Belle II Software light-2406-ragdoll
RecoPhoton.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 <framework/logging/Logger.h>
11
12#include <analysis/dataobjects/Particle.h>
13#include <mdst/dataobjects/ECLCluster.h>
14
15#include <analysis/ClusterUtility/ClusterUtils.h>
16
17#include <analysis/VertexFitting/TreeFitter/RecoPhoton.h>
18#include <analysis/VertexFitting/TreeFitter/FitParams.h>
19#include <analysis/VertexFitting/TreeFitter/ErrCode.h>
20
21#include <framework/gearbox/Const.h>
22#include <framework/geometry/B2Vector3.h>
23
24namespace TreeFitter {
25
26 RecoPhoton::RecoPhoton(Belle2::Particle* particle, const ParticleBase* mother) : RecoParticle(particle, mother),
27 m_dim(3),
28 m_init(false),
29 m_useEnergy(useEnergy(*particle)),
30 m_clusterPars(),
31 m_covariance(),
32 m_momentumScalingFactor(particle->getEffectiveMomentumScale())
33 {
34 initParams();
35 }
36
38 {
39 const int posindexmother = mother()->posIndex();
40
41 Eigen::Matrix<double, 1, 3> vertexToCluster = Eigen::Matrix<double, 1, 3>::Zero(1, 3);
42 for (unsigned int i = 0; i < 3; ++i) {
43 vertexToCluster(i) = m_clusterPars(i) - fitparams.getStateVector()(posindexmother + i);
44 }
45
46 const double distanceToMother = vertexToCluster.norm();
47 const double energy = m_momentumScalingFactor * m_clusterPars(3); // apply scaling factor to correct energy bias
48 const int momindex = momIndex();
49
50 for (unsigned int i = 0; i < 3; ++i) {
51 // px = E dx/|dx|
52 fitparams.getStateVector()(momindex + i) = energy * vertexToCluster(i) / distanceToMother;
53 }
54
55 return ErrCode(ErrCode::Status::success);
56 }
57
59 {
60 return ErrCode(ErrCode::Status::success);
61 }
62
64 {
65 bool rc = true;
66 const int pdg = particle.getPDGCode();
67 if (pdg &&
70 rc = false;
71 }
72 return rc;
73 }
74
76 {
77 const int momindex = momIndex();
78 const int posindex = mother()->posIndex();
79
80 const double factorE = 1000 * m_covariance(3, 3);
81 const double factorX = 1000; // ~ 10cm error on initial vertex
82
83 fitparams.getCovariance().block<4, 4>(momindex, momindex) =
84 Eigen::Matrix<double, 4, 4>::Identity(4, 4) * factorE;
85
86 fitparams.getCovariance().block<3, 3>(posindex, posindex) =
87 Eigen::Matrix<double, 3, 3>::Identity(3, 3) * factorX;
88
89 return ErrCode(ErrCode::Status::success);
90 }
91
93 {
94 const Belle2::ECLCluster* cluster = particle()->getECLCluster();
96 const Belle2::B2Vector3D centroid = cluster->getClusterPosition();
97 const double energy = cluster->getEnergy(clusterhypo);
98
99 m_init = true;
100 m_covariance = Eigen::Matrix<double, 4, 4>::Zero(4, 4);
102
103 TMatrixDSym cov_EPhiTheta = C.GetCovarianceMatrix3x3FromCluster(cluster);
104
105 Eigen::Matrix<double, 2, 2> covPhiTheta = Eigen::Matrix<double, 2, 2>::Zero(2, 2);
106
107 for (int row = 0; row < 2; ++row) {
108 // we go through all elements here instead of selfadjoint view later
109 for (int col = 0; col < 2; ++col) {
110 covPhiTheta(row, col) = cov_EPhiTheta[row + 1][col + 1];
111 }
112 }
113
114 // the in going x-E correlations are 0 so we don't fill them
115 const double R = cluster->getR();
116 const double theta = cluster->getPhi();
117 const double phi = cluster->getTheta();
118
119 const double st = std::sin(theta);
120 const double ct = std::cos(theta);
121 const double sp = std::sin(phi);
122 const double cp = std::cos(phi);
123
124 Eigen::Matrix<double, 2, 3> polarToCartesian = Eigen::Matrix<double, 2, 3>::Zero(2, 3);
125
126 // polarToCartesian({phi,theta} -> {x,y,z} )
127 polarToCartesian(0, 0) = -1. * R * st * sp; // dx/dphi
128 polarToCartesian(0, 1) = R * st * cp; // dy/dphi
129 polarToCartesian(0, 2) = 0; // dz/dphi
130
131 polarToCartesian(1, 0) = R * ct * cp; // dx/dtheta
132 polarToCartesian(1, 1) = R * ct * sp; // dy/dtheta
133 polarToCartesian(1, 2) = -1. * R * st; // dz/dtheta
134
135 m_covariance.block<3, 3>(0, 0) = polarToCartesian.transpose() * covPhiTheta * polarToCartesian;
136
137 m_covariance(3, 3) = cov_EPhiTheta[0][0];
138 m_clusterPars(0) = centroid.X();
139 m_clusterPars(1) = centroid.Y();
140 m_clusterPars(2) = centroid.Z();
141 m_clusterPars(3) = energy;
142
143 auto p_vec = particle()->getMomentum();
144 // find highest momentum, eliminate dim with highest mom
145 if ((std::abs(p_vec.X()) >= std::abs(p_vec.Y())) && (std::abs(p_vec.X()) >= std::abs(p_vec.Z()))) {
146 m_i1 = 0;
147 m_i2 = 1;
148 m_i3 = 2;
149 } else if ((std::abs(p_vec.Y()) >= std::abs(p_vec.X())) && (std::abs(p_vec.Y()) >= std::abs(p_vec.Z()))) {
150 m_i1 = 1;
151 m_i2 = 0;
152 m_i3 = 2;
153 } else if ((std::abs(p_vec.Z()) >= std::abs(p_vec.Y())) && (std::abs(p_vec.Z()) >= std::abs(p_vec.X()))) {
154 m_i1 = 2;
155 m_i2 = 1;
156 m_i3 = 0;
157 } else {
158 B2ERROR("Could not estimate highest momentum for photon constraint. Aborting this fit.\n px: "
159 << p_vec.X() << " py: " << p_vec.Y() << " pz: " << p_vec.Z() << " calculated from Ec: " << m_clusterPars(3));
160 return ErrCode(ErrCode::Status::photondimerror);
161 }
162
163 return ErrCode(ErrCode::Status::success);
164 }
165
167 {
168 const int momindex = momIndex();
169 const int posindex = mother()->posIndex();
186 const Eigen::Matrix<double, 1, 3> x_vertex = fitparams.getStateVector().segment(posindex, 3);
187 const Eigen::Matrix<double, 1, 3> p_vec = fitparams.getStateVector().segment(momindex, 3);
188
189 if (0 == p_vec[m_i1]) {
190 return ErrCode(ErrCode::photondimerror);
191 }
192
193 // p_vec[m_i1] must not be 0
194 const double elim = (m_clusterPars[m_i1] - x_vertex[m_i1]) / p_vec[m_i1];
195 const double mom = p_vec.norm();
196
197 // r'
198 Eigen::Matrix<double, 3, 1> residual3 = Eigen::Matrix<double, 3, 1>::Zero(3, 1);
199 residual3(0) = m_clusterPars[m_i2] - x_vertex[m_i2] - p_vec[m_i2] * elim;
200 residual3(1) = m_clusterPars[m_i3] - x_vertex[m_i3] - p_vec[m_i3] * elim;
201 residual3(2) = m_momentumScalingFactor * m_clusterPars[3] - mom; // scale measured energy by scaling factor
202
203 // dr'/dm | m:={xc,yc,zc,Ec} the measured quantities
204 Eigen::Matrix<double, 3, 4> P = Eigen::Matrix<double, 3, 4>::Zero(3, 4);
205 // deriving by the cluster pars
206 P(0, m_i2) = 1;
207 P(0, m_i1) = -p_vec[m_i2] / p_vec[m_i1];
208
209 P(1, m_i3) = 1;
210 P(1, m_i1) = -p_vec[m_i3] / p_vec[m_i1];
211 P(2, 3) = 1; // dE/dEc
212
213 p.getResiduals().segment(0, 3) = residual3;
214
215 p.getV() = P * m_covariance.selfadjointView<Eigen::Lower>() * P.transpose();
216
217 // dr'/dm | m:={x,y,z,px,py,pz,E}
218 // x := x_vertex (decay vertex of mother)
219 p.getH()(0, posindex + m_i1) = p_vec[m_i2] / p_vec[m_i1];
220 p.getH()(0, posindex + m_i2) = -1.0;
221 p.getH()(0, posindex + m_i3) = 0;
222
223 p.getH()(1, posindex + m_i1) = p_vec[m_i3] / p_vec[m_i1];
224 p.getH()(1, posindex + m_i2) = 0;
225 p.getH()(1, posindex + m_i3) = -1.0;
226
227 // elim already divided by p_vec[m_i1]
228 p.getH()(0, momindex + m_i1) = p_vec[m_i2] * elim / p_vec[m_i1];
229 p.getH()(0, momindex + m_i2) = -1. * elim;
230 p.getH()(0, momindex + m_i3) = 0;
231
232 p.getH()(1, momindex + m_i1) = p_vec[m_i3] * elim / p_vec[m_i1];
233 p.getH()(1, momindex + m_i2) = 0;
234 p.getH()(1, momindex + m_i3) = -1. * elim;
235
236 p.getH()(2, momindex + m_i1) = -1. * p_vec[m_i1] / mom;
237 p.getH()(2, momindex + m_i2) = -1. * p_vec[m_i2] / mom;
238 p.getH()(2, momindex + m_i3) = -1. * p_vec[m_i3] / mom;
239 // the photon does not store an energy in the state vector
240 // so no p.getH()(2, momindex + 3) here
241
242 return ErrCode(ErrCode::Status::success);
243 }
244
245}
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
Class to provide momentum-related information from ECLClusters.
Definition: ClusterUtils.h:36
const TMatrixDSym GetCovarianceMatrix3x3FromCluster(const ECLCluster *cluster)
Returns 3x3 covariance matrix (E, theta, phi)
The ParticleType class for identifying different particle types.
Definition: Const.h:408
static const ParticleType pi0
neutral pion particle
Definition: Const.h:674
static const ParticleType photon
photon particle
Definition: Const.h:673
ECL cluster data.
Definition: ECLCluster.h:27
EHypothesisBit
The hypothesis bits for this ECLCluster (Connected region (CR) is split using this hypothesis.
Definition: ECLCluster.h:31
Class to store reconstructed particles.
Definition: Particle.h:75
const ECLCluster * getECLCluster() const
Returns the pointer to the ECLCluster object that was used to create this Particle (if ParticleType =...
Definition: Particle.cc:891
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:454
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:560
ECLCluster::EHypothesisBit getECLClusterEHypothesisBit() const
Returns the ECLCluster EHypothesisBit for this Particle.
Definition: Particle.h:1001
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
Eigen::Matrix< double, -1, -1, 0, MAX_MATRIX_SIZE, MAX_MATRIX_SIZE > & getCovariance()
getter for the states covariance
Definition: FitParams.h:53
base class for all particles
Definition: ParticleBase.h:25
Belle2::Particle * particle() const
get basf2 particle
Definition: ParticleBase.h:92
virtual int posIndex() const
get vertex index (in statevector!)
Definition: ParticleBase.h:122
const ParticleBase * mother() const
getMother() / hasMother()
Definition: ParticleBase.h:98
class to store the projected residuals and the corresponding jacobian as well as the covariance matri...
Definition: Projection.h:18
base for RecoPhoton RecoTrack
Definition: RecoParticle.h:16
virtual int momIndex() const override
get momentum index
Definition: RecoParticle.h:42
Eigen::Matrix< double, 4, 4 > m_covariance
covariance (x_c,y_c,z_c,E_c) of measured pars
Definition: RecoPhoton.h:76
ErrCode initCovariance(FitParams &fitparams) const override
init covariance
Definition: RecoPhoton.cc:75
RecoPhoton(Belle2::Particle *bc, const ParticleBase *mother)
constructor
Definition: RecoPhoton.cc:26
ErrCode initParams()
update or init params
Definition: RecoPhoton.cc:92
const float m_momentumScalingFactor
scale the momentum / energy by this correction factor
Definition: RecoPhoton.h:86
int m_i3
another random index
Definition: RecoPhoton.h:83
int m_i1
index with the highest momentum.
Definition: RecoPhoton.h:79
virtual ErrCode initParticleWithMother(FitParams &fitparams) override
init particle with mother
Definition: RecoPhoton.cc:37
static bool useEnergy(const Belle2::Particle &cand)
has energy in fit params?
Definition: RecoPhoton.cc:63
int m_i2
random other index
Definition: RecoPhoton.h:81
virtual ErrCode initMotherlessParticle(FitParams &fitparams) override
init particle without mother
Definition: RecoPhoton.cc:58
bool m_init
was initialized*
Definition: RecoPhoton.h:67
ErrCode projectRecoConstraint(const FitParams &fitparams, Projection &p) const override
project photon constraint
Definition: RecoPhoton.cc:166
Eigen::Matrix< double, 1, 4 > m_clusterPars
constrains measured params (x_c, y_c, z_c, E_c)
Definition: RecoPhoton.h:73