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