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