Belle II Software development
FlightInfoVariables.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// Own header.
10#include <analysis/variables/FlightInfoVariables.h>
11
12// include VariableManager
13#include <analysis/VariableManager/Manager.h>
14
15#include <framework/database/DBObjPtr.h>
16#include <framework/logging/Logger.h>
17
18#include <analysis/dataobjects/Particle.h>
19#include <mdst/dataobjects/MCParticle.h>
20#include <mdst/dbobjects/BeamSpot.h>
21
22#include <TMatrixFSym.h>
23
24#include <boost/algorithm/string.hpp>
25
26namespace Belle2 {
32 namespace Variable {
33
34 //Function to check if a beam constrained rave fit has been performed
35 bool hasRAVEBeamConstrainedProductionVertex(const Particle* particle)
36 {
37 bool hasRAVEProdVertex = true;
38 std::vector<std::string> directions = {"x", "y", "z"};
39 for (auto ielement : directions) {
40 std::string prodVertPositionElement = boost::str(boost::format("prodVert%s") % boost::to_upper_copy(ielement));
41 hasRAVEProdVertex &= particle -> hasExtraInfo(prodVertPositionElement);
42 for (auto jelement : directions) {
43 std::string prodVertCovarianceElement = boost::str(boost::format("prodVertS%s%s") % ielement % jelement);
44 hasRAVEProdVertex &= particle -> hasExtraInfo(prodVertCovarianceElement);
45 }
46 }
47 return hasRAVEProdVertex;
48 }
49
50 // Helper function for flight distance and its uncertainty (provided as it is)
51 inline double getFlightInfoBtw(const Particle* particle, const Particle* daughter, double& outErr, const std::string& mode,
52 const bool motherToGranddaughter = false)
53 {
54 if (!particle || !daughter) {
55 outErr = Const::doubleNaN;
56 return Const::doubleNaN;
57 }
58 // check if the particle source is a composite particle or a V0
59 if ((particle->getParticleSource() != Particle::EParticleSourceObject::c_Composite and
60 particle->getParticleSource() != Particle::EParticleSourceObject::c_V0) or
61 (daughter->getParticleSource() != Particle::EParticleSourceObject::c_Composite and
62 daughter->getParticleSource() != Particle::EParticleSourceObject::c_V0)) {
63 B2WARNING("Attempting to calculate flight " << mode << " for neither composite particle nor V0");
64 outErr = Const::doubleNaN;
65 return Const::doubleNaN;
66 }
67 if (!(mode == "distance") && !(mode == "time")) {
68 B2WARNING("FlightInfo helper function called with mode '" << mode
69 << "'. Only 'distance' and 'time' are available.");
70 outErr = Const::doubleNaN;
71 return Const::doubleNaN;
72 }
73 // get TreeFitter values if they exist.
74 // Bypass this in case the variables are requested for the granddaughter with respect to the mother as
75 // TreeFitter will return the values of the granddaughter with respect to the daughter
76 if (!motherToGranddaughter) {
77 if (mode == "distance" &&
78 daughter->hasExtraInfo("decayLength") &&
79 daughter->hasExtraInfo("decayLengthErr")) {
80 outErr = daughter -> getExtraInfo("decayLengthErr");
81 return daughter -> getExtraInfo("decayLength");
82 }
83 if (mode == "time" &&
84 daughter->hasExtraInfo("lifeTime") &&
85 daughter->hasExtraInfo("lifeTimeErr")) {
86 outErr = daughter -> getExtraInfo("lifeTimeErr");
87 return daughter -> getExtraInfo("lifeTime");
88 }
89 }
90
91
92 double mumvtxX = particle->getX();
93 double mumvtxY = particle->getY();
94 double mumvtxZ = particle->getZ();
95 if (particle == daughter) {
96 if (hasRAVEBeamConstrainedProductionVertex(particle)) {
97 mumvtxX = particle->getExtraInfo("prodVertX");
98 mumvtxY = particle->getExtraInfo("prodVertY");
99 mumvtxZ = particle->getExtraInfo("prodVertZ");
100 } else {
101 //if no production vertex assume the particle originated at the ip
102 static DBObjPtr<BeamSpot> beamSpotDB;
103 mumvtxX = (beamSpotDB->getIPPosition()).X();
104 mumvtxY = (beamSpotDB->getIPPosition()).Y();
105 mumvtxZ = (beamSpotDB->getIPPosition()).Z();
106 }
107 }
108 //daughter vertex
109 double vtxX = daughter->getX();
110 double vtxY = daughter->getY();
111 double vtxZ = daughter->getZ();
112 // daughter MOMENTUM
113 double p = daughter->getP();
114 double pX = daughter->getPx();
115 double pY = daughter->getPy();
116 double pZ = daughter->getPz();
117
118 //versor of the daughter momentum
119 double nX = pX / p;
120 double nY = pY / p;
121 double nZ = pZ / p;
122 //mother vertex
123 double lX = vtxX - mumvtxX;
124 double lY = vtxY - mumvtxY;
125 double lZ = vtxZ - mumvtxZ;
126
127 //flight distance
128 double fD = lX * nX + lY * nY + lZ * nZ;
129 //flight time
130 double fT = daughter->getPDGMass() / Const::speedOfLight * fD / p;
131
132 //covariance matrix of momentum and vertex for the Dz
133 //ORDER = px,py,pz,E,x,y,z
134 TMatrixFSym dauCov = daughter->getMomentumVertexErrorMatrix();
135 TMatrixFSym mumCov = particle->getVertexErrorMatrix(); //order: x,y,z
136 if (particle == daughter) {
137 if (hasRAVEBeamConstrainedProductionVertex(particle)) {
138 std::vector<std::string> directions = {"x", "y", "z"};
139 for (unsigned int i = 0; i < directions.size(); i++) {
140 for (unsigned int j = 0; j < directions.size(); j++) {
141 mumCov[i][j] = particle->getExtraInfo(boost::str(boost::format("prodVertS%s%s") % directions[i] % directions[j]));
142 }
143 }
144 } else {
145 static DBObjPtr<BeamSpot> beamSpotDB;
146 mumCov = beamSpotDB->getCovVertex();
147 }
148 }
149 //compute total covariance matrix
150 //ORDER = px dau, py dau, pz dau, E dau, x dau, y dau, z dau, x mum, y mum, z mum
151
152 TMatrixFSym Cov(10);
153 for (int i = 0; i < 10; i++)
154 for (int j = 0; j < 10; j++)
155 if (i < 7 && j < 7)
156 Cov[i][j] = dauCov[i][j];
157 else if (i > 6 && j > 6)
158 Cov[i][j] = mumCov[i - 7][j - 7];
159 else
160 Cov[i][j] = 0;
161
162 if (mode == "distance") {
163 TMatrixF deriv(10, 1);
164 deriv[0][0] = (lX - nX * fD) / p; //px Daughter
165 deriv[1][0] = (lY - nY * fD) / p; //py Daughter
166 deriv[2][0] = (lZ - nZ * fD) / p; //pz Daughter
167 deriv[3][0] = 0; //E Daughter
168 deriv[4][0] = nX; //vtxX Daughter
169 deriv[5][0] = nY; //vtxY Daughter
170 deriv[6][0] = nZ; //vtxZ Daughter
171 deriv[7][0] = - nX; //vtxX Mother
172 deriv[8][0] = - nY; //vtxY Mother
173 deriv[9][0] = - nZ; //vtxZ Mother
174
175
176 TMatrixF tmp(10, 1);
177 tmp.Mult(Cov, deriv);
178
179 TMatrixF result(1, 1);
180 result.Mult(deriv.T(), tmp);
181
182 outErr = sqrt(result[0][0]);
183 return fD;
184 }
185 if (mode == "time") {
186 TMatrixF deriv(10, 1);
187 deriv[0][0] = (daughter->getPDGMass() / Const::speedOfLight * lX - 2 * pX * fT) / p / p; //px Daughter
188 deriv[1][0] = (daughter->getPDGMass() / Const::speedOfLight * lY - 2 * pY * fT) / p / p; //py Daughter
189 deriv[2][0] = (daughter->getPDGMass() / Const::speedOfLight * lZ - 2 * pZ * fT) / p / p; //pz Daughter
190 deriv[3][0] = 0; //E Daughter
191 deriv[4][0] = daughter->getPDGMass() / Const::speedOfLight * pX / p / p; //vtxX Daughter
192 deriv[5][0] = daughter->getPDGMass() / Const::speedOfLight * pY / p / p; //vtxY Daughter
193 deriv[6][0] = daughter->getPDGMass() / Const::speedOfLight * pZ / p / p; //vtxZ Daughter
194 deriv[7][0] = - daughter->getPDGMass() / Const::speedOfLight * pX / p / p; //vtxX Mother
195 deriv[8][0] = - daughter->getPDGMass() / Const::speedOfLight * pY / p / p; //vtxY Mother
196 deriv[9][0] = - daughter->getPDGMass() / Const::speedOfLight * pZ / p / p; //vtxZ Mother
197
198
199 TMatrixF tmp(10, 1);
200 tmp.Mult(Cov, deriv);
201
202 TMatrixF result(1, 1);
203 result.Mult(deriv.T(), tmp);
204 outErr = sqrt(result[0][0]);
205 return fT;
206 }
207 outErr = Const::doubleNaN;
208 return Const::doubleNaN;
209 }
210
211 // Helper function for MC flight time and distance
212 inline double getMCFlightInfoBtw(const MCParticle* mcparticle, const std::string& mode)
213 {
214
215 if (mcparticle == nullptr)
216 return Const::doubleNaN;
217
218
219 if (mode == "distance") {
220 ROOT::Math::XYZVector deltaVtx = mcparticle->getDecayVertex() - mcparticle->getProductionVertex();
221 double distance = deltaVtx.R();
222 if (distance < 0)
223 B2WARNING("Negative true flight distance, it's forbidden -> something went wrong.");
224
225 return distance;
226 }
227
228 if (mode == "time") {
229 double lifetime = mcparticle->getLifetime();
230 double mass = mcparticle->getMass();
231 double time = Const::doubleNaN;
232 if (mass == 0)
233 B2WARNING("you are asking for the proper time of a massless particle which is not allowed, returning -99.");
234 else {
235 double energy = mcparticle->getEnergy();
236 time = lifetime / energy * mass;
237 }
238
239 if (time < 0)
240 B2WARNING("Negative true proper time, it's forbidden -> something went wrong.");
241
242 return time;
243
244 }
245 B2WARNING("MCFlightInfo helper function called with mode '" << mode
246 << "'. Only 'distance' and 'time' are available.");
247 return Const::doubleNaN;
248
249 }
250
251 double flightDistance(const Particle* part)
252 {
253 double flightDistanceError = Const::doubleNaN;
254 return getFlightInfoBtw(part, part, flightDistanceError, "distance");
255 }
256
257 double flightTime(const Particle* part)
258 {
259 double flightTimeError = Const::doubleNaN;
260 return getFlightInfoBtw(part, part, flightTimeError, "time");
261 }
262
263 double flightDistanceErr(const Particle* part)
264 {
265 double flightDistanceError = Const::doubleNaN;
266 getFlightInfoBtw(part, part, flightDistanceError, "distance");
267 return flightDistanceError;
268 }
269
270 double flightTimeErr(const Particle* part)
271 {
272 double flightTimeError = Const::doubleNaN;
273 getFlightInfoBtw(part, part, flightTimeError, "time");
274 return flightTimeError;
275 }
276
277 inline double getVertexDistance(const Particle* particle, const Particle* daughter, double& vertexDistanceErr,
278 bool prodVertIsIP = false)
279 {
280 if (!particle || !daughter) {
281 vertexDistanceErr = Const::doubleNaN;
282 return Const::doubleNaN;
283 }
284
285 // production vertex
286 double prodVtxX = particle->getX();
287 double prodVtxY = particle->getY();
288 double prodVtxZ = particle->getZ();
289 if (particle == daughter || prodVertIsIP) {
290 if (particle->hasExtraInfo("prodVertX")) prodVtxX = particle->getExtraInfo("prodVertX");
291 if (particle->hasExtraInfo("prodVertY")) prodVtxY = particle->getExtraInfo("prodVertY");
292 if (particle->hasExtraInfo("prodVertZ")) prodVtxZ = particle->getExtraInfo("prodVertZ");
293 }
294
295 // decay vertex
296 double vtxX = daughter->getX();
297 double vtxY = daughter->getY();
298 double vtxZ = daughter->getZ();
299
300 // difference between vertices
301 double lX = vtxX - prodVtxX;
302 double lY = vtxY - prodVtxY;
303 double lZ = vtxZ - prodVtxZ;
304
305 // vertex distance
306 double lD = sqrt(lX * lX + lY * lY + lZ * lZ);
307
308 // covariance matrix of both vertices
309 TMatrixFSym decCov = daughter->getVertexErrorMatrix();
310 TMatrixFSym prodCov = particle->getVertexErrorMatrix();
311 if (particle == daughter || prodVertIsIP) {
312 if (particle->hasExtraInfo("prodVertSxx")) prodCov[0][0] = particle->getExtraInfo("prodVertSxx");
313 if (particle->hasExtraInfo("prodVertSxy")) prodCov[0][1] = particle->getExtraInfo("prodVertSxy");
314 if (particle->hasExtraInfo("prodVertSxz")) prodCov[0][2] = particle->getExtraInfo("prodVertSxz");
315 if (particle->hasExtraInfo("prodVertSyx")) prodCov[1][0] = particle->getExtraInfo("prodVertSyx");
316 if (particle->hasExtraInfo("prodVertSyy")) prodCov[1][1] = particle->getExtraInfo("prodVertSyy");
317 if (particle->hasExtraInfo("prodVertSyz")) prodCov[1][2] = particle->getExtraInfo("prodVertSyz");
318 if (particle->hasExtraInfo("prodVertSzx")) prodCov[2][0] = particle->getExtraInfo("prodVertSzx");
319 if (particle->hasExtraInfo("prodVertSzy")) prodCov[2][1] = particle->getExtraInfo("prodVertSzy");
320 if (particle->hasExtraInfo("prodVertSzz")) prodCov[2][2] = particle->getExtraInfo("prodVertSzz");
321 }
322
323 // compute total covariance matrix
324 // ORDER = prod x, prod y, prod z, decay x, decay y, decay z
325
326 TMatrixFSym Cov(6);
327 for (int i = 0; i < 6; i++)
328 for (int j = 0; j < 6; j++)
329 if (i < 3 && j < 3)
330 Cov[i][j] = prodCov[i][j];
331 else if (i > 2 && j > 2)
332 Cov[i][j] = decCov[i - 3][j - 3];
333 else
334 Cov[i][j] = 0;
335
336 TMatrixF deriv(6, 1);
337 deriv[0][0] = - lX / lD; // prodVtxX
338 deriv[1][0] = - lY / lD; // prodVtxY
339 deriv[2][0] = - lZ / lD; // prodVtxZ
340 deriv[3][0] = lX / lD; // vtxX
341 deriv[4][0] = lY / lD; // vtxY
342 deriv[5][0] = lZ / lD; // vtxZ
343
344
345 TMatrixF tmp(6, 1);
346 tmp.Mult(Cov, deriv);
347
348 TMatrixF result(1, 1);
349 result.Mult(deriv.T(), tmp);
350
351 vertexDistanceErr = sqrt(result[0][0]);
352 return lD;
353 }
354
355 double vertexDistance(const Particle* part)
356 {
357 double vertexDistanceError = Const::doubleNaN;
358 if (!part->hasExtraInfo("prodVertX") || !part->hasExtraInfo("prodVertY") || !part->hasExtraInfo("prodVertZ")) {
359 return Const::doubleNaN;
360 }
361 return getVertexDistance(part, part, vertexDistanceError);
362 }
363
364 double vertexDistanceErr(const Particle* part)
365 {
366 double vertexDistanceError = Const::doubleNaN;
367 if (!part->hasExtraInfo("prodVertX") || !part->hasExtraInfo("prodVertY") || !part->hasExtraInfo("prodVertZ")) {
368 return Const::doubleNaN;
369 }
370 getVertexDistance(part, part, vertexDistanceError);
371 return vertexDistanceError;
372 }
373
374 double vertexDistanceSignificance(const Particle* part)
375 {
376 double vertexDistanceError = Const::doubleNaN;
377 if (!part->hasExtraInfo("prodVertX") || !part->hasExtraInfo("prodVertY") || !part->hasExtraInfo("prodVertZ")) {
378 return Const::doubleNaN;
379 }
380 return getVertexDistance(part, part, vertexDistanceError) / vertexDistanceError;
381 }
382
383 double flightTimeOfDaughter(const Particle* particle, const std::vector<double>& daughters)
384 {
385 if (!particle)
386 return Const::doubleNaN; // Initial particle is NULL
387
388 long daughterNumber = -1;
389 if (daughters.size() > 0) {
390 daughterNumber = std::lround(daughters[0]);
391 } else {
392 B2ERROR("At least one integer, the index of the daughter, must be provided to flightTimeOfDaughter!");
393 return Const::doubleNaN;
394 }
395 int nDaughters = static_cast<int>(particle->getNDaughters());
396 if (daughterNumber >= nDaughters) {
397 B2ERROR("The daughter index provided to flightTimeOfDaughter is larger than the number of daughters of this particle!");
398 return Const::doubleNaN;
399 }
400
401 long grandDaughterNumber = -1;
402 if (daughters.size() == 2) {
403 grandDaughterNumber = std::lround(daughters[1]);
404 }
405
406 const Particle* daughter = particle->getDaughter(daughterNumber);
407 double flightTimeError;
408 if (grandDaughterNumber > -1) {
409 if (grandDaughterNumber < (int)daughter->getNDaughters()) {
410 return getFlightInfoBtw(particle, daughter->getDaughter(grandDaughterNumber), flightTimeError, "time", true);
411 } else {
412 B2ERROR("The granddaughter index provided to flightTimeOfDaughter is too large!");
413 return Const::doubleNaN;
414 }
415 } else {
416 return getFlightInfoBtw(particle, daughter, flightTimeError, "time");
417 }
418 }
419
420 // Flight time uncertainty
421 double flightTimeOfDaughterErr(const Particle* particle, const std::vector<double>& daughters)
422 {
423 if (!particle)
424 return Const::doubleNaN; // Initial particle is NULL
425
426 long daughterNumber = -1;
427 if (daughters.size() > 0) {
428 daughterNumber = std::lround(daughters[0]);
429 } else {
430 B2ERROR("At least one integer, the index of the daughter, must be provided to flightTimeOfDaughterErr!");
431 return Const::doubleNaN;
432 }
433 int nDaughters = static_cast<int>(particle->getNDaughters());
434 if (daughterNumber >= nDaughters) {
435 B2ERROR("The daughter index provided to flightTimeOfDaughterErr is larger than the number of daughters of this particle!");
436 return Const::doubleNaN;
437 }
438
439 long grandDaughterNumber = -1;
440 if (daughters.size() == 2) {
441 grandDaughterNumber = std::lround(daughters[1]);
442 }
443
444 const Particle* daughter = particle->getDaughter(daughterNumber);
445 double flightTimeError = Const::doubleNaN;;
446 if (grandDaughterNumber > -1) {
447 if (grandDaughterNumber < (int)daughter->getNDaughters()) {
448 getFlightInfoBtw(particle, daughter->getDaughter(grandDaughterNumber), flightTimeError, "time", true);
449 } else {
450 B2ERROR("The granddaughter index provided to flightTimeOfDaughterErr is too large!");
451 }
452 } else {
453 getFlightInfoBtw(particle, daughter, flightTimeError, "time");
454 }
455 return flightTimeError;
456 }
457
458 double flightDistanceOfDaughter(const Particle* particle, const std::vector<double>& daughters)
459 {
460 if (!particle)
461 return Const::doubleNaN; // Initial particle is NULL
462
463 long daughterNumber = -1;
464 if (daughters.size() > 0) {
465 daughterNumber = std::lround(daughters[0]);
466 } else {
467 B2ERROR("At least one integer, the index of the daughter, must be provided to flightDistanceOfDaughter!");
468 return Const::doubleNaN;
469 }
470 int nDaughters = static_cast<int>(particle->getNDaughters());
471 if (daughterNumber >= nDaughters) {
472 B2ERROR("The daughter index provided to flightDistanceOfDaughter is larger than the number of daughters of this particle!");
473 return Const::doubleNaN;
474 }
475
476 long grandDaughterNumber = -1;
477 if (daughters.size() == 2) {
478 grandDaughterNumber = std::lround(daughters[1]);
479 }
480
481 const Particle* daughter = particle->getDaughter(daughterNumber);
482 double flightDistanceError;
483 if (grandDaughterNumber > -1) {
484 if (grandDaughterNumber < (int)daughter->getNDaughters()) {
485 return getFlightInfoBtw(particle, daughter->getDaughter(grandDaughterNumber), flightDistanceError, "distance", true);
486 } else {
487 B2ERROR("The granddaughter index provided to flightDistanceOfDaughter is too large!");
488 return Const::doubleNaN;
489 }
490 } else {
491 return getFlightInfoBtw(particle, daughter, flightDistanceError, "distance");
492 }
493 }
494
495 // Flight distance uncertainty
496 double flightDistanceOfDaughterErr(const Particle* particle, const std::vector<double>& daughters)
497 {
498 if (!particle)
499 return Const::doubleNaN; // Initial particle is NULL
500
501 long daughterNumber = -1;
502 if (daughters.size() > 0) {
503 daughterNumber = std::lround(daughters[0]);
504 } else {
505 B2ERROR("At least one integer, the index of the daughter, must be provided to flightDistanceOfDaughterErr!");
506 return Const::doubleNaN;
507 }
508 int nDaughters = static_cast<int>(particle->getNDaughters());
509 if (daughterNumber >= nDaughters) {
510 B2ERROR("The daughter index provided to flightDistanceOfDaughterErr is larger than the number of daughters of this particle!");
511 return Const::doubleNaN;
512 }
513
514 long grandDaughterNumber = -1;
515 if (daughters.size() == 2) {
516 grandDaughterNumber = std::lround(daughters[1]);
517 }
518
519 const Particle* daughter = particle->getDaughter(daughterNumber);
520 double flightDistanceError = Const::doubleNaN;
521 if (grandDaughterNumber > -1) {
522 if (grandDaughterNumber < (int)daughter->getNDaughters()) {
523 getFlightInfoBtw(particle, daughter->getDaughter(grandDaughterNumber), flightDistanceError, "distance", true);
524 } else {
525 B2ERROR("The granddaughter index provided to flightDistanceOfDaughterErr is too large!");
526 }
527 } else {
528 getFlightInfoBtw(particle, daughter, flightDistanceError, "distance");
529 }
530 return flightDistanceError;
531 }
532
533 // Distance between mother and daughter vertices
534 double vertexDistanceOfDaughter(const Particle* particle, const std::vector<double>& arguments)
535 {
536 if (!particle)
537 return Const::doubleNaN; // Initial particle is NULL
538 long daughterNumber = -1;
539 if (arguments.size() > 0) {
540 daughterNumber = std::lround(arguments[0]);
541 } else {
542 B2ERROR("At least one integer, the index of the daughter, must be provided to vertexDistanceOfDaughter!");
543 return Const::doubleNaN;
544 }
545 int nDaughters = static_cast<int>(particle->getNDaughters());
546 if (daughterNumber >= nDaughters) {
547 B2ERROR("The daughter index provided to vertexDistanceOfDaughter is larger than the number of daughters of this particle!");
548 return Const::doubleNaN;
549 }
550
551 bool prodVertIsIP = true;
552 if (arguments.size() == 2) {
553 prodVertIsIP = false;
554 }
555
556 const Particle* daughter = particle->getDaughter(daughterNumber);
557 double vertexDistanceError;
558 return getVertexDistance(particle, daughter, vertexDistanceError, prodVertIsIP);
559 }
560
561 // Uncertainty on distance between mother and daughter vertices
562 double vertexDistanceOfDaughterErr(const Particle* particle, const std::vector<double>& arguments)
563 {
564 if (!particle)
565 return Const::doubleNaN; // Initial particle is NULL
566 long daughterNumber = -1;
567 if (arguments.size() > 0) {
568 daughterNumber = std::lround(arguments[0]);
569 } else {
570 B2ERROR("At least one integer, the index of the daughter, must be provided to vertexDistanceOfDaughterErr!");
571 return Const::doubleNaN;
572 }
573 int nDaughters = static_cast<int>(particle->getNDaughters());
574 if (daughterNumber >= nDaughters) {
575 B2ERROR("The daughter index provided to vertexDistanceOfDaughterErr is larger than the number of daughters of this particle!");
576 return Const::doubleNaN;
577 }
578
579 bool prodVertIsIP = true;
580 if (arguments.size() == 2) {
581 prodVertIsIP = false;
582 }
583
584 const Particle* daughter = particle->getDaughter(daughterNumber);
585 double vertexDistanceError;
586 getVertexDistance(particle, daughter, vertexDistanceError, prodVertIsIP);
587 return vertexDistanceError;
588 }
589
590 // Significance of distance between mother and daughter vertices
591 double vertexDistanceOfDaughterSignificance(const Particle* particle, const std::vector<double>& arguments)
592 {
593 if (!particle)
594 return Const::doubleNaN; // Initial particle is NULL
595 long daughterNumber = -1;
596 if (arguments.size() > 0) {
597 daughterNumber = std::lround(arguments[0]);
598 } else {
599 B2ERROR("At least one integer, the index of the daughter, must be provided to vertexDistanceOfDaughterSignificance!");
600 return Const::doubleNaN;
601 }
602 int nDaughters = static_cast<int>(particle->getNDaughters());
603 if (daughterNumber >= nDaughters) {
604 B2ERROR("The daughter index provided to vertexDistanceOfDaughterSignificance is larger than the number of daughters of this particle!");
605 return Const::doubleNaN;
606 }
607
608 bool prodVertIsIP = true;
609 if (arguments.size() == 2) {
610 prodVertIsIP = false;
611 }
612
613 const Particle* daughter = particle->getDaughter(daughterNumber);
614 double vertexDistanceError;
615 return getVertexDistance(particle, daughter, vertexDistanceError, prodVertIsIP) / vertexDistanceError;
616 }
617
618 // MC variables
619
620 double mcFlightDistance(const Particle* particle)
621 {
622 if (particle == nullptr)
623 return Const::doubleNaN; // Initial particle is NULL
624
625 const MCParticle* mcparticle = particle->getMCParticle();
626
627 if (mcparticle == nullptr)
628 return Const::doubleNaN; // Initial particle is NULL
629
630 return getMCFlightInfoBtw(mcparticle, "distance");
631 }
632
633 double mcFlightTime(const Particle* particle)
634 {
635 if (particle == nullptr)
636 return Const::doubleNaN; // Initial particle is NULL
637
638 const MCParticle* mcparticle = particle->getMCParticle();
639
640 if (mcparticle == nullptr)
641 return Const::doubleNaN; // Initial particle is NULL
642
643
644 return getMCFlightInfoBtw(mcparticle, "time");
645 }
646
647 double mcFlightDistanceOfDaughter(const Particle* particle, const std::vector<double>& daughters)
648 {
649 if (!particle)
650 return Const::doubleNaN; // Initial particle is NULL
651
652 long daughterNumber = -1;
653 if (daughters.size() > 0) {
654 daughterNumber = std::lround(daughters[0]);
655 } else {
656 B2ERROR("At least one integer, the index of the daughter, must be provided to mcFlightDistanceOfDaughter!");
657 return Const::doubleNaN;
658 }
659 int nDaughters = static_cast<int>(particle->getNDaughters());
660 if (daughterNumber >= nDaughters) {
661 B2ERROR("The daughter index provided to mcFlightDistanceOfDaughter is larger than the number of daughters of this particle!");
662 return Const::doubleNaN;
663 }
664
665 long grandDaughterNumber = -1;
666 if (daughters.size() == 2) {
667 grandDaughterNumber = std::lround(daughters[1]);
668 }
669
670 const Particle* daughterReco = particle->getDaughter(daughterNumber);
671 //get the MC DAUGHTER
672 const MCParticle* daughter = daughterReco->getMCParticle();
673
674 double flightDistanceMC;
675 if (grandDaughterNumber > -1 && grandDaughterNumber < (int)daughterReco->getNDaughters()) {
676 // Compute value between mother and granddaughter
677 const MCParticle* gdaughter = daughterReco->getDaughter(grandDaughterNumber)->getMCParticle();
678 flightDistanceMC = getMCFlightInfoBtw(gdaughter, "distance");
679 } else {
680 // Compute value between mother and daughter
681 flightDistanceMC = getMCFlightInfoBtw(daughter, "distance");
682 }
683 return flightDistanceMC;
684 }
685
686 double mcFlightTimeOfDaughter(const Particle* particle, const std::vector<double>& daughters)
687 {
688 if (!particle)
689 return Const::doubleNaN; // Initial particle is NULL
690
691 long daughterNumber = -1;
692 if (daughters.size() > 0) {
693 daughterNumber = std::lround(daughters[0]);
694 } else {
695 B2ERROR("At least one integer, the index of the daughter, must be provided to mcFlightTimeOfDaughter!");
696 return Const::doubleNaN;
697 }
698 int nDaughters = static_cast<int>(particle->getNDaughters());
699 if (daughterNumber >= nDaughters) {
700 B2ERROR("The daughter index provided to mcFlightTimeOfDaughter is larger than the number of daughters of this particle!");
701 return Const::doubleNaN;
702 }
703
704 long grandDaughterNumber = -1;
705 if (daughters.size() == 2) {
706 grandDaughterNumber = std::lround(daughters[1]);
707 }
708
709 const Particle* daughterReco = particle->getDaughter(daughterNumber);
710 //get the MC DAUGHTER
711 const MCParticle* daughter = daughterReco->getMCParticle();
712
713 double flightTimeMC;
714 if (grandDaughterNumber > -1 && grandDaughterNumber < (int)daughterReco->getNDaughters()) {
715 // Compute value between mother and granddaughter
716 const MCParticle* gdaughter = daughterReco->getDaughter(grandDaughterNumber)->getMCParticle();
717 flightTimeMC = getMCFlightInfoBtw(gdaughter, "time");
718 } else {
719 flightTimeMC = getMCFlightInfoBtw(daughter, "time");
720 }
721 return flightTimeMC;
722 }
723
724
725 VARIABLE_GROUP("Flight Information");
726 REGISTER_VARIABLE("flightTime", flightTime,
727 "Returns the flight time of particle. If a treeFit has been performed the flight time calculated by TreeFitter is returned. Otherwise if a beam constrained rave fit has been performed the production vertex set by rave and the decay vertex are used to calculate the flight time. If neither fit has been performed the i.p. is taken to be the production vertex.\n\n",
728 "ns");
729 REGISTER_VARIABLE("flightDistance", flightDistance,
730 "Returns the flight distance of particle. If a treeFit has been performed the flight distance calculated by TreeFitter is returned. Otherwise if a beam constrained rave fit has been performed the production vertex set by rave and the decay vertex are used to calculate the flight distance. If neither fit has been performed the i.p. is taken to be the production vertex.\n\n",
731 "cm");
732 REGISTER_VARIABLE("flightTimeErr", flightTimeErr,
733 "Returns the flight time error of particle. If a treeFit has been performed the flight time error calculated by TreeFitter is returned. Otherwise if a beam constrained rave fit has been performed the production vertex set by rave and the decay vertex are used to calculate the flight time error. If neither fit has been performed the i.p. is taken to be the production vertex.\n\n",
734 "ns");
735 REGISTER_VARIABLE("flightDistanceErr", flightDistanceErr,
736 "Returns the flight distance error of particle. If a treeFit has been performed the flight distance error calculated by TreeFitter is returned. Otherwise if a beam constrained rave fit has been performed the production vertex set by rave and the decay vertex are used to calculate the flight distance error. If neither fit has been performed the i.p. is taken to be the production vertex.\n\n",
737 "cm");
738 // Daughters
739 REGISTER_VARIABLE("flightTimeOfDaughter(daughterN, gdaughterN = -1)", flightTimeOfDaughter,
740 "Returns the flight time between mother and daughter particle with daughterN index. If a treeFit has been performed the value calculated by treeFitter is returned. Otherwise the value is calculated using the decay vertices of the mother and daughter particle. If a second index granddaughterM is given the value is calculated between the mother and the Mth grandaughter (Mth daughter of Nth daughter).\n\n",
741 "ns");
742 REGISTER_VARIABLE("flightTimeOfDaughterErr(daughterN, gdaughterN = -1)", flightTimeOfDaughterErr,
743 "Returns the flight time error between mother and daughter particle with daughterN index. If a treeFit has been performed the value calculated by treeFitter is returned. Otherwise the value is calculated using the decay vertices of the mother and daughter particle. If a second index granddaughterM is given the value is calculated between the mother and the Mth grandaughter (Mth daughter of Nth daughter).\n\n",
744 "ns");
745 REGISTER_VARIABLE("flightDistanceOfDaughter(daughterN, gdaughterN = -1)", flightDistanceOfDaughter,
746 "Returns the flight distance between mother and daughter particle with daughterN index. If a treeFit has been performed the value calculated by treeFitter is returned. Otherwise the value is calculated using the decay vertices of the mother and daughter particle. If a second index granddaughterM is given the value is calculated between the mother and the Mth grandaughter (Mth daughter of Nth daughter).\n\n",
747 "cm");
748 REGISTER_VARIABLE("flightDistanceOfDaughterErr(daughterN, gdaughterN = -1)", flightDistanceOfDaughterErr,
749 "Returns the flight distance error between mother and daughter particle with daughterN index. If a treeFit has been performed the value calculated by treeFitter is returned. Otherwise the value is calculated using the decay vertices of the mother and daughter particle. If a second index granddaughterM is given the value is calculated between the mother and the Mth grandaughter (Mth daughter of Nth daughter).\n\n",
750 "cm");
751 // MC Info
752 REGISTER_VARIABLE("mcFlightDistance", mcFlightDistance,
753 "Returns the MC flight distance of the particle\n\n", "cm");
754 REGISTER_VARIABLE("mcFlightTime", mcFlightTime,
755 "Returns the MC flight time of the particle\n\n", "ns");
756 REGISTER_VARIABLE("mcFlightDistanceOfDaughter(daughterN, gdaughterN = -1)", mcFlightDistanceOfDaughter,
757 "Returns the MC flight distance between mother and daughter particle using generated info\n\n", "cm");
758 REGISTER_VARIABLE("mcFlightTimeOfDaughter(daughterN, gdaughterN = -1)", mcFlightTimeOfDaughter,
759 "Returns the MC flight time between mother and daughter particle using generated info\n\n", "ns");
760 //Vertex Distance
761 REGISTER_VARIABLE("vertexDistance", vertexDistance,
762 "Returns the distance between the production and decay vertex of a particle. Returns NaN if particle has no production or decay vertex.\n\n",
763 "cm");
764 REGISTER_VARIABLE("vertexDistanceErr", vertexDistanceErr,
765 "Returns the uncertainty on the distance between the production and decay vertex of a particle. Returns NaN if particle has no production or decay vertex.\n\n",
766 "cm");
767 REGISTER_VARIABLE("vertexDistanceSignificance", vertexDistanceSignificance,
768 "Returns the distance between the production and decay vertex of a particle in units of the uncertainty on this value, i.e. the significance of the vertex separation.");
769 REGISTER_VARIABLE("vertexDistanceOfDaughter(daughterN[, option])", vertexDistanceOfDaughter,
770 "If any integer is provided as second argument it returns the distance between the decay vertices of the particle and of its daughter with index daughterN.\n"
771 "Otherwise, it is assumed that the particle has a production vertex (typically the IP) which is used to calculate the distance to the daughter's decay vertex.\n"
772 "Returns NaN in case anything goes wrong.\n\n", "cm");
773 REGISTER_VARIABLE("vertexDistanceOfDaughterErr(daughterN[, option])", vertexDistanceOfDaughterErr,
774 "If any integer is provided as second argument it returns the uncertainty on the distance between the decay vertices of the particle and of its daughter with index daughterN.\n"
775 "Otherwise, it is assumed that the particle has a production vertex (typically the IP) with a corresponding covariance matrix to calculate the uncertainty on the distance to the daughter's decay vertex.\n"
776 "Returns NaN in case anything goes wrong.\n\n", "cm");
777 REGISTER_VARIABLE("vertexDistanceOfDaughterSignificance(daughterN[, option)", vertexDistanceOfDaughterSignificance,
778 "If any integer is provided as second argument it returns the distance between the decay vertices of the particle and of its daughter with index daughterN in units of the uncertainty on this value.\n"
779 "Otherwise, it is assumed that the particle has a production vertex (typically the IP) with a corresponding covariance matrix and the significance of the separation to this vertex is calculated.");
780 }
782} // Belle2 namespace
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
static const double doubleNaN
quiet_NaN
Definition: Const.h:703
const MCParticle * getMCParticle() const
Returns the pointer to the MCParticle object that was used to create this Particle (ParticleType == c...
Definition: Particle.cc:942
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition: Particle.cc:631
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Abstract base class for different kinds of events.