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