Belle II Software  release-05-01-25
ClusterUtils.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2017-2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Torben Ferber (torben.ferber@desy.de) *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <analysis/ClusterUtility/ClusterUtils.h>
12 #include <framework/logging/Logger.h>
13 
14 using namespace Belle2;
15 
16 // -----------------------------------------------------------------------------
17 ClusterUtils::ClusterUtils() = default;
18 
19 // -----------------------------------------------------------------------------
21 {
22  // Use the geometry origin (0,0,0) and *not* the IP position
23  return Get4MomentumFromCluster(cluster, TVector3(0.0, 0.0, 0.0), hypo);
24 }
25 
27 {
28 
29  // Use the default vertex from the beam parameters if none is given.
30  return Get4MomentumFromCluster(cluster, GetIPPosition(), hypo);
31 }
32 
33 const TLorentzVector ClusterUtils::Get4MomentumFromCluster(const ECLCluster* cluster, const TVector3& vertex,
35 {
36 
37  // Get particle direction from vertex and reconstructed cluster position.
38  TVector3 direction = cluster->getClusterPosition() - vertex;
39 
40  // Always ignore mass here (even for neutral hadrons) therefore the magnitude
41  // of the momentum is equal to the cluster energy under this hypo.
42  const double E = cluster->getEnergy(hypo); //must not be changed or clusterE getters will be wrong
43  const double px = E * sin(direction.Theta()) * cos(direction.Phi());
44  const double py = E * sin(direction.Theta()) * sin(direction.Phi());
45  const double pz = E * cos(direction.Theta());
46 
47  const TLorentzVector l(px, py, pz, E);
48  return l;
49 }
50 
51 // -----------------------------------------------------------------------------
52 
54 {
55 
57 }
58 
59 const TMatrixDSym ClusterUtils::GetCovarianceMatrix4x4FromCluster(const ECLCluster* cluster, const TVector3& vertex,
60  const TMatrixDSym& covmatvertex, ECLCluster::EHypothesisBit hypo)
61 {
62 
63  // Get the covariance matrix (theta, phi, energy) from the ECL cluster.
64  TMatrixDSym covmatecl = cluster->getCovarianceMatrix3x3();
65 
66  // Combine into the total covariance matrix from the ECL cluster (0..2) and the IP(3..5) as two 3x3 block matrices.
67  TMatrixDSym covmatcombined(6);
68  for (int i = 0; i < 3; i++) {
69  for (int j = 0; j <= i ; j++) {
70  covmatcombined(i, j) = covmatcombined(j, i) = covmatecl(i, j);
71  covmatcombined(i + 3, j + 3) = covmatcombined(j + 3, i + 3) = covmatvertex(i, j);
72  }
73  }
74 
75  // Calculate the Jacobi matrix J = dpx/dE dpx/dphi dpx/dtheta dpx/dx dpy/dy dpz/dz
76  // ...
77  // dpx/dE dE/dphi dE/dtheta dpx/dx dE/dy dE/dz
78  TMatrixD jacobian(4, 6);
79 
80  const double R = cluster->getR();
81  const double energy = cluster->getEnergy(hypo);
82  const double theta = cluster->getTheta();
83  const double phi = cluster->getPhi();
84 
85  const double st = sin(theta);
86  const double ct = cos(theta);
87  const double sp = sin(phi);
88  const double cp = cos(phi);
89 
90  const double clx = R * st * cp; // cluster position x
91  const double cly = R * st * sp; // cluster position y
92  const double clz = R * ct; // cluster position z
93 
94  const double vx = vertex.X(); // vertex position x
95  const double vy = vertex.Y(); // vertex position y
96  const double vz = vertex.Z(); // vertex position z
97 
98  const double dx = clx - vx; // direction vector x
99  const double dy = cly - vy; // direction vector y
100  const double dz = clz - vz; // direction vector z
101  const double dx2 = dx * dx; // direction vector x squared
102  const double dy2 = dy * dy; // direction vector y squared
103  const double dz2 = dz * dz; // direction vector z squared
104 
105  const double r2 = (dx * dx + dy * dy + dz * dz);
106  const double r12 = sqrt(r2);
107  const double r32 = pow(r2, 1.5);
108 
109  // px = E * dx/r
110  jacobian(0, 0) = dx / r12; // dpx/dE
111  jacobian(0, 1) = -energy * R * ((dx * dy * cp) + ((dy2 + dz2) * sp)) * st / (r32); // dpx/dphi
112  jacobian(0, 2) = energy * R * (((dy2 + dz2) * cp * ct) - (dx * dy * ct * sp) + (dx * dz * st)) / (r32); // dpx/dtheta
113  jacobian(0, 3) = -energy * (dy2 + dz2) / (r32); // dpx/dvx
114  jacobian(0, 4) = energy * dx * dy / (r32); // dpx/dvy
115  jacobian(0, 5) = energy * dx * dz / (r32); //dpx/dvz
116 
117  // py = E * dy/r
118  jacobian(1, 0) = dy / r12; // dpy/dE
119  jacobian(1, 1) = energy * R * (((dx2 + dz2) * cp) + (dx * dy * sp)) * st / (r32); // dpy/dphi
120  jacobian(1, 2) = energy * R * ((-dx * dy * cp * ct) + ((dx2 + dz2) * ct * sp) + (dy * dz * st)) / (r32); // dpz/dtheta
121  jacobian(1, 3) = energy * dx * dy / (r32); // dpy/dvx
122  jacobian(1, 4) = -energy * (dx2 + dz2) / (r32); // dpy/dvy
123  jacobian(1, 5) = energy * dy * dz / (r32); //dpy/dvz
124 
125  // pz = E * dz/r
126  jacobian(2, 0) = dz / r12; // dpz/dE
127  jacobian(2, 1) = energy * R * dz * (-dy * cp + dx * sp) * st / (r32); // dpz/dphi
128  jacobian(2, 2) = -energy * R * ((dx * dz * cp * ct) + (dy * dz * ct * sp) + (dx2 + dy2) * st) / (r32); // dpz/dtheta
129  jacobian(2, 3) = energy * dx * dz / (r32); // dpz/dvx
130  jacobian(2, 4) = energy * dy * dz / (r32); // dpz/dvy
131  jacobian(2, 5) = -energy * (dx2 + dy2) / (r32); //dpz/ydvz
132 
133  // E
134  jacobian(3, 0) = 1.0; // dE/dE
135  jacobian(3, 1) = 0.0; // dE/dphi
136  jacobian(3, 2) = 0.0; // dE/dtheta
137  jacobian(3, 3) = 0.0; // dE/dvx
138  jacobian(3, 4) = 0.0; // dE/dvy
139  jacobian(3, 5) = 0.0; // dE/dvz
140 
141  TMatrixDSym covmatCart(4);
142  covmatCart = covmatcombined.Similarity(jacobian);
143 
144  return covmatCart;
145 }
146 
147 // -----------------------------------------------------------------------------
148 
150 {
151 
153 }
154 
155 const TMatrixDSym ClusterUtils::GetCovarianceMatrix7x7FromCluster(const ECLCluster* cluster, const TVector3& vertex,
156  const TMatrixDSym& covmatvertex, ECLCluster::EHypothesisBit hypo)
157 {
158 
159  TMatrixDSym covmat4x4 = GetCovarianceMatrix4x4FromCluster(cluster, vertex, covmatvertex, hypo);
160 
161  TMatrixDSym covmatCart(7);
162 
163  // px, py, pz, E
164  for (int i = 0; i < 4; i++) {
165  for (int j = 0; j <= i ; j++) {
166  covmatCart(i, j) = covmatCart(j, i) = covmat4x4(i, j);
167  }
168  }
169 
170  // x, y, z
171  for (int i = 0; i < 3; i++) {
172  for (int j = 0; j <= i ; j++) {
173  covmatCart(i + 4, j + 4) = covmatCart(j + 4, i + 4) = covmatvertex(i, j);
174  }
175  }
176 
177  return covmatCart;
178 }
179 
180 
181 // -----------------------------------------------------------------------------
183 {
184  if (!m_beamSpotDB) {
185  B2WARNING("Beamspot not available, using (0, 0, 0) as IP position instead.");
186  return TVector3(0.0, 0.0, 0.0);
187  } else return m_beamSpotDB->getIPPosition();
188 }
189 
190 // -----------------------------------------------------------------------------
192 {
193  if (!m_beamSpotDB) {
194  B2WARNING("Beam parameters not available, using ((1, 0, 0), (0, 1, 0), (0, 0, 1)) as IP covariance matrix instead.");
195 
196  TMatrixDSym covmat(3);
197  for (int i = 0; i < 3; ++i) {
198  covmat(i, i) = 1.0; // 1.0*1.0 cm^2
199  }
200  return covmat;
201  } else return m_beamSpotDB->getCovVertex();
202 }
Belle2::ECLCluster
ECL cluster data.
Definition: ECLCluster.h:39
Belle2::ECLCluster::EHypothesisBit
EHypothesisBit
The hypothesis bits for this ECLCluster (Connected region (CR) is split using this hypothesis.
Definition: ECLCluster.h:43
Belle2::ClusterUtils::Get4MomentumFromCluster
const TLorentzVector Get4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns four momentum vector.
Definition: ClusterUtils.cc:26
Belle2::ClusterUtils::GetCovarianceMatrix4x4FromCluster
const TMatrixDSym GetCovarianceMatrix4x4FromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns 4x4 covariance matrix (px, py, pz, E)
Definition: ClusterUtils.cc:53
Belle2::ClusterUtils::GetIPPosition
const TVector3 GetIPPosition()
Returns default IP position from beam parameters.
Definition: ClusterUtils.cc:182
Belle2::ClusterUtils::GetCovarianceMatrix7x7FromCluster
const TMatrixDSym GetCovarianceMatrix7x7FromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns 7x7 covariance matrix (px, py, pz, E, x, y, z)
Definition: ClusterUtils.cc:149
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ClusterUtils::GetCluster4MomentumFromCluster
const TLorentzVector GetCluster4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns cluster four momentum vector.
Definition: ClusterUtils.cc:20
Belle2::ClusterUtils::m_beamSpotDB
DBObjPtr< BeamSpot > m_beamSpotDB
Beam spot database object.
Definition: ClusterUtils.h:111
Belle2::ClusterUtils::ClusterUtils
ClusterUtils()
Constructor.
Belle2::ClusterUtils::GetIPPositionCovarianceMatrix
const TMatrixDSym GetIPPositionCovarianceMatrix()
Returns default IP position covariance matrix from beam parameters.
Definition: ClusterUtils.cc:191