Belle II Software development
VertexVariables.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/variables/VertexVariables.h>
10
11#include <analysis/dataobjects/Particle.h>
12#include <analysis/utility/ReferenceFrame.h>
13
14#include <framework/database/DBObjPtr.h>
15#include <framework/logging/Logger.h>
16#include <framework/utilities/Conversion.h>
17
18#include <TMatrixFSym.h>
19
20#include <mdst/dbobjects/BeamSpot.h>
21#include <mdst/dataobjects/MCParticle.h>
22#include <mdst/dataobjects/Track.h>
23#include <mdst/dataobjects/TrackFitResult.h>
24
25
26namespace Belle2 {
31 class Particle;
32
33 namespace Variable {
34
35 // Generate Vertex information
36
37 double mcDecayVertexX(const Particle* part)
38 {
39 auto* mcparticle = part->getMCParticle();
40 if (!mcparticle) return Const::doubleNaN;
41 return mcparticle->getDecayVertex().X();
42 }
43
44 double mcDecayVertexY(const Particle* part)
45 {
46 auto* mcparticle = part->getMCParticle();
47 if (!mcparticle) return Const::doubleNaN;
48 return mcparticle->getDecayVertex().Y();
49 }
50
51 double mcDecayVertexZ(const Particle* part)
52 {
53 auto* mcparticle = part->getMCParticle();
54 if (!mcparticle) return Const::doubleNaN;
55 return mcparticle->getDecayVertex().Z();
56 }
57
58 double mcDecayVertexRho(const Particle* part)
59 {
60 auto* mcparticle = part->getMCParticle();
61 if (!mcparticle) return Const::doubleNaN;
62 return mcparticle->getDecayVertex().Rho();
63 }
64
65 B2Vector3D getMcDecayVertexFromIP(const MCParticle* mcparticle)
66 {
67 static DBObjPtr<BeamSpot> beamSpotDB;
68 const auto& frame = ReferenceFrame::GetCurrent();
69 return frame.getVertex(mcparticle->getDecayVertex() - ROOT::Math::XYZVector(beamSpotDB->getIPPosition()));
70 }
71
72 double mcDecayVertexFromIPX(const Particle* part)
73 {
74 auto* mcparticle = part->getMCParticle();
75 if (!mcparticle) return Const::doubleNaN;
76 return getMcDecayVertexFromIP(mcparticle).X();
77 }
78
79 double mcDecayVertexFromIPY(const Particle* part)
80 {
81 auto* mcparticle = part->getMCParticle();
82 if (!mcparticle) return Const::doubleNaN;
83 return getMcDecayVertexFromIP(mcparticle).Y();
84 }
85
86 double mcDecayVertexFromIPZ(const Particle* part)
87 {
88 auto* mcparticle = part->getMCParticle();
89 if (!mcparticle) return Const::doubleNaN;
90 return getMcDecayVertexFromIP(mcparticle).Z();
91 }
92
93 double mcDecayVertexFromIPRho(const Particle* part)
94 {
95 auto* mcparticle = part->getMCParticle();
96 if (!mcparticle) return Const::doubleNaN;
97 return getMcDecayVertexFromIP(mcparticle).Perp();
98 }
99
100 double mcDecayVertexFromIPDistance(const Particle* part)
101 {
102 auto* mcparticle = part->getMCParticle();
103 if (!mcparticle) return Const::doubleNaN;
104 return getMcDecayVertexFromIP(mcparticle).Mag();
105 }
106
107 double mcProductionVertexX(const Particle* part)
108 {
109 auto* mcparticle = part->getMCParticle();
110 if (!mcparticle) return Const::doubleNaN;
111 return mcparticle->getProductionVertex().X();
112 }
113
114 double mcProductionVertexY(const Particle* part)
115 {
116 auto* mcparticle = part->getMCParticle();
117 if (!mcparticle) return Const::doubleNaN;
118 return mcparticle->getProductionVertex().Y();
119 }
120
121 double mcProductionVertexZ(const Particle* part)
122 {
123 auto* mcparticle = part->getMCParticle();
124 if (!mcparticle) return Const::doubleNaN;
125 return mcparticle->getProductionVertex().Z();
126 }
127
128 B2Vector3D getMcProductionVertexFromIP(const MCParticle* mcparticle)
129 {
130 static DBObjPtr<BeamSpot> beamSpotDB;
131 const auto& frame = ReferenceFrame::GetCurrent();
132 return frame.getVertex(mcparticle->getProductionVertex() - ROOT::Math::XYZVector(beamSpotDB->getIPPosition()));
133 }
134
135 double mcProductionVertexFromIPX(const Particle* part)
136 {
137 auto* mcparticle = part->getMCParticle();
138 if (!mcparticle) return Const::doubleNaN;
139 return getMcProductionVertexFromIP(mcparticle).X();
140 }
141
142 double mcProductionVertexFromIPY(const Particle* part)
143 {
144 auto* mcparticle = part->getMCParticle();
145 if (!mcparticle) return Const::doubleNaN;
146 return getMcProductionVertexFromIP(mcparticle).Y();
147 }
148
149 double mcProductionVertexFromIPZ(const Particle* part)
150 {
151 auto* mcparticle = part->getMCParticle();
152 if (!mcparticle) return Const::doubleNaN;
153 return getMcProductionVertexFromIP(mcparticle).Z();
154 }
155
156 // vertex or POCA with respect to origin ------------------------------
157
158 double particleX(const Particle* part)
159 {
160 const auto& frame = ReferenceFrame::GetCurrent();
161 return frame.getVertex(part).X();
162 }
163
164 double particleY(const Particle* part)
165 {
166 const auto& frame = ReferenceFrame::GetCurrent();
167 return frame.getVertex(part).Y();
168 }
169
170 double particleZ(const Particle* part)
171 {
172 const auto& frame = ReferenceFrame::GetCurrent();
173 return frame.getVertex(part).Z();
174 }
175
176 inline double getParticleUncertaintyByIndex(const Particle* part, unsigned int index)
177 {
178 if (!part) {
179 B2FATAL("The particle provide does not exist.");
180 }
181 const auto& errMatrix = part->getVertexErrorMatrix();
182 return std::sqrt(errMatrix(index, index));
183 }
184
185 double particleDXUncertainty(const Particle* part)
186 {
187 // uncertainty on x (with respect to the origin)
188 return getParticleUncertaintyByIndex(part, 0);
189 }
190
191 double particleDYUncertainty(const Particle* part)
192 {
193 // uncertainty on y (with respect to the origin)
194 return getParticleUncertaintyByIndex(part, 1);
195 }
196
197 double particleDZUncertainty(const Particle* part)
198 {
199 // uncertainty on z (with respect to the origin)
200 return getParticleUncertaintyByIndex(part, 2);
201 }
202
203 //----------------------------------------------------------------------------------
204 // vertex or POCA with respect to measured IP
205
206 B2Vector3D getVertexD(const Particle* part)
207 {
208 static DBObjPtr<BeamSpot> beamSpotDB;
209 const auto& frame = ReferenceFrame::GetCurrent();
210 auto trackFit = part->getTrackFitResult();
211 if (!trackFit)
212 return frame.getVertex(part->getVertex() - ROOT::Math::XYZVector(beamSpotDB->getIPPosition()));
213
214 UncertainHelix helix = trackFit->getUncertainHelix();
215 helix.passiveMoveBy(ROOT::Math::XYZVector(beamSpotDB->getIPPosition()));
216 return frame.getVertex(helix.getPerigee());
217 }
218
219
220 double particleDX(const Particle* part)
221 {
222 return getVertexD(part).X();
223 }
224
225 double particleDY(const Particle* part)
226 {
227 return getVertexD(part).Y();
228 }
229
230 double particleDZ(const Particle* part)
231 {
232 return getVertexD(part).Z();
233 }
234
235 double particleDRho(const Particle* part)
236 {
237 return getVertexD(part).Perp();
238 }
239
240 double particleDPhi(const Particle* part)
241 {
242 return getVertexD(part).Phi();
243 }
244
245 double particleDCosTheta(const Particle* part)
246 {
247 return getVertexD(part).CosTheta();
248 }
249
250 double particleDistance(const Particle* part)
251 {
252 return getVertexD(part).Mag();
253 }
254
255 double particleDistanceSignificance(const Particle* part)
256 {
257 // significance is defined as s = r/sigma_r, therefore:
258 // s &= \frac{r}{\sqrt{ \sum_{ij} \frac{\partial r}{x_i} V_{ij} \frac{\partial r}{x_j}}}
259 // &= \frac{r^2}{\sqrt{\vec{x}V\vec{x}}}
260 // where:
261 // r &= \sqrt{\vec{x}*\vec{x}}
262 // and V_{ij} is the covariance matrix
263 static DBObjPtr<BeamSpot> beamSpotDB;
264 const auto& frame = ReferenceFrame::GetCurrent();
265 const B2Vector3D& vertex = frame.getVertex(part->getVertex() - ROOT::Math::XYZVector(beamSpotDB->getIPPosition()));
266 const TMatrixFSym& vertexErr = frame.getVertexErrorMatrix(static_cast<TMatrixDSym>(part->getVertexErrorMatrix()) +
267 beamSpotDB->getCovVertex());
268 const double denominator = vertex * B2Vector3D(vertexErr * vertex);
269 if (denominator <= 0) return Const::doubleNaN;
270
271 return vertex.Mag2() / std::sqrt(denominator);
272 }
273
274 // Production vertex position
275
276 double particleProductionX(const Particle* part)
277 {
278 if (!part->hasExtraInfo("prodVertX")) return Const::doubleNaN;
279 return part->getExtraInfo("prodVertX");
280 }
281
282 double particleProductionY(const Particle* part)
283 {
284 if (!part->hasExtraInfo("prodVertY")) return Const::doubleNaN;
285 return part->getExtraInfo("prodVertY");
286 }
287
288 double particleProductionZ(const Particle* part)
289 {
290 if (!part->hasExtraInfo("prodVertZ")) return Const::doubleNaN;
291 return part->getExtraInfo("prodVertZ");
292 }
293
294 // Production vertex covariance matrix
295 double particleProductionCovElement(const Particle* part, const std::vector<double>& indices)
296 {
297 if (indices.size() != 2) {
298 B2FATAL("Number of arguments of prodVertexCov function is incorrect!");
299 }
300
301 int ielement = std::lround(indices[0]);
302 int jelement = std::lround(indices[1]);
303
304 if (std::min(ielement, jelement) < 0 || std::max(ielement, jelement) > 2) {
305 B2ERROR("Range of indexes of prodVertexCov function is incorrect!");
306 }
307
308 const std::vector<char> names = {'x', 'y', 'z'};
309 const std::string prodVertS = Form("prodVertS%c%c", names[ielement], names[jelement]);
310
311 if (!part->hasExtraInfo(prodVertS)) return Const::doubleNaN;
312 return part->getExtraInfo(prodVertS);
313 }
314
315 double particleProductionXErr(const Particle* part)
316 {
317 if (!part->hasExtraInfo("prodVertSxx")) return Const::doubleNaN;
318 return std::sqrt(part->getExtraInfo("prodVertSxx"));
319 }
320
321 double particleProductionYErr(const Particle* part)
322 {
323 if (!part->hasExtraInfo("prodVertSyy")) return Const::doubleNaN;
324 return std::sqrt(part->getExtraInfo("prodVertSyy"));
325 }
326
327 double particleProductionZErr(const Particle* part)
328 {
329 if (!part->hasExtraInfo("prodVertSzz")) return Const::doubleNaN;
330 return std::sqrt(part->getExtraInfo("prodVertSzz"));
331 }
332
333 VARIABLE_GROUP("Vertex Information");
334 // Generated vertex information
335 REGISTER_VARIABLE("mcDecayVertexX", mcDecayVertexX,
336 "Returns the x position of the decay vertex of the matched generated particle. Returns nan if the particle has no matched generated particle.\n\n",
337 "cm");
338 REGISTER_VARIABLE("mcDecayVertexY", mcDecayVertexY,
339 "Returns the y position of the decay vertex of the matched generated particle. Returns nan if the particle has no matched generated particle.\n\n",
340 "cm");
341 REGISTER_VARIABLE("mcDecayVertexZ", mcDecayVertexZ,
342 "Returns the z position of the decay vertex of the matched generated particle. Returns nan if the particle has no matched generated particle.\n\n",
343 "cm");
344 REGISTER_VARIABLE("mcDecayVertexRho", mcDecayVertexRho,
345 "Returns the transverse position of the decay vertex of the matched generated particle. Returns nan if the particle has no matched generated particle.\n\n",
346 "cm");
347 REGISTER_VARIABLE("mcDecayVertexFromIPX", mcDecayVertexFromIPX,
348 "Returns the x position of the decay vertex of the matched generated particle wrt the IP. Returns nan if the particle has no matched generated particle.\n\n",
349 "cm");
350 REGISTER_VARIABLE("mcDecayVertexFromIPY", mcDecayVertexFromIPY,
351 "Returns the y position of the decay vertex of the matched generated particle wrt the IP. Returns nan if the particle has no matched generated particle.\n\n",
352 "cm");
353 REGISTER_VARIABLE("mcDecayVertexFromIPZ", mcDecayVertexFromIPZ,
354 "Returns the z position of the decay vertex of the matched generated particle wrt the IP. Returns nan if the particle has no matched generated particle.\n\n",
355 "cm");
356 REGISTER_VARIABLE("mcDecayVertexFromIPRho", mcDecayVertexFromIPRho,
357 "Returns the transverse position of the decay vertex of the matched generated particle wrt the IP. Returns nan if the particle has no matched generated particle.\n\n",
358 "cm");
359 REGISTER_VARIABLE("mcDecayVertexFromIPDistance", mcDecayVertexFromIPDistance,
360 "Returns the distance of the decay vertex of the matched generated particle from the IP. Returns nan if the particle has no matched generated particle.\n\n",
361 "cm");
362 REGISTER_VARIABLE("mcProductionVertexX", mcProductionVertexX,
363 "Returns the x position of production vertex of matched generated particle. Returns nan if the particle has no matched generated particle.\n\n",
364 "cm");
365 REGISTER_VARIABLE("mcProductionVertexY", mcProductionVertexY,
366 "Returns the y position of production vertex of matched generated particle. Returns nan if the particle has no matched generated particle.\n\n",
367 "cm");
368 REGISTER_VARIABLE("mcProductionVertexZ", mcProductionVertexZ,
369 "Returns the z position of production vertex of matched generated particle. Returns nan if the particle has no matched generated particle.\n\n",
370 "cm");
371 REGISTER_VARIABLE("mcProductionVertexFromIPX", mcProductionVertexFromIPX,
372 "Returns the x position of the production vertex of the matched generated particle wrt the IP. Returns nan if the particle has no matched generated particle.\n\n",
373 "cm");
374 REGISTER_VARIABLE("mcProductionVertexFromIPY", mcProductionVertexFromIPY,
375 "Returns the y position of the production vertex of the matched generated particle wrt the IP. Returns nan if the particle has no matched generated particle.\n\n",
376 "cm");
377 REGISTER_VARIABLE("mcProductionVertexFromIPZ", mcProductionVertexFromIPZ,
378 "Returns the z position of the production vertex of the matched generated particle wrt the IP. Returns nan if the particle has no matched generated particle.\n\n",
379 "cm");
380
381 // Decay vertex position
382 REGISTER_VARIABLE("distance", particleDistance,
383 R"DOC(3D distance between the IP and the particle decay vertex, if available.
384
385In case the particle has been created from a track, the distance is defined between the POCA and IP.
386If the particle is built from an ECL cluster, the decay vertex is set to the nominal IP.
387If the particle is created from a KLM cluster, the distance is calculated between the IP and the cluster itself.
388
389)DOC", "cm");
390
391 REGISTER_VARIABLE("significanceOfDistance", particleDistanceSignificance,
392 "significance of distance from vertex or POCA to interaction point(-1 in case of numerical problems)");
393 REGISTER_VARIABLE("dx", particleDX, "vertex or POCA in case of tracks x with respect to IP\n\n", "cm");
394 REGISTER_VARIABLE("dy", particleDY, "vertex or POCA in case of tracks y with respect to IP\n\n", "cm");
395 REGISTER_VARIABLE("dz", particleDZ, "vertex or POCA in case of tracks z with respect to IP\n\n", "cm");
396 REGISTER_VARIABLE("x", particleX,
397 "x coordinate of vertex in case of composite particle, or point of closest approach (POCA) in case of a track\n\n", "cm");
398 REGISTER_VARIABLE("y", particleY,
399 "y coordinate of vertex in case of composite particle, or point of closest approach (POCA) in case of a track\n\n", "cm");
400 REGISTER_VARIABLE("z", particleZ,
401 "z coordinate of vertex in case of composite particle, or point of closest approach (POCA) in case of a track\n\n", "cm");
402 REGISTER_VARIABLE("x_uncertainty", particleDXUncertainty, "uncertainty on x (measured with respect to the origin)\n\n", "cm");
403 REGISTER_VARIABLE("y_uncertainty", particleDYUncertainty, "uncertainty on y (measured with respect to the origin)\n\n", "cm");
404 REGISTER_VARIABLE("z_uncertainty", particleDZUncertainty, "uncertainty on z (measured with respect to the origin)\n\n", "cm");
405 REGISTER_VARIABLE("dr", particleDRho, "transverse distance with respect to IP for a vertex; track abs(d0) relative to IP for a track.\n\n",
406 "cm");
407 REGISTER_VARIABLE("dphi", particleDPhi, "vertex azimuthal angle of the vertex or POCA in degrees with respect to IP\n\n", "rad");
408 REGISTER_VARIABLE("dcosTheta", particleDCosTheta, "vertex or POCA polar angle with respect to IP");
409 // Production vertex position
410 REGISTER_VARIABLE("prodVertexX", particleProductionX,
411 "Returns the x position of particle production vertex. Returns NaN if particle has no production vertex.\n\n", "cm");
412 REGISTER_VARIABLE("prodVertexY", particleProductionY,
413 "Returns the y position of particle production vertex.\n\n", "cm");
414 REGISTER_VARIABLE("prodVertexZ", particleProductionZ,
415 "Returns the z position of particle production vertex.\n\n", "cm");
416 // Production vertex covariance matrix
417 REGISTER_VARIABLE("prodVertexCov(i,j)", particleProductionCovElement,
418 "Returns the ij covariance matrix component of particle production vertex, arguments i,j should be 0, 1 or 2. Returns NaN if particle has no production covariance matrix.\n\n",
419 ":math:`\\text{cm}^2`");
420 REGISTER_VARIABLE("prodVertexXErr", particleProductionXErr,
421 "Returns the x position uncertainty of particle production vertex. Returns NaN if particle has no production vertex.\n\n", "cm");
422 REGISTER_VARIABLE("prodVertexYErr", particleProductionYErr,
423 "Returns the y position uncertainty of particle production vertex.\n\n", "cm");
424 REGISTER_VARIABLE("prodVertexZErr", particleProductionZErr,
425 "Returns the z position uncertainty of particle production vertex.\n\n", "cm");
426
427 }
429}
DataType Phi() const
The azimuth angle.
Definition: B2Vector3.h:151
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:435
DataType CosTheta() const
Cosine of the polar angle.
Definition: B2Vector3.h:155
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
DataType Mag() const
The magnitude (rho in spherical coordinate system).
Definition: B2Vector3.h:159
DataType Perp() const
The transverse component (R in cylindrical coordinate system).
Definition: B2Vector3.h:200
static const double doubleNaN
quiet_NaN
Definition: Const.h:703
static const ReferenceFrame & GetCurrent()
Get current rest frame.
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
Abstract base class for different kinds of events.