Belle II Software light-2406-ragdoll
ClusterUtils.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <analysis/ClusterUtility/ClusterUtils.h>
10#include <framework/logging/Logger.h>
11
12using namespace Belle2;
13using namespace ROOT::Math;
14
15// -----------------------------------------------------------------------------
17
18// -----------------------------------------------------------------------------
20{
21 // Use the geometry origin (0,0,0) and *not* the IP position
22 return Get4MomentumFromCluster(cluster, XYZVector(0.0, 0.0, 0.0), hypo);
23}
24
26{
27
28 // Use the default vertex from the beam parameters if none is given.
29 return Get4MomentumFromCluster(cluster, GetIPPosition(), hypo);
30}
31
32const PxPyPzEVector ClusterUtils::Get4MomentumFromCluster(const ECLCluster* cluster, const ROOT::Math::XYZVector& vertex,
34{
35
36 // Get particle direction from vertex and reconstructed cluster position.
37 XYZVector direction(cluster->getClusterPosition().X() - vertex.X(),
38 cluster->getClusterPosition().Y() - vertex.Y(),
39 cluster->getClusterPosition().Z() - vertex.Z());
40
41 // Always ignore mass here (even for neutral hadrons) therefore the magnitude
42 // of the momentum is equal to the cluster energy under this hypo.
43 const double E = cluster->getEnergy(hypo); // must not be changed or clusterE getters will be wrong
44 const double px = E * sin(direction.Theta()) * cos(direction.Phi());
45 const double py = E * sin(direction.Theta()) * sin(direction.Phi());
46 const double pz = E * cos(direction.Theta());
47
48 const PxPyPzEVector l(px, py, pz, E);
49 return l;
50}
51
52// -----------------------------------------------------------------------------
53
55{
56
57 return GetJacobiMatrix4x6FromCluster(cluster, GetIPPosition(), hypo);
58}
59
60const TMatrixD ClusterUtils::GetJacobiMatrix4x6FromCluster(const ECLCluster* cluster, const ROOT::Math::XYZVector& vertex,
62{
63
64 // Calculate the Jacobi matrix J = dpx/dE dpx/dphi dpx/dtheta dpx/dx dpy/dy dpz/dz
65 // ...
66 // dpx/dE dE/dphi dE/dtheta dpx/dx dE/dy dE/dz
67 TMatrixD jacobian(4, 6);
68
69 const double R = cluster->getR();
70 const double energy = cluster->getEnergy(hypo);
71 const double theta = cluster->getTheta();
72 const double phi = cluster->getPhi();
73
74 const double st = sin(theta);
75 const double ct = cos(theta);
76 const double sp = sin(phi);
77 const double cp = cos(phi);
78
79 const double clx = R * st * cp; // cluster position x
80 const double cly = R * st * sp; // cluster position y
81 const double clz = R * ct; // cluster position z
82
83 const double vx = vertex.X(); // vertex position x
84 const double vy = vertex.Y(); // vertex position y
85 const double vz = vertex.Z(); // vertex position z
86
87 const double dx = clx - vx; // direction vector x
88 const double dy = cly - vy; // direction vector y
89 const double dz = clz - vz; // direction vector z
90 const double dx2 = dx * dx; // direction vector x squared
91 const double dy2 = dy * dy; // direction vector y squared
92 const double dz2 = dz * dz; // direction vector z squared
93
94 const double r2 = (dx * dx + dy * dy + dz * dz);
95 const double r12 = sqrt(r2);
96 const double r32 = pow(r2, 1.5);
97
98 // px = E * dx/r
99 jacobian(0, 0) = dx / r12; // dpx/dE
100 jacobian(0, 1) = -energy * R * ((dx * dy * cp) + ((dy2 + dz2) * sp)) * st / (r32); // dpx/dphi
101 jacobian(0, 2) = energy * R * (((dy2 + dz2) * cp * ct) - (dx * dy * ct * sp) + (dx * dz * st)) / (r32); // dpx/dtheta
102 jacobian(0, 3) = -energy * (dy2 + dz2) / (r32); // dpx/dvx
103 jacobian(0, 4) = energy * dx * dy / (r32); // dpx/dvy
104 jacobian(0, 5) = energy * dx * dz / (r32); // dpx/dvz
105
106 // py = E * dy/r
107 jacobian(1, 0) = dy / r12; // dpy/dE
108 jacobian(1, 1) = energy * R * (((dx2 + dz2) * cp) + (dx * dy * sp)) * st / (r32); // dpy/dphi
109 jacobian(1, 2) = energy * R * ((-dx * dy * cp * ct) + ((dx2 + dz2) * ct * sp) + (dy * dz * st)) / (r32); // dpz/dtheta
110 jacobian(1, 3) = energy * dx * dy / (r32); // dpy/dvx
111 jacobian(1, 4) = -energy * (dx2 + dz2) / (r32); // dpy/dvy
112 jacobian(1, 5) = energy * dy * dz / (r32); // dpy/dvz
113
114 // pz = E * dz/r
115 jacobian(2, 0) = dz / r12; // dpz/dE
116 jacobian(2, 1) = energy * R * dz * (-dy * cp + dx * sp) * st / (r32); // dpz/dphi
117 jacobian(2, 2) = -energy * R * ((dx * dz * cp * ct) + (dy * dz * ct * sp) + (dx2 + dy2) * st) / (r32); // dpz/dtheta
118 jacobian(2, 3) = energy * dx * dz / (r32); // dpz/dvx
119 jacobian(2, 4) = energy * dy * dz / (r32); // dpz/dvy
120 jacobian(2, 5) = -energy * (dx2 + dy2) / (r32); // dpz/ydvz
121
122 // E
123 jacobian(3, 0) = 1.0; // dE/dE
124 jacobian(3, 1) = 0.0; // dE/dphi
125 jacobian(3, 2) = 0.0; // dE/dtheta
126 jacobian(3, 3) = 0.0; // dE/dvx
127 jacobian(3, 4) = 0.0; // dE/dvy
128 jacobian(3, 5) = 0.0; // dE/dvz
129
130 // jacobian.Print();
131 return jacobian;
132}
133
135{
136 // Get the covariance matrix (energy, phi, theta) from the ECL cluster.
137 TMatrixDSym covmatecl = cluster->getCovarianceMatrix3x3();
138
142 if (m_photonEnergyResolutionDB.isValid() and cluster->hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
143 double pEnergy = cluster->getEnergy(ECLCluster::EHypothesisBit::c_nPhotons);
144 double pTheta = cluster->getTheta();
145 double pPhi = cluster->getPhi();
146
147 double energyCovarianceElement = m_photonEnergyResolutionDB->getRelativeEnergyResolution(pEnergy, pTheta, pPhi);
148 if (energyCovarianceElement != -1) {
149 covmatecl(0, 0) = energyCovarianceElement * pEnergy * energyCovarianceElement * pEnergy;
150 }
151 }
152
153 return covmatecl;
154}
155
156const TMatrixDSym ClusterUtils::GetCovarianceMatrix4x4FromCluster(const ECLCluster* cluster, const TMatrixD& jacobiMatrix)
157{
158
160}
161
162const TMatrixDSym ClusterUtils::GetCovarianceMatrix4x4FromCluster(const ECLCluster* cluster, const TMatrixDSym& covmatvertex,
163 const TMatrixD& jacobiMatrix)
164{
165
166 // Get the covariance matrix (energy, phi, theta)
167 TMatrixDSym covmatecl = GetCovarianceMatrix3x3FromCluster(cluster);
168
169 // Combine into the total covariance matrix from the ECL cluster (0..2) and the IP(3..5) as two 3x3 block matrices.
170 TMatrixDSym covmatcombined(6);
171 for (int i = 0; i < 3; i++) {
172 for (int j = 0; j <= i; j++) {
173 covmatcombined(i, j) = covmatcombined(j, i) = covmatecl(i, j);
174 covmatcombined(i + 3, j + 3) = covmatcombined(j + 3, i + 3) = covmatvertex(i, j);
175 }
176 }
177
178 TMatrixDSym covmatCart(4);
179 covmatCart = covmatcombined.Similarity(jacobiMatrix);
180 return covmatCart;
181}
182
183// -----------------------------------------------------------------------------
184
185const TMatrixDSym ClusterUtils::GetCovarianceMatrix7x7FromCluster(const ECLCluster* cluster, const TMatrixD& jacobiMatrix)
186{
187
189}
190
192 const TMatrixDSym& covmatvertex, const TMatrixD& jacobiMatrix)
193{
194 TMatrixDSym covmat4x4 = GetCovarianceMatrix4x4FromCluster(cluster, covmatvertex, jacobiMatrix);
195
196 TMatrixDSym covmatCart(7);
197
198 // px, py, pz, E
199 for (int i = 0; i < 4; i++) {
200 for (int j = 0; j <= i; j++) {
201 covmatCart(i, j) = covmatCart(j, i) = covmat4x4(i, j);
202 }
203 }
204
205 // x, y, z
206 for (int i = 0; i < 3; i++) {
207 for (int j = 0; j <= i; j++) {
208 covmatCart(i + 4, j + 4) = covmatCart(j + 4, i + 4) = covmatvertex(i, j);
209 }
210 }
211
212 return covmatCart;
213}
214
215// -----------------------------------------------------------------------------
217{
218 if (!m_beamSpotDB) {
219 B2WARNING("Beamspot not available, using (0, 0, 0) as IP position instead.");
220 return XYZVector(0.0, 0.0, 0.0);
221 } else
222 return XYZVector(m_beamSpotDB->getIPPosition().X(), m_beamSpotDB->getIPPosition().Y(), m_beamSpotDB->getIPPosition().Z());
223}
224
225// -----------------------------------------------------------------------------
227{
228 if (!m_beamSpotDB) {
229 B2WARNING("Beam parameters not available, using ((1, 0, 0), (0, 1, 0), (0, 0, 1)) as IP covariance matrix instead.");
230
231 TMatrixDSym covmat(3);
232 for (int i = 0; i < 3; ++i) {
233 covmat(i, i) = 1.0; // 1.0*1.0 cm^2
234 }
235 return covmat;
236 } else
237 return m_beamSpotDB->getCovVertex();
238}
const TMatrixDSym GetCovarianceMatrix7x7FromCluster(const ECLCluster *cluster, const TMatrixD &jacobiMatrix)
Returns 7x7 covariance matrix (px, py, pz, E, x, y, z)
const ROOT::Math::PxPyPzEVector Get4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns four momentum vector.
Definition: ClusterUtils.cc:25
OptionalDBObjPtr< ECLPhotonEnergyResolution > m_photonEnergyResolutionDB
Photon energy resolution database object.
Definition: ClusterUtils.h:130
const TMatrixD GetJacobiMatrix4x6FromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns 4x6 Jacobi matrix (px, py, pz, E)
Definition: ClusterUtils.cc:54
const ROOT::Math::PxPyPzEVector GetCluster4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns cluster four momentum vector.
Definition: ClusterUtils.cc:19
const TMatrixDSym GetCovarianceMatrix4x4FromCluster(const ECLCluster *cluster, const TMatrixD &jacobiMatrix)
Returns 4x4 covariance matrix (px, py, pz, E)
DBObjPtr< BeamSpot > m_beamSpotDB
Beam spot database object.
Definition: ClusterUtils.h:125
const TMatrixDSym GetCovarianceMatrix3x3FromCluster(const ECLCluster *cluster)
Returns 3x3 covariance matrix (E, theta, phi)
const TMatrixDSym GetIPPositionCovarianceMatrix()
Returns default IP position covariance matrix from beam parameters.
ClusterUtils()
Constructor.
const ROOT::Math::XYZVector GetIPPosition()
Returns default IP position from beam parameters.
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
@ c_nPhotons
CR is split into n photons (N1)
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24