Belle II Software  release-06-01-15
Variables.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 include
10 #include <analysis/variables/Variables.h>
11 #include <analysis/VariableManager/Manager.h>
12 #include <analysis/utility/PCmsLabTransform.h>
13 #include <analysis/utility/ReferenceFrame.h>
14 
15 #include <analysis/utility/MCMatching.h>
16 #include <analysis/ClusterUtility/ClusterUtils.h>
17 
18 // framework - DataStore
19 #include <framework/datastore/StoreArray.h>
20 #include <framework/datastore/StoreObjPtr.h>
21 
22 // dataobjects
23 #include <analysis/dataobjects/Particle.h>
24 #include <analysis/dataobjects/EventExtraInfo.h>
25 #include <analysis/dataobjects/EventShapeContainer.h>
26 
27 #include <mdst/dataobjects/MCParticle.h>
28 #include <mdst/dataobjects/Track.h>
29 #include <mdst/dataobjects/ECLCluster.h>
30 
31 #include <mdst/dbobjects/BeamSpot.h>
32 
33 // framework aux
34 #include <framework/logging/Logger.h>
35 #include <framework/geometry/BFieldManager.h>
36 #include <framework/gearbox/Const.h>
37 
38 #include <TLorentzVector.h>
39 #include <TRandom.h>
40 #include <TVectorF.h>
41 #include <TVector3.h>
42 
43 #include <iostream>
44 #include <cmath>
45 
46 using namespace std;
47 
48 namespace Belle2 {
53  namespace Variable {
54 
55 // momentum -------------------------------------------
56 
57  double particleP(const Particle* part)
58  {
59  const auto& frame = ReferenceFrame::GetCurrent();
60  return frame.getMomentum(part).P();
61  }
62 
63  double particleE(const Particle* part)
64  {
65  const auto& frame = ReferenceFrame::GetCurrent();
66  return frame.getMomentum(part).E();
67  }
68 
69  double particleClusterEUncertainty(const Particle* part)
70  {
71  const ECLCluster* cluster = part->getECLCluster();
72  if (cluster) {
73  const auto EPhiThetaCov = cluster->getCovarianceMatrix3x3();
74  return std::sqrt(EPhiThetaCov[0][0]);
75  }
76  return std::numeric_limits<double>::quiet_NaN();
77  }
78 
79  double particlePx(const Particle* part)
80  {
81  const auto& frame = ReferenceFrame::GetCurrent();
82  return frame.getMomentum(part).Px();
83  }
84 
85  double particlePy(const Particle* part)
86  {
87  const auto& frame = ReferenceFrame::GetCurrent();
88  return frame.getMomentum(part).Py();
89  }
90 
91  double particlePz(const Particle* part)
92  {
93  const auto& frame = ReferenceFrame::GetCurrent();
94  return frame.getMomentum(part).Pz();
95  }
96 
97  double particlePt(const Particle* part)
98  {
99  const auto& frame = ReferenceFrame::GetCurrent();
100  return frame.getMomentum(part).Pt();
101  }
102 
103  double covMatrixElement(const Particle* part, const std::vector<double>& element)
104  {
105  int elementI = int(std::lround(element[0]));
106  int elementJ = int(std::lround(element[1]));
107 
108  if (elementI < 0 || elementI > 6) {
109  B2WARNING("Requested particle's momentumVertex covariance matrix element is out of boundaries [0 - 6]: i = " << elementI);
110  return std::numeric_limits<double>::quiet_NaN();
111  }
112  if (elementJ < 0 || elementJ > 6) {
113  B2WARNING("Requested particle's momentumVertex covariance matrix element is out of boundaries [0 - 6]: j = " << elementJ);
114  return std::numeric_limits<double>::quiet_NaN();
115  }
116 
117  return part->getMomentumVertexErrorMatrix()(elementI, elementJ);
118  }
119 
120  double particleEUncertainty(const Particle* part)
121  {
122  const auto& frame = ReferenceFrame::GetCurrent();
123 
124  double errorSquared = frame.getMomentumErrorMatrix(part)(3, 3);
125 
126  if (errorSquared >= 0.0)
127  return sqrt(errorSquared);
128  else
129  return std::numeric_limits<double>::quiet_NaN();
130  }
131 
132  double particlePErr(const Particle* part)
133  {
134  TMatrixD jacobianRot(3, 3);
135  jacobianRot.Zero();
136 
137  double cosPhi = cos(particlePhi(part));
138  double sinPhi = sin(particlePhi(part));
139  double cosTheta = particleCosTheta(part);
140  double sinTheta = sin(acos(cosTheta));
141  double p = particleP(part);
142 
143  jacobianRot(0, 0) = sinTheta * cosPhi;
144  jacobianRot(0, 1) = sinTheta * sinPhi;
145  jacobianRot(1, 0) = cosTheta * cosPhi / p;
146  jacobianRot(1, 1) = cosTheta * sinPhi / p;
147  jacobianRot(0, 2) = cosTheta;
148  jacobianRot(2, 0) = -sinPhi / sinTheta / p;
149  jacobianRot(1, 2) = -sinTheta / p;
150  jacobianRot(2, 1) = cosPhi / sinTheta / p;
151 
152  const auto& frame = ReferenceFrame::GetCurrent();
153 
154  double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2, " ").Similarity(jacobianRot)(0, 0);
155 
156  if (errorSquared >= 0.0)
157  return sqrt(errorSquared);
158  else
159  return std::numeric_limits<double>::quiet_NaN();
160  }
161 
162  double particlePxErr(const Particle* part)
163  {
164  const auto& frame = ReferenceFrame::GetCurrent();
165 
166  double errorSquared = frame.getMomentumErrorMatrix(part)(0, 0);
167 
168  if (errorSquared >= 0.0)
169  return sqrt(errorSquared);
170  else
171  return std::numeric_limits<double>::quiet_NaN();
172  }
173 
174  double particlePyErr(const Particle* part)
175  {
176  const auto& frame = ReferenceFrame::GetCurrent();
177  double errorSquared = frame.getMomentumErrorMatrix(part)(1, 1);
178 
179  if (errorSquared >= 0.0)
180  return sqrt(errorSquared);
181  else
182  return std::numeric_limits<double>::quiet_NaN();
183  }
184 
185  double particlePzErr(const Particle* part)
186  {
187  const auto& frame = ReferenceFrame::GetCurrent();
188  double errorSquared = frame.getMomentumErrorMatrix(part)(2, 2);
189 
190  if (errorSquared >= 0.0)
191  return sqrt(errorSquared);
192  else
193  return std::numeric_limits<double>::quiet_NaN();
194  }
195 
196  double particlePtErr(const Particle* part)
197  {
198  TMatrixD jacobianRot(3, 3);
199  jacobianRot.Zero();
200 
201  double px = particlePx(part);
202  double py = particlePy(part);
203  double pt = particlePt(part);
204 
205  jacobianRot(0, 0) = px / pt;
206  jacobianRot(0, 1) = py / pt;
207  jacobianRot(1, 0) = -py / (pt * pt);
208  jacobianRot(1, 1) = px / (pt * pt);
209  jacobianRot(2, 2) = 1;
210 
211  const auto& frame = ReferenceFrame::GetCurrent();
212  double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2, " ").Similarity(jacobianRot)(0, 0);
213 
214  if (errorSquared >= 0.0)
215  return sqrt(errorSquared);
216  else
217  return std::numeric_limits<double>::quiet_NaN();
218  }
219 
220  double momentumDeviationChi2(const Particle* part)
221  {
222  double result = std::numeric_limits<double>::quiet_NaN();
223 
224  // check if error matrix is set
225  if (part->getPValue() < 0.0)
226  return result;
227 
228  // check if mc match exists
229  const MCParticle* mcp = part->getMCParticle();
230  if (mcp == nullptr)
231  return result;
232 
233  result = 0.0;
234  result += TMath::Power(part->getPx() - mcp->getMomentum().Px(), 2.0) / part->getMomentumVertexErrorMatrix()(0, 0);
235  result += TMath::Power(part->getPy() - mcp->getMomentum().Py(), 2.0) / part->getMomentumVertexErrorMatrix()(1, 1);
236  result += TMath::Power(part->getPz() - mcp->getMomentum().Pz(), 2.0) / part->getMomentumVertexErrorMatrix()(2, 2);
237 
238  return result;
239  }
240 
241  double particleTheta(const Particle* part)
242  {
243  const auto& frame = ReferenceFrame::GetCurrent();
244  return acos(frame.getMomentum(part).CosTheta());
245  }
246 
247  double particleThetaErr(const Particle* part)
248  {
249  TMatrixD jacobianRot(3, 3);
250  jacobianRot.Zero();
251 
252  double cosPhi = cos(particlePhi(part));
253  double sinPhi = sin(particlePhi(part));
254  double cosTheta = particleCosTheta(part);
255  double sinTheta = sin(acos(cosTheta));
256  double p = particleP(part);
257 
258  jacobianRot(0, 0) = sinTheta * cosPhi;
259  jacobianRot(0, 1) = sinTheta * sinPhi;
260  jacobianRot(1, 0) = cosTheta * cosPhi / p;
261  jacobianRot(1, 1) = cosTheta * sinPhi / p;
262  jacobianRot(0, 2) = cosTheta;
263  jacobianRot(2, 0) = -sinPhi / sinTheta / p;
264  jacobianRot(1, 2) = -sinTheta / p;
265  jacobianRot(2, 1) = cosPhi / sinTheta / p;
266 
267  const auto& frame = ReferenceFrame::GetCurrent();
268  double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2, " ").Similarity(jacobianRot)(1, 1);
269 
270  if (errorSquared >= 0.0)
271  return sqrt(errorSquared);
272  else
273  return std::numeric_limits<double>::quiet_NaN();
274  }
275 
276  double particleCosTheta(const Particle* part)
277  {
278  const auto& frame = ReferenceFrame::GetCurrent();
279  return frame.getMomentum(part).CosTheta();
280  }
281 
282  double particleCosThetaErr(const Particle* part)
283  {
284  return fabs(particleThetaErr(part) * sin(particleTheta(part)));
285  }
286 
287  double particlePhi(const Particle* part)
288  {
289  const auto& frame = ReferenceFrame::GetCurrent();
290  return frame.getMomentum(part).Phi();
291  }
292 
293  double particlePhiErr(const Particle* part)
294  {
295  TMatrixD jacobianRot(3, 3);
296  jacobianRot.Zero();
297 
298  double px = particlePx(part);
299  double py = particlePy(part);
300  double pt = particlePt(part);
301 
302  jacobianRot(0, 0) = px / pt;
303  jacobianRot(0, 1) = py / pt;
304  jacobianRot(1, 0) = -py / (pt * pt);
305  jacobianRot(1, 1) = px / (pt * pt);
306  jacobianRot(2, 2) = 1;
307 
308  const auto& frame = ReferenceFrame::GetCurrent();
309  double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2, " ").Similarity(jacobianRot)(1, 1);
310 
311  if (errorSquared >= 0.0)
312  return sqrt(errorSquared);
313  else
314  return std::numeric_limits<double>::quiet_NaN();
315  }
316 
317  double particleXp(const Particle* part)
318  {
319  PCmsLabTransform T;
320  TLorentzVector p4 = part -> get4Vector();
321  TLorentzVector p4CMS = T.rotateLabToCms() * p4;
322  float s = T.getCMSEnergy();
323  float M = part->getMass();
324  return p4CMS.P() / TMath::Sqrt(s * s / 4 - M * M);
325  }
326 
327  double particlePDGCode(const Particle* part)
328  {
329  return part->getPDGCode();
330  }
331 
332  double cosAngleBetweenMomentumAndVertexVectorInXYPlane(const Particle* part)
333  {
334  static DBObjPtr<BeamSpot> beamSpotDB;
335  double px = part->getMomentum().Px();
336  double py = part->getMomentum().Py();
337 
338  double xV = part->getVertex().X();
339  double yV = part->getVertex().Y();
340 
341  double xIP = (beamSpotDB->getIPPosition()).X();
342  double yIP = (beamSpotDB->getIPPosition()).Y();
343 
344  double x = xV - xIP;
345  double y = yV - yIP;
346 
347  double cosangle = (px * x + py * y) / (sqrt(px * px + py * py) * sqrt(x * x + y * y));
348  return cosangle;
349  }
350 
351  double cosAngleBetweenMomentumAndVertexVector(const Particle* part)
352  {
353  static DBObjPtr<BeamSpot> beamSpotDB;
354  return std::cos((part->getVertex() - beamSpotDB->getIPPosition()).Angle(part->getMomentum()));
355  }
356 
357  double cosThetaBetweenParticleAndNominalB(const Particle* part)
358  {
359 
360  int particlePDG = abs(part->getPDGCode());
361  if (particlePDG != 511 and particlePDG != 521)
362  B2FATAL("The Variables cosThetaBetweenParticleAndNominalB is only meant to be used on B mesons!");
363 
364  PCmsLabTransform T;
365  double e_Beam = T.getCMSEnergy() / 2.0; // GeV
366  double m_B = part->getPDGMass();
367  // if this is a continuum run, use an approximate Y(4S) CMS energy
368  if (e_Beam * e_Beam - m_B * m_B < 0) {
369  e_Beam = 1.0579400E+1 / 2.0; // GeV
370  }
371  double p_B = std::sqrt(e_Beam * e_Beam - m_B * m_B);
372 
373  TLorentzVector p = T.rotateLabToCms() * part->get4Vector();
374  double e_d = p.E();
375  double m_d = p.M();
376  double p_d = p.Rho();
377 
378  double theta_BY = (2 * e_Beam * e_d - m_B * m_B - m_d * m_d)
379  / (2 * p_B * p_d);
380  return theta_BY;
381  }
382 
383  double cosToThrustOfEvent(const Particle* part)
384  {
385  StoreObjPtr<EventShapeContainer> evtShape;
386  if (!evtShape) {
387  B2WARNING("Cannot find thrust of event information, did you forget to load the event shape calculation?");
388  return std::numeric_limits<float>::quiet_NaN();
389  }
390  PCmsLabTransform T;
391  TVector3 th = evtShape->getThrustAxis();
392  TVector3 particleMomentum = (T.rotateLabToCms() * part -> get4Vector()).Vect();
393  return std::cos(th.Angle(particleMomentum));
394  }
395 
396  double ImpactXY(const Particle* particle)
397  {
398  static DBObjPtr<BeamSpot> beamSpotDB;
399 
400  TVector3 mom = particle->getMomentum();
401 
402  TVector3 r = particle->getVertex() - beamSpotDB->getIPPosition();
403 
404  TVector3 Bfield = BFieldManager::getInstance().getFieldInTesla(beamSpotDB->getIPPosition());
405 
406  TVector3 curvature = - Bfield * Const::speedOfLight * particle->getCharge(); //Curvature of the track
407  double T = TMath::Sqrt(mom.Perp2() - 2 * curvature * r.Cross(mom) + curvature.Mag2() * r.Perp2());
408 
409  return TMath::Abs((-2 * r.Cross(mom).z() + curvature.Mag() * r.Perp2()) / (T + mom.Perp()));
410  }
411 
412  double ArmenterosLongitudinalMomentumAsymmetry(const Particle* part)
413  {
414  const auto& frame = ReferenceFrame::GetCurrent();
415  int nDaughters = part -> getNDaughters();
416  if (nDaughters != 2)
417  B2FATAL("You are trying to use an Armenteros variable. The mother particle is required to have exactly two daughters");
418 
419  const auto& daughters = part -> getDaughters();
420  TVector3 motherMomentum = frame.getMomentum(part).Vect();
421  TVector3 daughter1Momentum = frame.getMomentum(daughters[0]).Vect();
422  TVector3 daughter2Momentum = frame.getMomentum(daughters[1]).Vect();
423 
424  int daughter1Charge = daughters[0] -> getCharge();
425  int daughter2Charge = daughters[1] -> getCharge();
426  double daughter1Ql = daughter1Momentum.Dot(motherMomentum) / motherMomentum.Mag();
427  double daughter2Ql = daughter2Momentum.Dot(motherMomentum) / motherMomentum.Mag();
428 
429  double Arm_alpha;
430  if (daughter2Charge > daughter1Charge)
431  Arm_alpha = (daughter2Ql - daughter1Ql) / (daughter2Ql + daughter1Ql);
432  else
433  Arm_alpha = (daughter1Ql - daughter2Ql) / (daughter1Ql + daughter2Ql);
434 
435  return Arm_alpha;
436  }
437 
438  double ArmenterosDaughter1Qt(const Particle* part)
439  {
440  const auto& frame = ReferenceFrame::GetCurrent();
441  int nDaughters = part -> getNDaughters();
442  if (nDaughters != 2)
443  B2FATAL("You are trying to use an Armenteros variable. The mother particle is required to have exactly two daughters.");
444 
445  const auto& daughters = part -> getDaughters();
446  TVector3 motherMomentum = frame.getMomentum(part).Vect();
447  TVector3 daughter1Momentum = frame.getMomentum(daughters[0]).Vect();
448  double qt = daughter1Momentum.Perp(motherMomentum);
449 
450  return qt;
451  }
452 
453  double ArmenterosDaughter2Qt(const Particle* part)
454  {
455  const auto& frame = ReferenceFrame::GetCurrent();
456  int nDaughters = part -> getNDaughters();
457  if (nDaughters != 2)
458  B2FATAL("You are trying to use an Armenteros variable. The mother particle is required to have exactly two daughters.");
459 
460  const auto& daughters = part -> getDaughters();
461  TVector3 motherMomentum = frame.getMomentum(part).Vect();
462  TVector3 daughter2Momentum = frame.getMomentum(daughters[1]).Vect();
463  double qt = daughter2Momentum.Perp(motherMomentum);
464 
465  return qt;
466  }
467 
468 
469 // mass ------------------------------------------------------------
470 
471  double particleMass(const Particle* part)
472  {
473  return part->getMass();
474  }
475 
476  double particleDMass(const Particle* part)
477  {
478  return part->getMass() - part->getPDGMass();
479  }
480 
481  double particleInvariantMassFromDaughters(const Particle* part)
482  {
483  const std::vector<Particle*> daughters = part->getDaughters();
484  if (daughters.size() > 0) {
485  TLorentzVector sum;
486  for (auto daughter : daughters)
487  sum += daughter->get4Vector();
488 
489  return sum.M();
490  } else {
491  return part->getMass(); // !
492  }
493  }
494 
495  double particleInvariantMassLambda(const Particle* part)
496  {
497  const std::vector<Particle*> daughters = part->getDaughters();
498  if (daughters.size() == 2) {
499  TLorentzVector dt1;
500  TLorentzVector dt2;
501  TLorentzVector dtsum;
502  double mpi = Const::pionMass;
503  double mpr = Const::protonMass;
504  dt1 = daughters[0]->get4Vector();
505  dt2 = daughters[1]->get4Vector();
506  double E1 = hypot(mpi, dt1.P());
507  double E2 = hypot(mpr, dt2.P());
508  dtsum = dt1 + dt2;
509  return sqrt((E1 + E2) * (E1 + E2) - dtsum.P() * dtsum.P());
510 
511  } else {
512  return part->getMass();
513  }
514  }
515 
516  double particleInvariantMassError(const Particle* part)
517  {
518  float invMass = part->getMass();
519  TMatrixFSym covarianceMatrix = part->getMomentumErrorMatrix();
520 
521  TVectorF jacobian(Particle::c_DimMomentum);
522  jacobian[0] = -1.0 * part->getPx() / invMass;
523  jacobian[1] = -1.0 * part->getPy() / invMass;
524  jacobian[2] = -1.0 * part->getPz() / invMass;
525  jacobian[3] = 1.0 * part->getEnergy() / invMass;
526 
527  double result = jacobian * (covarianceMatrix * jacobian);
528 
529  if (result < 0.0)
530  return std::numeric_limits<double>::quiet_NaN();
531 
532  return TMath::Sqrt(result);
533  }
534 
535  double particleInvariantMassSignificance(const Particle* part)
536  {
537  return particleDMass(part) / particleInvariantMassError(part);
538  }
539 
540  double particleMassSquared(const Particle* part)
541  {
542  TLorentzVector p4 = part->get4Vector();
543  return p4.M2();
544  }
545 
546  double b2bTheta(const Particle* part)
547  {
548  PCmsLabTransform T;
549  TLorentzVector pcms = T.rotateLabToCms() * part->get4Vector();
550  TLorentzVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
551  TLorentzVector b2blab = T.rotateCmsToLab() * b2bcms;
552  return b2blab.Theta();
553  }
554 
555  double b2bPhi(const Particle* part)
556  {
557  PCmsLabTransform T;
558  TLorentzVector pcms = T.rotateLabToCms() * part->get4Vector();
559  TLorentzVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
560  TLorentzVector b2blab = T.rotateCmsToLab() * b2bcms;
561  return b2blab.Phi();
562  }
563 
564  double b2bClusterTheta(const Particle* part)
565  {
566  // get associated ECLCluster
567  const ECLCluster* cluster = part->getECLCluster();
568  if (!cluster) return std::numeric_limits<float>::quiet_NaN();
569  const ECLCluster::EHypothesisBit clusterHypothesis = part->getECLClusterEHypothesisBit();
570 
571  // get 4 momentum from cluster
572  ClusterUtils clutls;
573  TLorentzVector p4Cluster = clutls.Get4MomentumFromCluster(cluster, clusterHypothesis);
574 
575  // find the vector that balances this in the CMS
576  PCmsLabTransform T;
577  TLorentzVector pcms = T.rotateLabToCms() * p4Cluster;
578  TLorentzVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
579  TLorentzVector b2blab = T.rotateCmsToLab() * b2bcms;
580  return b2blab.Theta();
581  }
582 
583  double b2bClusterPhi(const Particle* part)
584  {
585  // get associated ECLCluster
586  const ECLCluster* cluster = part->getECLCluster();
587  if (!cluster) return std::numeric_limits<float>::quiet_NaN();
588  const ECLCluster::EHypothesisBit clusterHypothesis = part->getECLClusterEHypothesisBit();
589 
590  // get 4 momentum from cluster
591  ClusterUtils clutls;
592  TLorentzVector p4Cluster = clutls.Get4MomentumFromCluster(cluster, clusterHypothesis);
593 
594  // find the vector that balances this in the CMS
595  PCmsLabTransform T;
596  TLorentzVector pcms = T.rotateLabToCms() * p4Cluster;
597  TLorentzVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
598  TLorentzVector b2blab = T.rotateCmsToLab() * b2bcms;
599  return b2blab.Phi();
600  }
601 
602 
603 // released energy --------------------------------------------------
604 
605  double particleQ(const Particle* part)
606  {
607  float m = part->getMass();
608  for (unsigned i = 0; i < part->getNDaughters(); i++) {
609  const Particle* child = part->getDaughter(i);
610  if (child)
611  m -= child->getMass();
612  }
613  return m;
614  }
615 
616  double particleDQ(const Particle* part)
617  {
618  float m = part->getMass() - part->getPDGMass();
619  for (unsigned i = 0; i < part->getNDaughters(); i++) {
620  const Particle* child = part->getDaughter(i);
621  if (child)
622  m -= (child->getMass() - child->getPDGMass());
623  }
624  return m;
625  }
626 
627 // Mbc and deltaE
628 
629  double particleMbc(const Particle* part)
630  {
631  PCmsLabTransform T;
632  TLorentzVector vec = T.rotateLabToCms() * part->get4Vector();
633  double E = T.getCMSEnergy() / 2;
634  double m2 = E * E - vec.Vect().Mag2();
635  double mbc = m2 >= 0 ? sqrt(m2) : std::numeric_limits<double>::quiet_NaN();
636  return mbc;
637  }
638 
639  double particleDeltaE(const Particle* part)
640  {
641  PCmsLabTransform T;
642  TLorentzVector vec = T.rotateLabToCms() * part->get4Vector();
643  return vec.E() - T.getCMSEnergy() / 2;
644  }
645 
646 // other ------------------------------------------------------------
647 
648 
649  void printParticleInternal(const Particle* p, int depth)
650  {
651  stringstream s("");
652  for (int i = 0; i < depth; i++) {
653  s << " ";
654  }
655  s << p->getPDGCode();
656  const MCParticle* mcp = p->getMCParticle();
657  if (mcp) {
658  unsigned int flags = MCMatching::getMCErrors(p, mcp);
659  s << " -> MC: " << mcp->getPDG() << ", mcErrors: " << flags << " ("
660  << MCMatching::explainFlags(flags) << ")";
661  s << ", mc-index " << mcp->getIndex();
662  } else {
663  s << " (no MC match)";
664  }
665  s << ", mdst-source " << p->getMdstSource();
666  B2INFO(s.str());
667  for (const auto* daughter : p->getDaughters()) {
668  printParticleInternal(daughter, depth + 1);
669  }
670  }
671 
672  double printParticle(const Particle* p)
673  {
674  printParticleInternal(p, 0);
675  return 0.0;
676  }
677 
678 
679  double particleMCMomentumTransfer2(const Particle* part)
680  {
681  // for B meson MC particles only
682  const MCParticle* mcB = part->getMCParticle();
683 
684  if (!mcB)
685  return std::numeric_limits<double>::quiet_NaN();
686 
687  TLorentzVector pB = mcB->get4Vector();
688 
689  std::vector<MCParticle*> mcDaug = mcB->getDaughters();
690 
691  if (mcDaug.empty())
692  return std::numeric_limits<double>::quiet_NaN();
693 
694  // B -> X l nu
695  // q = pB - pX
696  TLorentzVector pX;
697 
698  for (auto mcTemp : mcDaug) {
699  if (abs(mcTemp->getPDG()) <= 16)
700  continue;
701 
702  pX += mcTemp->get4Vector();
703  }
704 
705  TLorentzVector q = pB - pX;
706 
707  return q.Mag2();
708  }
709 
710 // Recoil Kinematics related ---------------------------------------------
711  double recoilPx(const Particle* particle)
712  {
713  PCmsLabTransform T;
714 
715  // Initial state (e+e- momentum in LAB)
716  TLorentzVector pIN = T.getBeamFourMomentum();
717 
718  // Use requested frame for final calculation
719  const auto& frame = ReferenceFrame::GetCurrent();
720  return frame.getMomentum(pIN - particle->get4Vector()).Px();
721  }
722 
723  double recoilPy(const Particle* particle)
724  {
725  PCmsLabTransform T;
726 
727  // Initial state (e+e- momentum in LAB)
728  TLorentzVector pIN = T.getBeamFourMomentum();
729 
730  // Use requested frame for final calculation
731  const auto& frame = ReferenceFrame::GetCurrent();
732  return frame.getMomentum(pIN - particle->get4Vector()).Py();
733  }
734 
735  double recoilPz(const Particle* particle)
736  {
737  PCmsLabTransform T;
738 
739  // Initial state (e+e- momentum in LAB)
740  TLorentzVector pIN = T.getBeamFourMomentum();
741 
742  // Use requested frame for final calculation
743  const auto& frame = ReferenceFrame::GetCurrent();
744  return frame.getMomentum(pIN - particle->get4Vector()).Pz();
745  }
746 
747  double recoilMomentum(const Particle* particle)
748  {
749  PCmsLabTransform T;
750 
751  // Initial state (e+e- momentum in LAB)
752  TLorentzVector pIN = T.getBeamFourMomentum();
753 
754  // Use requested frame for final calculation
755  const auto& frame = ReferenceFrame::GetCurrent();
756  return frame.getMomentum(pIN - particle->get4Vector()).P();
757  }
758 
759  double recoilMomentumTheta(const Particle* particle)
760  {
761  PCmsLabTransform T;
762 
763  // Initial state (e+e- momentum in LAB)
764  TLorentzVector pIN = T.getBeamFourMomentum();
765 
766  // Use requested frame for final calculation
767  const auto& frame = ReferenceFrame::GetCurrent();
768  return frame.getMomentum(pIN - particle->get4Vector()).Theta();
769  }
770 
771  double recoilMomentumPhi(const Particle* particle)
772  {
773  PCmsLabTransform T;
774 
775  // Initial state (e+e- momentum in LAB)
776  TLorentzVector pIN = T.getBeamFourMomentum();
777 
778  // Use requested frame for final calculation
779  const auto& frame = ReferenceFrame::GetCurrent();
780  return frame.getMomentum(pIN - particle->get4Vector()).Phi();
781  }
782 
783  double recoilEnergy(const Particle* particle)
784  {
785  PCmsLabTransform T;
786 
787  // Initial state (e+e- momentum in LAB)
788  TLorentzVector pIN = T.getBeamFourMomentum();
789 
790  // Use requested frame for final calculation
791  const auto& frame = ReferenceFrame::GetCurrent();
792  return frame.getMomentum(pIN - particle->get4Vector()).E();
793  }
794 
795  double recoilMass(const Particle* particle)
796  {
797  PCmsLabTransform T;
798 
799  // Initial state (e+e- momentum in LAB)
800  TLorentzVector pIN = T.getBeamFourMomentum();
801 
802  // Use requested frame for final calculation
803  const auto& frame = ReferenceFrame::GetCurrent();
804  return frame.getMomentum(pIN - particle->get4Vector()).M();
805  }
806 
807  double recoilMassSquared(const Particle* particle)
808  {
809  PCmsLabTransform T;
810 
811  // Initial state (e+e- momentum in LAB)
812  TLorentzVector pIN = T.getBeamFourMomentum();
813 
814  // Use requested frame for final calculation
815  const auto& frame = ReferenceFrame::GetCurrent();
816  return frame.getMomentum(pIN - particle->get4Vector()).M2();
817  }
818 
819  double m2RecoilSignalSide(const Particle* part)
820  {
821  PCmsLabTransform T;
822  double beamEnergy = T.getCMSEnergy() / 2.;
823  if (part->getNDaughters() != 2) return std::numeric_limits<double>::quiet_NaN();
824  TLorentzVector tagVec = T.rotateLabToCms()
825  * part->getDaughter(0)->get4Vector();
826  TLorentzVector sigVec = T.rotateLabToCms()
827  * part->getDaughter(1)->get4Vector();
828  tagVec.SetE(-beamEnergy);
829  return (-tagVec - sigVec).M2();
830  }
831 
832  double recoilMCDecayType(const Particle* particle)
833  {
834  auto* mcp = particle->getMCParticle();
835 
836  if (!mcp)
837  return std::numeric_limits<double>::quiet_NaN();
838 
839  MCParticle* mcMother = mcp->getMother();
840 
841  if (!mcMother)
842  return std::numeric_limits<double>::quiet_NaN();
843 
844  std::vector<MCParticle*> daughters = mcMother->getDaughters();
845 
846  if (daughters.size() != 2)
847  return std::numeric_limits<double>::quiet_NaN();
848 
849  MCParticle* recoilMC = nullptr;
850  if (daughters[0]->getArrayIndex() == mcp->getArrayIndex())
851  recoilMC = daughters[1];
852  else
853  recoilMC = daughters[0];
854 
855  if (!recoilMC->hasStatus(MCParticle::c_PrimaryParticle))
856  return std::numeric_limits<double>::quiet_NaN();
857 
858  int decayType = 0;
859  checkMCParticleDecay(recoilMC, decayType, false);
860 
861  if (decayType == 0)
862  checkMCParticleDecay(recoilMC, decayType, true);
863 
864  return decayType;
865  }
866 
867  void checkMCParticleDecay(MCParticle* mcp, int& decayType, bool recursive)
868  {
869  int nHadronicParticles = 0;
870  int nPrimaryParticleDaughters = 0;
871  std::vector<MCParticle*> daughters = mcp->getDaughters();
872 
873  // Are any of the daughters primary particles? How many of them are hadrons?
874  for (auto& daughter : daughters) {
875  if (!daughter->hasStatus(MCParticle::c_PrimaryParticle))
876  continue;
877 
878  nPrimaryParticleDaughters++;
879  if (abs(daughter->getPDG()) > Const::photon.getPDGCode())
880  nHadronicParticles++;
881  }
882 
883  if (nPrimaryParticleDaughters > 1) {
884  for (auto& daughter : daughters) {
885  if (!daughter->hasStatus(MCParticle::c_PrimaryParticle))
886  continue;
887 
888  if (abs(daughter->getPDG()) == 12 or abs(daughter->getPDG()) == 14 or abs(daughter->getPDG()) == 16) {
889  if (!recursive) {
890  if (nHadronicParticles == 0) {
891  decayType = 1.0;
892  break;
893  } else {
894  decayType = 2.0;
895  break;
896  }
897  } else {
898  decayType = 3.0;
899  break;
900  }
901  }
902 
903  else if (recursive)
904  checkMCParticleDecay(daughter, decayType, recursive);
905  }
906  }
907  }
908 
909  double nRemainingTracksInEvent(const Particle* particle)
910  {
911 
912  StoreArray<Track> tracks;
913  int event_tracks = tracks.getEntries();
914 
915  int par_tracks = 0;
916  const auto& daughters = particle->getFinalStateDaughters();
917  for (const auto& daughter : daughters) {
918  int pdg = abs(daughter->getPDGCode());
919  if (pdg == Const::electron.getPDGCode() or pdg == Const::muon.getPDGCode() or pdg == Const::pion.getPDGCode()
920  or pdg == Const::kaon.getPDGCode() or pdg == Const::proton.getPDGCode())
921  par_tracks++;
922  }
923  return event_tracks - par_tracks;
924  }
925 
926  double trackMatchType(const Particle* particle)
927  {
928  // Particle does not contain a ECL Cluster
929  double result = std::numeric_limits<double>::quiet_NaN();
930 
931  const ECLCluster* cluster = particle->getECLCluster();
932  if (cluster) {
933  // No associated track is default
934  result = 0;
935  if (cluster->isTrack()) {
936  // There is a track match
937  result = 1.0;
938  }
939  }
940  return result;
941  }
942 
943  double False(const Particle*)
944  {
945  return 0;
946  }
947 
948  double True(const Particle*)
949  {
950  return 1;
951  }
952 
953  double infinity(const Particle*)
954  {
955  double inf = std::numeric_limits<double>::infinity();
956  return inf;
957  }
958 
959  double random(const Particle*)
960  {
961  return gRandom->Uniform(0, 1);
962  }
963 
964  double eventRandom(const Particle*)
965  {
966  std::string key = "__eventRandom";
967  StoreObjPtr<EventExtraInfo> eventExtraInfo;
968  if (not eventExtraInfo.isValid())
969  eventExtraInfo.create();
970  if (eventExtraInfo->hasExtraInfo(key)) {
971  return eventExtraInfo->getExtraInfo(key);
972  } else {
973  double value = gRandom->Uniform(0, 1);
974  eventExtraInfo->addExtraInfo(key, value);
975  return value;
976  }
977  }
978 
979  VARIABLE_GROUP("Kinematics");
980  REGISTER_VARIABLE("p", particleP, "momentum magnitude");
981  REGISTER_VARIABLE("E", particleE, "energy");
982 
983  REGISTER_VARIABLE("E_uncertainty", particleEUncertainty, "energy uncertainty (sqrt(sigma2))");
984  REGISTER_VARIABLE("ECLClusterE_uncertainty", particleClusterEUncertainty,
985  "energy uncertainty as given by the underlying ECL cluster.");
986  REGISTER_VARIABLE("px", particlePx, "momentum component x");
987  REGISTER_VARIABLE("py", particlePy, "momentum component y");
988  REGISTER_VARIABLE("pz", particlePz, "momentum component z");
989  REGISTER_VARIABLE("pt", particlePt, "transverse momentum");
990  REGISTER_VARIABLE("xp", particleXp,
991  "scaled momentum: the momentum of the particle in the CMS as a fraction of its maximum available momentum in the collision");
992  REGISTER_VARIABLE("pErr", particlePErr, "error of momentum magnitude");
993  REGISTER_VARIABLE("pxErr", particlePxErr, "error of momentum component x");
994  REGISTER_VARIABLE("pyErr", particlePyErr, "error of momentum component y");
995  REGISTER_VARIABLE("pzErr", particlePzErr, "error of momentum component z");
996  REGISTER_VARIABLE("ptErr", particlePtErr, "error of transverse momentum");
997  REGISTER_VARIABLE("momVertCovM(i,j)", covMatrixElement,
998  "returns the (i,j)-th element of the MomentumVertex Covariance Matrix (7x7).\n"
999  "Order of elements in the covariance matrix is: px, py, pz, E, x, y, z.");
1000  REGISTER_VARIABLE("momDevChi2", momentumDeviationChi2,
1001  "momentum deviation chi^2 value calculated as"
1002  "chi^2 = sum_i (p_i - mc(p_i))^2/sigma(p_i)^2, where sum runs over i = px, py, pz and"
1003  "mc(p_i) is the mc truth value and sigma(p_i) is the estimated error of i-th component of momentum vector")
1004  REGISTER_VARIABLE("theta", particleTheta, "polar angle in radians");
1005  REGISTER_VARIABLE("thetaErr", particleThetaErr, "error of polar angle in radians");
1006  REGISTER_VARIABLE("cosTheta", particleCosTheta, "momentum cosine of polar angle");
1007  REGISTER_VARIABLE("cosThetaErr", particleCosThetaErr, "error of momentum cosine of polar angle");
1008  REGISTER_VARIABLE("phi", particlePhi, "momentum azimuthal angle in radians");
1009  REGISTER_VARIABLE("phiErr", particlePhiErr, "error of momentum azimuthal angle in radians");
1010  REGISTER_VARIABLE("PDG", particlePDGCode, "PDG code");
1011 
1012  REGISTER_VARIABLE("cosAngleBetweenMomentumAndVertexVectorInXYPlane",
1013  cosAngleBetweenMomentumAndVertexVectorInXYPlane,
1014  "cosine of the angle between momentum and vertex vector (vector connecting ip and fitted vertex) of this particle in xy-plane");
1015  REGISTER_VARIABLE("cosAngleBetweenMomentumAndVertexVector",
1016  cosAngleBetweenMomentumAndVertexVector,
1017  "cosine of the angle between momentum and vertex vector (vector connecting ip and fitted vertex) of this particle");
1018  REGISTER_VARIABLE("cosThetaBetweenParticleAndNominalB",
1019  cosThetaBetweenParticleAndNominalB,
1020  "cosine of the angle in CMS between momentum the particle and a nominal B particle. It is somewhere between -1 and 1 if only a massless particle like a neutrino is missing in the reconstruction.");
1021  REGISTER_VARIABLE("cosToThrustOfEvent", cosToThrustOfEvent,
1022  "Returns the cosine of the angle between the particle and the thrust axis of the event, as calculate by the EventShapeCalculator module. buildEventShape() must be run before calling this variable")
1023 
1024  REGISTER_VARIABLE("ImpactXY" , ImpactXY , "The impact parameter of the given particle in the xy plane");
1025 
1026  REGISTER_VARIABLE("M", particleMass, R"DOC(
1027 The particle's mass.
1028 
1029 Note that this is context-dependent variable and can take different values depending on the situation. This should be the "best" value possible with the information provided.
1030 
1031 - If this particle is track- or cluster- based, then this is the value of the mass hypothesis.
1032 - If this particle is an MC particle then this is the mass of that particle.
1033 - If this particle is composite, then *initially* this takes the value of the invariant mass of the daughters.
1034 - If this particle is composite and a *mass or vertex fit* has been performed then this may be updated by the fit.
1035 
1036  * You will see a difference between this mass and the :b2:var:`InvM`.
1037  )DOC");
1038  REGISTER_VARIABLE("dM", particleDMass, "mass minus nominal mass");
1039  REGISTER_VARIABLE("Q", particleQ, "released energy in decay");
1040  REGISTER_VARIABLE("dQ", particleDQ,
1041  "released energy in decay minus nominal one");
1042  REGISTER_VARIABLE("Mbc", particleMbc, "beam constrained mass");
1043  REGISTER_VARIABLE("deltaE", particleDeltaE, "energy difference");
1044  REGISTER_VARIABLE("M2", particleMassSquared,
1045  "The particle's mass squared.");
1046 
1047  REGISTER_VARIABLE("InvM", particleInvariantMassFromDaughters,
1048  "invariant mass (determined from particle's daughter 4-momentum vectors). If this particle has no daughters, defaults to :b2:var:`M`.");
1049  REGISTER_VARIABLE("InvMLambda", particleInvariantMassLambda,
1050  "Invariant mass (determined from particle's daughter 4-momentum vectors), assuming the first daughter is a pion and the second daughter is a proton.\n"
1051  "If the particle has not 2 daughters, it returns just the mass value.");
1052 
1053  REGISTER_VARIABLE("ErrM", particleInvariantMassError,
1054  "uncertainty of invariant mass");
1055  REGISTER_VARIABLE("SigM", particleInvariantMassSignificance,
1056  "signed deviation of particle's invariant mass from its nominal mass in units of the uncertainty on the invariant mass (dM/ErrM)");
1057 
1058  REGISTER_VARIABLE("pxRecoil", recoilPx,
1059  "component x of 3-momentum recoiling against given Particle");
1060  REGISTER_VARIABLE("pyRecoil", recoilPy,
1061  "component y of 3-momentum recoiling against given Particle");
1062  REGISTER_VARIABLE("pzRecoil", recoilPz,
1063  "component z of 3-momentum recoiling against given Particle");
1064 
1065  REGISTER_VARIABLE("pRecoil", recoilMomentum,
1066  "magnitude of 3 - momentum recoiling against given Particle");
1067  REGISTER_VARIABLE("pRecoilTheta", recoilMomentumTheta,
1068  "Polar angle of a particle's missing momentum");
1069  REGISTER_VARIABLE("pRecoilPhi", recoilMomentumPhi,
1070  "Azimuthal angle of a particle's missing momentum");
1071  REGISTER_VARIABLE("eRecoil", recoilEnergy,
1072  "energy recoiling against given Particle");
1073  REGISTER_VARIABLE("mRecoil", recoilMass,
1074  "Invariant mass of the system recoiling against given Particle");
1075  REGISTER_VARIABLE("m2Recoil", recoilMassSquared,
1076  "invariant mass squared of the system recoiling against given Particle");
1077  REGISTER_VARIABLE("m2RecoilSignalSide", m2RecoilSignalSide,
1078  "Squared recoil mass of the signal side which is calculated in the CMS frame under the assumption that the signal and tag side are produced back to back and the tag side energy equals the beam energy. The variable must be applied to the Upsilon and the tag side must be the first, the signal side the second daughter ");
1079 
1080  REGISTER_VARIABLE("b2bTheta", b2bTheta,
1081  "Polar angle in the lab system that is back-to-back to the particle in the CMS. Useful for low multiplicity studies.")
1082  REGISTER_VARIABLE("b2bPhi", b2bPhi,
1083  "Azimuthal angle in the lab system that is back-to-back to the particle in the CMS. Useful for low multiplicity studies.")
1084  REGISTER_VARIABLE("b2bClusterTheta", b2bClusterTheta,
1085  "Polar angle in the lab system that is back-to-back to the particle's associated ECLCluster in the CMS. Returns NAN if no cluster is found. Useful for low multiplicity studies.")
1086  REGISTER_VARIABLE("b2bClusterPhi", b2bClusterPhi,
1087  "Azimuthal angle in the lab system that is back-to-back to the particle's associated ECLCluster in the CMS. Returns NAN if no cluster is found. Useful for low multiplicity studies.")
1088  REGISTER_VARIABLE("ArmenterosLongitudinalMomentumAsymmetry", ArmenterosLongitudinalMomentumAsymmetry,
1089  "Longitudinal momentum asymmetry of V0's daughters.\n"
1090  "The mother (V0) is required to have exactly two daughters");
1091  REGISTER_VARIABLE("ArmenterosDaughter1Qt", ArmenterosDaughter1Qt,
1092  "Transverse momentum of the first daughter with respect to the V0 mother.\n"
1093  "The mother is required to have exactly two daughters");
1094  REGISTER_VARIABLE("ArmenterosDaughter2Qt", ArmenterosDaughter2Qt,
1095  "Transverse momentum of the second daughter with respect to the V0 mother.\n"
1096  "The mother is required to have exactly two daughters");
1097 
1098  VARIABLE_GROUP("Miscellaneous");
1099  REGISTER_VARIABLE("nRemainingTracksInEvent", nRemainingTracksInEvent,
1100  "Number of tracks in the event - Number of tracks( = charged FSPs) of particle.");
1101  REGISTER_VARIABLE("trackMatchType", trackMatchType, R"DOC(
1102 
1103  * -1 particle has no ECL cluster
1104  * 0 particle has no associated track
1105  * 1 there is a matched track called connected - region(CR) track match
1106  )DOC");
1107  MAKE_DEPRECATED("trackMatchType", false, "light-minos-2012", R"DOC(
1108  Use better variables like `trackNECLClusters`, `clusterTrackMatch`, and `nECLClusterTrackMatches`.)DOC");
1109 
1110  REGISTER_VARIABLE("decayTypeRecoil", recoilMCDecayType,
1111  "type of the particle decay(no related mcparticle = -1, hadronic = 0, direct leptonic = 1, direct semileptonic = 2,"
1112  "lower level leptonic = 3.");
1113 
1114  REGISTER_VARIABLE("printParticle", printParticle,
1115  "For debugging, print Particle and daughter PDG codes, plus MC match. Returns 0.");
1116  REGISTER_VARIABLE("mcMomTransfer2", particleMCMomentumTransfer2,
1117  "Return the true momentum transfer to lepton pair in a B(semi -) leptonic B meson decay.");
1118  REGISTER_VARIABLE("False", False,
1119  "returns always 0, used for testing and debugging.");
1120  REGISTER_VARIABLE("True", True,
1121  "returns always 1, used for testing and debugging.");
1122  REGISTER_VARIABLE("infinity", infinity,
1123  "returns std::numeric_limits<double>::infinity()");
1124  REGISTER_VARIABLE("random", random,
1125  "return a random number between 0 and 1 for each candidate. Can be used, e.g. for picking a random"
1126  "candidate in the best candidate selection.");
1127  REGISTER_VARIABLE("eventRandom", eventRandom,
1128  "[Eventbased] Returns a random number between 0 and 1 for this event. Can be used, e.g. for applying an event prescale.");
1129 
1130  }
1132 }
MAKE_DEPRECATED
#define MAKE_DEPRECATED(name, make_fatal, version, description)
Registers a variable as deprecated.
Definition: Manager.h:345
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17