Belle II Software  release-08-01-10
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 header.
10 #include <analysis/variables/Variables.h>
11 
12 // include VariableManager
13 #include <analysis/VariableManager/Manager.h>
14 
15 #include <analysis/utility/PCmsLabTransform.h>
16 #include <analysis/utility/ReferenceFrame.h>
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/ParticleList.h>
27 #include <framework/dataobjects/EventExtraInfo.h>
28 #include <analysis/dataobjects/EventShapeContainer.h>
29 
30 #include <mdst/dataobjects/MCParticle.h>
31 #include <mdst/dataobjects/Track.h>
32 #include <mdst/dataobjects/ECLCluster.h>
33 #include <mdst/dataobjects/V0.h>
34 
35 #include <mdst/dbobjects/BeamSpot.h>
36 
37 // framework aux
38 #include <framework/logging/Logger.h>
39 #include <framework/geometry/B2Vector3.h>
40 #include <framework/geometry/BFieldManager.h>
41 #include <framework/gearbox/Const.h>
42 #include <framework/utilities/Conversion.h>
43 
44 #include <Math/Vector4D.h>
45 #include <TRandom.h>
46 #include <TVectorF.h>
47 
48 #include <iostream>
49 #include <cmath>
50 #include <boost/algorithm/string.hpp>
51 
52 
53 using namespace std;
54 
55 namespace Belle2 {
60  namespace Variable {
61 
62 // momentum -------------------------------------------
63 
64  double particleP(const Particle* part)
65  {
66  const auto& frame = ReferenceFrame::GetCurrent();
67  return frame.getMomentum(part).P();
68  }
69 
70  double particleE(const Particle* part)
71  {
72  const auto& frame = ReferenceFrame::GetCurrent();
73  return frame.getMomentum(part).E();
74  }
75 
76  double particleClusterEUncertainty(const Particle* part)
77  {
78  const ECLCluster* cluster = part->getECLCluster();
79  if (cluster) {
80  const auto EPhiThetaCov = cluster->getCovarianceMatrix3x3();
81  return std::sqrt(EPhiThetaCov[0][0]);
82  }
83  return Const::doubleNaN;
84  }
85 
86  double particlePx(const Particle* part)
87  {
88  const auto& frame = ReferenceFrame::GetCurrent();
89  return frame.getMomentum(part).Px();
90  }
91 
92  double particlePy(const Particle* part)
93  {
94  const auto& frame = ReferenceFrame::GetCurrent();
95  return frame.getMomentum(part).Py();
96  }
97 
98  double particlePz(const Particle* part)
99  {
100  const auto& frame = ReferenceFrame::GetCurrent();
101  return frame.getMomentum(part).Pz();
102  }
103 
104  double particlePt(const Particle* part)
105  {
106  const auto& frame = ReferenceFrame::GetCurrent();
107  return frame.getMomentum(part).Pt();
108  }
109 
110  double covMatrixElement(const Particle* part, const std::vector<double>& element)
111  {
112  int elementI = int(std::lround(element[0]));
113  int elementJ = int(std::lround(element[1]));
114 
115  if (elementI < 0 || elementI > 6) {
116  B2WARNING("Requested particle's momentumVertex covariance matrix element is out of boundaries [0 - 6]:" << LogVar("i", elementI));
117  return Const::doubleNaN;
118  }
119  if (elementJ < 0 || elementJ > 6) {
120  B2WARNING("Requested particle's momentumVertex covariance matrix element is out of boundaries [0 - 6]:" << LogVar("j", elementJ));
121  return Const::doubleNaN;
122  }
123 
124  return part->getMomentumVertexErrorMatrix()(elementI, elementJ);
125  }
126 
127  double particleEUncertainty(const Particle* part)
128  {
129  const auto& frame = ReferenceFrame::GetCurrent();
130 
131  double errorSquared = frame.getMomentumErrorMatrix(part)(3, 3);
132 
133  if (errorSquared >= 0.0)
134  return sqrt(errorSquared);
135  else
136  return Const::doubleNaN;
137  }
138 
139  double particlePErr(const Particle* part)
140  {
141  TMatrixD jacobianRot(3, 3);
142  jacobianRot.Zero();
143 
144  double cosPhi = cos(particlePhi(part));
145  double sinPhi = sin(particlePhi(part));
146  double cosTheta = particleCosTheta(part);
147  double sinTheta = sin(acos(cosTheta));
148  double p = particleP(part);
149 
150  jacobianRot(0, 0) = sinTheta * cosPhi;
151  jacobianRot(0, 1) = sinTheta * sinPhi;
152  jacobianRot(1, 0) = cosTheta * cosPhi / p;
153  jacobianRot(1, 1) = cosTheta * sinPhi / p;
154  jacobianRot(0, 2) = cosTheta;
155  jacobianRot(2, 0) = -sinPhi / sinTheta / p;
156  jacobianRot(1, 2) = -sinTheta / p;
157  jacobianRot(2, 1) = cosPhi / sinTheta / p;
158 
159  const auto& frame = ReferenceFrame::GetCurrent();
160 
161  double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2, " ").Similarity(jacobianRot)(0, 0);
162 
163  if (errorSquared >= 0.0)
164  return sqrt(errorSquared);
165  else
166  return Const::doubleNaN;
167  }
168 
169  double particlePxErr(const Particle* part)
170  {
171  const auto& frame = ReferenceFrame::GetCurrent();
172 
173  double errorSquared = frame.getMomentumErrorMatrix(part)(0, 0);
174 
175  if (errorSquared >= 0.0)
176  return sqrt(errorSquared);
177  else
178  return Const::doubleNaN;
179  }
180 
181  double particlePyErr(const Particle* part)
182  {
183  const auto& frame = ReferenceFrame::GetCurrent();
184  double errorSquared = frame.getMomentumErrorMatrix(part)(1, 1);
185 
186  if (errorSquared >= 0.0)
187  return sqrt(errorSquared);
188  else
189  return Const::doubleNaN;
190  }
191 
192  double particlePzErr(const Particle* part)
193  {
194  const auto& frame = ReferenceFrame::GetCurrent();
195  double errorSquared = frame.getMomentumErrorMatrix(part)(2, 2);
196 
197  if (errorSquared >= 0.0)
198  return sqrt(errorSquared);
199  else
200  return Const::doubleNaN;
201  }
202 
203  double particlePtErr(const Particle* part)
204  {
205  TMatrixD jacobianRot(3, 3);
206  jacobianRot.Zero();
207 
208  double px = particlePx(part);
209  double py = particlePy(part);
210  double pt = particlePt(part);
211 
212  jacobianRot(0, 0) = px / pt;
213  jacobianRot(0, 1) = py / pt;
214  jacobianRot(1, 0) = -py / (pt * pt);
215  jacobianRot(1, 1) = px / (pt * pt);
216  jacobianRot(2, 2) = 1;
217 
218  const auto& frame = ReferenceFrame::GetCurrent();
219  double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2, " ").Similarity(jacobianRot)(0, 0);
220 
221  if (errorSquared >= 0.0)
222  return sqrt(errorSquared);
223  else
224  return Const::doubleNaN;
225  }
226 
227  double momentumDeviationChi2(const Particle* part)
228  {
229  double result = Const::doubleNaN;
230 
231  // check if error matrix is set
232  if (part->getPValue() < 0.0)
233  return result;
234 
235  // check if mc match exists
236  const MCParticle* mcp = part->getMCParticle();
237  if (mcp == nullptr)
238  return result;
239 
240  result = 0.0;
241  result += TMath::Power(part->getPx() - mcp->getMomentum().X(), 2.0) / part->getMomentumVertexErrorMatrix()(0, 0);
242  result += TMath::Power(part->getPy() - mcp->getMomentum().Y(), 2.0) / part->getMomentumVertexErrorMatrix()(1, 1);
243  result += TMath::Power(part->getPz() - mcp->getMomentum().Z(), 2.0) / part->getMomentumVertexErrorMatrix()(2, 2);
244 
245  return result;
246  }
247 
248  double particleTheta(const Particle* part)
249  {
250  const auto& frame = ReferenceFrame::GetCurrent();
251  return frame.getMomentum(part).Theta();
252  }
253 
254  double particleThetaErr(const Particle* part)
255  {
256  TMatrixD jacobianRot(3, 3);
257  jacobianRot.Zero();
258 
259  double cosPhi = cos(particlePhi(part));
260  double sinPhi = sin(particlePhi(part));
261  double cosTheta = particleCosTheta(part);
262  double sinTheta = sin(acos(cosTheta));
263  double p = particleP(part);
264 
265  jacobianRot(0, 0) = sinTheta * cosPhi;
266  jacobianRot(0, 1) = sinTheta * sinPhi;
267  jacobianRot(1, 0) = cosTheta * cosPhi / p;
268  jacobianRot(1, 1) = cosTheta * sinPhi / p;
269  jacobianRot(0, 2) = cosTheta;
270  jacobianRot(2, 0) = -sinPhi / sinTheta / p;
271  jacobianRot(1, 2) = -sinTheta / p;
272  jacobianRot(2, 1) = cosPhi / sinTheta / p;
273 
274  const auto& frame = ReferenceFrame::GetCurrent();
275  double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2, " ").Similarity(jacobianRot)(1, 1);
276 
277  if (errorSquared >= 0.0)
278  return sqrt(errorSquared);
279  else
280  return Const::doubleNaN;
281  }
282 
283  double particleCosTheta(const Particle* part)
284  {
285  const auto& frame = ReferenceFrame::GetCurrent();
286  return cos(frame.getMomentum(part).Theta());
287  }
288 
289  double particleCosThetaErr(const Particle* part)
290  {
291  return fabs(particleThetaErr(part) * sin(particleTheta(part)));
292  }
293 
294  double particlePhi(const Particle* part)
295  {
296  const auto& frame = ReferenceFrame::GetCurrent();
297  return frame.getMomentum(part).Phi();
298  }
299 
300  double particlePhiErr(const Particle* part)
301  {
302  TMatrixD jacobianRot(3, 3);
303  jacobianRot.Zero();
304 
305  double px = particlePx(part);
306  double py = particlePy(part);
307  double pt = particlePt(part);
308 
309  jacobianRot(0, 0) = px / pt;
310  jacobianRot(0, 1) = py / pt;
311  jacobianRot(1, 0) = -py / (pt * pt);
312  jacobianRot(1, 1) = px / (pt * pt);
313  jacobianRot(2, 2) = 1;
314 
315  const auto& frame = ReferenceFrame::GetCurrent();
316  double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2, " ").Similarity(jacobianRot)(1, 1);
317 
318  if (errorSquared >= 0.0)
319  return sqrt(errorSquared);
320  else
321  return Const::doubleNaN;
322  }
323 
324  double particleXp(const Particle* part)
325  {
326  PCmsLabTransform T;
327  ROOT::Math::PxPyPzEVector p4 = part -> get4Vector();
328  ROOT::Math::PxPyPzEVector p4CMS = T.rotateLabToCms() * p4;
329  float s = T.getCMSEnergy();
330  float M = part->getMass();
331  return p4CMS.P() / TMath::Sqrt(s * s / 4 - M * M);
332  }
333 
334  int particlePDGCode(const Particle* part)
335  {
336  return part->getPDGCode();
337  }
338 
339  double cosAngleBetweenMomentumAndVertexVectorInXYPlane(const Particle* part)
340  {
341  static DBObjPtr<BeamSpot> beamSpotDB;
342  double px = part->getPx();
343  double py = part->getPy();
344 
345  double xV = part->getX();
346  double yV = part->getY();
347 
348  double xIP = (beamSpotDB->getIPPosition()).X();
349  double yIP = (beamSpotDB->getIPPosition()).Y();
350 
351  double x = xV - xIP;
352  double y = yV - yIP;
353 
354  double cosangle = (px * x + py * y) / (sqrt(px * px + py * py) * sqrt(x * x + y * y));
355  return cosangle;
356  }
357 
358  double cosAngleBetweenMomentumAndVertexVector(const Particle* part)
359  {
360  static DBObjPtr<BeamSpot> beamSpotDB;
361  return cos((B2Vector3D(part->getVertex()) - beamSpotDB->getIPPosition()).Angle(B2Vector3D(part->getMomentum())));
362  }
363 
364  double cosThetaBetweenParticleAndNominalB(const Particle* part)
365  {
366 
367  int particlePDG = abs(part->getPDGCode());
368  if (particlePDG != 511 and particlePDG != 521)
369  B2FATAL("The Variables cosThetaBetweenParticleAndNominalB is only meant to be used on B mesons!");
370 
371  PCmsLabTransform T;
372  double e_Beam = T.getCMSEnergy() / 2.0; // GeV
373  double m_B = part->getPDGMass();
374  // if this is a continuum run, use an approximate Y(4S) CMS energy
375  if (e_Beam * e_Beam - m_B * m_B < 0) {
376  e_Beam = 1.0579400E+1 / 2.0; // GeV
377  }
378  double p_B = std::sqrt(e_Beam * e_Beam - m_B * m_B);
379 
380  ROOT::Math::PxPyPzEVector p = T.rotateLabToCms() * part->get4Vector();
381  double e_d = p.E();
382  double m_d = p.M();
383  double p_d = p.P();
384 
385  double theta_BY = (2 * e_Beam * e_d - m_B * m_B - m_d * m_d)
386  / (2 * p_B * p_d);
387  return theta_BY;
388  }
389 
390  double cosToThrustOfEvent(const Particle* part)
391  {
392  StoreObjPtr<EventShapeContainer> evtShape;
393  if (!evtShape) {
394  B2WARNING("Cannot find thrust of event information, did you forget to load the event shape calculation?");
395  return Const::doubleNaN;
396  }
397  PCmsLabTransform T;
398  B2Vector3D th = evtShape->getThrustAxis();
399  B2Vector3D particleMomentum = (T.rotateLabToCms() * part -> get4Vector()).Vect();
400  return std::cos(th.Angle(particleMomentum));
401  }
402 
403  double ImpactXY(const Particle* particle)
404  {
405  static DBObjPtr<BeamSpot> beamSpotDB;
406 
407  ROOT::Math::XYZVector mom = particle->getMomentum();
408 
409  ROOT::Math::XYZVector r = particle->getVertex() - ROOT::Math::XYZVector(beamSpotDB->getIPPosition());
410 
411  ROOT::Math::XYZVector Bfield = BFieldManager::getInstance().getFieldInTesla(ROOT::Math::XYZVector(beamSpotDB->getIPPosition()));
412 
413  ROOT::Math::XYZVector curvature = - Bfield * Const::speedOfLight * particle->getCharge(); //Curvature of the track
414  double T = TMath::Sqrt(mom.Perp2() - 2.0 * curvature.Dot(r.Cross(mom)) + curvature.Mag2() * r.Perp2());
415 
416  return TMath::Abs((-2 * r.Cross(mom).Z() + curvature.R() * r.Perp2()) / (T + mom.Rho()));
417  }
418 
419  double ArmenterosLongitudinalMomentumAsymmetry(const Particle* part)
420  {
421  const auto& frame = ReferenceFrame::GetCurrent();
422  int nDaughters = part -> getNDaughters();
423  if (nDaughters != 2)
424  B2FATAL("You are trying to use an Armenteros variable. The mother particle is required to have exactly two daughters");
425 
426  const auto& daughters = part -> getDaughters();
427  B2Vector3D motherMomentum = frame.getMomentum(part).Vect();
428  B2Vector3D daughter1Momentum = frame.getMomentum(daughters[0]).Vect();
429  B2Vector3D daughter2Momentum = frame.getMomentum(daughters[1]).Vect();
430 
431  int daughter1Charge = daughters[0] -> getCharge();
432  int daughter2Charge = daughters[1] -> getCharge();
433  double daughter1Ql = daughter1Momentum.Dot(motherMomentum) / motherMomentum.Mag();
434  double daughter2Ql = daughter2Momentum.Dot(motherMomentum) / motherMomentum.Mag();
435 
436  double Arm_alpha;
437  if (daughter2Charge > daughter1Charge)
438  Arm_alpha = (daughter2Ql - daughter1Ql) / (daughter2Ql + daughter1Ql);
439  else
440  Arm_alpha = (daughter1Ql - daughter2Ql) / (daughter1Ql + daughter2Ql);
441 
442  return Arm_alpha;
443  }
444 
445  double ArmenterosDaughter1Qt(const Particle* part)
446  {
447  const auto& frame = ReferenceFrame::GetCurrent();
448  int nDaughters = part -> getNDaughters();
449  if (nDaughters != 2)
450  B2FATAL("You are trying to use an Armenteros variable. The mother particle is required to have exactly two daughters.");
451 
452  const auto& daughters = part -> getDaughters();
453  B2Vector3D motherMomentum = frame.getMomentum(part).Vect();
454  B2Vector3D daughter1Momentum = frame.getMomentum(daughters[0]).Vect();
455  double qt = daughter1Momentum.Perp(motherMomentum);
456 
457  return qt;
458  }
459 
460  double ArmenterosDaughter2Qt(const Particle* part)
461  {
462  const auto& frame = ReferenceFrame::GetCurrent();
463  int nDaughters = part -> getNDaughters();
464  if (nDaughters != 2)
465  B2FATAL("You are trying to use an Armenteros variable. The mother particle is required to have exactly two daughters.");
466 
467  const auto& daughters = part -> getDaughters();
468  B2Vector3D motherMomentum = frame.getMomentum(part).Vect();
469  B2Vector3D daughter2Momentum = frame.getMomentum(daughters[1]).Vect();
470  double qt = daughter2Momentum.Perp(motherMomentum);
471 
472  return qt;
473  }
474 
475 
476 // mass ------------------------------------------------------------
477 
478  double particleMass(const Particle* part)
479  {
480  return part->getMass();
481  }
482 
483  double particleDMass(const Particle* part)
484  {
485  return part->getMass() - part->getPDGMass();
486  }
487 
488  double particleInvariantMassFromDaughters(const Particle* part)
489  {
490  const std::vector<Particle*> daughters = part->getDaughters();
491  if (daughters.size() > 0) {
492  ROOT::Math::PxPyPzEVector sum;
493  for (auto daughter : daughters)
494  sum += daughter->get4Vector();
495 
496  return sum.M();
497  } else {
498  return part->getMass(); // !
499  }
500  }
501 
502  double particleInvariantMassFromDaughtersDisplaced(const Particle* part)
503  {
504  ROOT::Math::XYZVector vertex = part->getVertex();
505  if (part->getParticleSource() != Particle::EParticleSourceObject::c_V0
506  && vertex.Rho() < 0.5) return particleInvariantMassFromDaughters(part);
507 
508  const std::vector<Particle*> daughters = part->getDaughters();
509  if (daughters.size() == 0) return particleInvariantMassFromDaughters(part);
510 
511  const double bField = BFieldManager::getFieldInTesla(vertex).Z();
512  ROOT::Math::PxPyPzMVector sum;
513  for (auto daughter : daughters) {
514  const TrackFitResult* tfr = daughter->getTrackFitResult();
515  if (!tfr) {
516  sum += daughter->get4Vector();
517  continue;
518  }
519  Helix helix = tfr->getHelix();
520  helix.passiveMoveBy(vertex);
521  double scalingFactor = daughter->getEffectiveMomentumScale();
522  double momX = scalingFactor * helix.getMomentumX(bField);
523  double momY = scalingFactor * helix.getMomentumY(bField);
524  double momZ = scalingFactor * helix.getMomentumZ(bField);
525  float mPDG = daughter->getPDGMass();
526  sum += ROOT::Math::PxPyPzMVector(momX, momY, momZ, mPDG);
527  }
528  return sum.M();
529  }
530 
531  double particleInvariantMassLambda(const Particle* part)
532  {
533  const std::vector<Particle*> daughters = part->getDaughters();
534  if (daughters.size() == 2) {
535  ROOT::Math::PxPyPzEVector dt1;
536  ROOT::Math::PxPyPzEVector dt2;
537  ROOT::Math::PxPyPzEVector dtsum;
538  double mpi = Const::pionMass;
539  double mpr = Const::protonMass;
540  dt1 = daughters[0]->get4Vector();
541  dt2 = daughters[1]->get4Vector();
542  double E1 = hypot(mpi, dt1.P());
543  double E2 = hypot(mpr, dt2.P());
544  dtsum = dt1 + dt2;
545  return sqrt((E1 + E2) * (E1 + E2) - dtsum.P() * dtsum.P());
546 
547  } else {
548  return part->getMass();
549  }
550  }
551 
552  double particleInvariantMassError(const Particle* part)
553  {
554  float invMass = part->getMass();
555  TMatrixFSym covarianceMatrix = part->getMomentumErrorMatrix();
556 
557  TVectorF jacobian(Particle::c_DimMomentum);
558  jacobian[0] = -1.0 * part->getPx() / invMass;
559  jacobian[1] = -1.0 * part->getPy() / invMass;
560  jacobian[2] = -1.0 * part->getPz() / invMass;
561  jacobian[3] = 1.0 * part->getEnergy() / invMass;
562 
563  double result = jacobian * (covarianceMatrix * jacobian);
564 
565  if (result < 0.0)
566  return Const::doubleNaN;
567 
568  return TMath::Sqrt(result);
569  }
570 
571  double particleInvariantMassSignificance(const Particle* part)
572  {
573  return particleDMass(part) / particleInvariantMassError(part);
574  }
575 
576  double particleMassSquared(const Particle* part)
577  {
578  ROOT::Math::PxPyPzEVector p4 = part->get4Vector();
579  return p4.M2();
580  }
581 
582  double b2bTheta(const Particle* part)
583  {
584  PCmsLabTransform T;
585  ROOT::Math::PxPyPzEVector pcms = T.rotateLabToCms() * part->get4Vector();
586  ROOT::Math::PxPyPzEVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
587  ROOT::Math::PxPyPzEVector b2blab = T.rotateCmsToLab() * b2bcms;
588  return b2blab.Theta();
589  }
590 
591  double b2bPhi(const Particle* part)
592  {
593  PCmsLabTransform T;
594  ROOT::Math::PxPyPzEVector pcms = T.rotateLabToCms() * part->get4Vector();
595  ROOT::Math::PxPyPzEVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
596  ROOT::Math::PxPyPzEVector b2blab = T.rotateCmsToLab() * b2bcms;
597  return b2blab.Phi();
598  }
599 
600  double b2bClusterTheta(const Particle* part)
601  {
602  // get associated ECLCluster
603  const ECLCluster* cluster = part->getECLCluster();
604  if (!cluster) return Const::doubleNaN;
605  const ECLCluster::EHypothesisBit clusterHypothesis = part->getECLClusterEHypothesisBit();
606 
607  // get 4 momentum from cluster
608  ClusterUtils clutls;
609  ROOT::Math::PxPyPzEVector p4Cluster = clutls.Get4MomentumFromCluster(cluster, clusterHypothesis);
610 
611  // find the vector that balances this in the CMS
612  PCmsLabTransform T;
613  ROOT::Math::PxPyPzEVector pcms = T.rotateLabToCms() * p4Cluster;
614  ROOT::Math::PxPyPzEVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
615  ROOT::Math::PxPyPzEVector b2blab = T.rotateCmsToLab() * b2bcms;
616  return b2blab.Theta();
617  }
618 
619  double b2bClusterPhi(const Particle* part)
620  {
621  // get associated ECLCluster
622  const ECLCluster* cluster = part->getECLCluster();
623  if (!cluster) return Const::doubleNaN;
624  const ECLCluster::EHypothesisBit clusterHypothesis = part->getECLClusterEHypothesisBit();
625 
626  // get 4 momentum from cluster
627  ClusterUtils clutls;
628  ROOT::Math::PxPyPzEVector p4Cluster = clutls.Get4MomentumFromCluster(cluster, clusterHypothesis);
629 
630  // find the vector that balances this in the CMS
631  PCmsLabTransform T;
632  ROOT::Math::PxPyPzEVector pcms = T.rotateLabToCms() * p4Cluster;
633  ROOT::Math::PxPyPzEVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
634  ROOT::Math::PxPyPzEVector b2blab = T.rotateCmsToLab() * b2bcms;
635  return b2blab.Phi();
636  }
637 
638 
639 // released energy --------------------------------------------------
640 
641  double particleQ(const Particle* part)
642  {
643  float m = part->getMass();
644  for (unsigned i = 0; i < part->getNDaughters(); i++) {
645  const Particle* child = part->getDaughter(i);
646  if (child)
647  m -= child->getMass();
648  }
649  return m;
650  }
651 
652  double particleDQ(const Particle* part)
653  {
654  float m = part->getMass() - part->getPDGMass();
655  for (unsigned i = 0; i < part->getNDaughters(); i++) {
656  const Particle* child = part->getDaughter(i);
657  if (child)
658  m -= (child->getMass() - child->getPDGMass());
659  }
660  return m;
661  }
662 
663 // Mbc and deltaE
664 
665  double particleMbc(const Particle* part)
666  {
667  PCmsLabTransform T;
668  ROOT::Math::PxPyPzEVector vec = T.rotateLabToCms() * part->get4Vector();
669  double E = T.getCMSEnergy() / 2;
670  double m2 = E * E - vec.P2();
671  double mbc = m2 >= 0 ? sqrt(m2) : Const::doubleNaN;
672  return mbc;
673  }
674 
675  double particleDeltaE(const Particle* part)
676  {
677  PCmsLabTransform T;
678  ROOT::Math::PxPyPzEVector vec = T.rotateLabToCms() * part->get4Vector();
679  return vec.E() - T.getCMSEnergy() / 2;
680  }
681 
682 // other ------------------------------------------------------------
683 
684 
685  void printParticleInternal(const Particle* p, int depth)
686  {
687  stringstream s("");
688  for (int i = 0; i < depth; i++) {
689  s << " ";
690  }
691  s << p->getPDGCode();
692  const MCParticle* mcp = p->getMCParticle();
693  if (mcp) {
694  unsigned int flags = MCMatching::getMCErrors(p, mcp);
695  s << " -> MC: " << mcp->getPDG() << ", mcErrors: " << flags << " ("
696  << MCMatching::explainFlags(flags) << ")";
697  s << ", mc-index " << mcp->getIndex();
698  } else {
699  s << " (no MC match)";
700  }
701  s << ", mdst-source " << p->getMdstSource();
702  B2INFO(s.str());
703  for (const auto* daughter : p->getDaughters()) {
704  printParticleInternal(daughter, depth + 1);
705  }
706  }
707 
708  bool printParticle(const Particle* p)
709  {
710  printParticleInternal(p, 0);
711  return 0;
712  }
713 
714 
715  double particleMCMomentumTransfer2(const Particle* part)
716  {
717  // for B meson MC particles only
718  const MCParticle* mcB = part->getMCParticle();
719 
720  if (!mcB)
721  return Const::doubleNaN;
722 
723  ROOT::Math::PxPyPzEVector pB = mcB->get4Vector();
724 
725  std::vector<MCParticle*> mcDaug = mcB->getDaughters();
726 
727  if (mcDaug.empty())
728  return Const::doubleNaN;
729 
730  // B -> X l nu
731  // q = pB - pX
732  ROOT::Math::PxPyPzEVector pX;
733 
734  for (auto mcTemp : mcDaug) {
735  if (abs(mcTemp->getPDG()) <= 16)
736  continue;
737 
738  pX += mcTemp->get4Vector();
739  }
740 
741  ROOT::Math::PxPyPzEVector q = pB - pX;
742 
743  return q.M2();
744  }
745 
746 // Recoil Kinematics related ---------------------------------------------
747  double recoilPx(const Particle* particle)
748  {
749  PCmsLabTransform T;
750 
751  // Initial state (e+e- momentum in LAB)
752  ROOT::Math::PxPyPzEVector pIN = T.getBeamFourMomentum();
753 
754  // Use requested frame for final calculation
755  const auto& frame = ReferenceFrame::GetCurrent();
756  return frame.getMomentum(pIN - particle->get4Vector()).Px();
757  }
758 
759  double recoilPy(const Particle* particle)
760  {
761  PCmsLabTransform T;
762 
763  // Initial state (e+e- momentum in LAB)
764  ROOT::Math::PxPyPzEVector pIN = T.getBeamFourMomentum();
765 
766  // Use requested frame for final calculation
767  const auto& frame = ReferenceFrame::GetCurrent();
768  return frame.getMomentum(pIN - particle->get4Vector()).Py();
769  }
770 
771  double recoilPz(const Particle* particle)
772  {
773  PCmsLabTransform T;
774 
775  // Initial state (e+e- momentum in LAB)
776  ROOT::Math::PxPyPzEVector pIN = T.getBeamFourMomentum();
777 
778  // Use requested frame for final calculation
779  const auto& frame = ReferenceFrame::GetCurrent();
780  return frame.getMomentum(pIN - particle->get4Vector()).Pz();
781  }
782 
783  double recoilMomentum(const Particle* particle)
784  {
785  PCmsLabTransform T;
786 
787  // Initial state (e+e- momentum in LAB)
788  ROOT::Math::PxPyPzEVector pIN = T.getBeamFourMomentum();
789 
790  // Use requested frame for final calculation
791  const auto& frame = ReferenceFrame::GetCurrent();
792  return frame.getMomentum(pIN - particle->get4Vector()).P();
793  }
794 
795  double recoilMomentumTheta(const Particle* particle)
796  {
797  PCmsLabTransform T;
798 
799  // Initial state (e+e- momentum in LAB)
800  ROOT::Math::PxPyPzEVector pIN = T.getBeamFourMomentum();
801 
802  // Use requested frame for final calculation
803  const auto& frame = ReferenceFrame::GetCurrent();
804  return frame.getMomentum(pIN - particle->get4Vector()).Theta();
805  }
806 
807  double recoilMomentumPhi(const Particle* particle)
808  {
809  PCmsLabTransform T;
810 
811  // Initial state (e+e- momentum in LAB)
812  ROOT::Math::PxPyPzEVector pIN = T.getBeamFourMomentum();
813 
814  // Use requested frame for final calculation
815  const auto& frame = ReferenceFrame::GetCurrent();
816  return frame.getMomentum(pIN - particle->get4Vector()).Phi();
817  }
818 
819  double recoilEnergy(const Particle* particle)
820  {
821  PCmsLabTransform T;
822 
823  // Initial state (e+e- momentum in LAB)
824  ROOT::Math::PxPyPzEVector pIN = T.getBeamFourMomentum();
825 
826  // Use requested frame for final calculation
827  const auto& frame = ReferenceFrame::GetCurrent();
828  return frame.getMomentum(pIN - particle->get4Vector()).E();
829  }
830 
831  double recoilMass(const Particle* particle)
832  {
833  PCmsLabTransform T;
834 
835  // Initial state (e+e- momentum in LAB)
836  ROOT::Math::PxPyPzEVector pIN = T.getBeamFourMomentum();
837 
838  // Use requested frame for final calculation
839  const auto& frame = ReferenceFrame::GetCurrent();
840  return frame.getMomentum(pIN - particle->get4Vector()).M();
841  }
842 
843  double recoilMassSquared(const Particle* particle)
844  {
845  PCmsLabTransform T;
846 
847  // Initial state (e+e- momentum in LAB)
848  ROOT::Math::PxPyPzEVector pIN = T.getBeamFourMomentum();
849 
850  // Use requested frame for final calculation
851  const auto& frame = ReferenceFrame::GetCurrent();
852  return frame.getMomentum(pIN - particle->get4Vector()).M2();
853  }
854 
855  double m2RecoilSignalSide(const Particle* part)
856  {
857  PCmsLabTransform T;
858  double beamEnergy = T.getCMSEnergy() / 2.;
859  if (part->getNDaughters() != 2) return Const::doubleNaN;
860  ROOT::Math::PxPyPzEVector tagVec = T.rotateLabToCms() * part->getDaughter(0)->get4Vector();
861  ROOT::Math::PxPyPzEVector sigVec = T.rotateLabToCms() * part->getDaughter(1)->get4Vector();
862  tagVec.SetE(-beamEnergy);
863  return (-tagVec - sigVec).M2();
864  }
865 
866  double recoilMCDecayType(const Particle* particle)
867  {
868  auto* mcp = particle->getMCParticle();
869 
870  if (!mcp)
871  return Const::doubleNaN;
872 
873  MCParticle* mcMother = mcp->getMother();
874 
875  if (!mcMother)
876  return Const::doubleNaN;
877 
878  std::vector<MCParticle*> daughters = mcMother->getDaughters();
879 
880  if (daughters.size() != 2)
881  return Const::doubleNaN;
882 
883  MCParticle* recoilMC = nullptr;
884  if (daughters[0]->getArrayIndex() == mcp->getArrayIndex())
885  recoilMC = daughters[1];
886  else
887  recoilMC = daughters[0];
888 
889  if (!recoilMC->hasStatus(MCParticle::c_PrimaryParticle))
890  return Const::doubleNaN;
891 
892  int decayType = 0;
893  checkMCParticleDecay(recoilMC, decayType, false);
894 
895  if (decayType == 0)
896  checkMCParticleDecay(recoilMC, decayType, true);
897 
898  return decayType;
899  }
900 
901  void checkMCParticleDecay(MCParticle* mcp, int& decayType, bool recursive)
902  {
903  int nHadronicParticles = 0;
904  int nPrimaryParticleDaughters = 0;
905  std::vector<MCParticle*> daughters = mcp->getDaughters();
906 
907  // Are any of the daughters primary particles? How many of them are hadrons?
908  for (auto& daughter : daughters) {
909  if (!daughter->hasStatus(MCParticle::c_PrimaryParticle))
910  continue;
911 
912  nPrimaryParticleDaughters++;
913  if (abs(daughter->getPDG()) > Const::photon.getPDGCode())
914  nHadronicParticles++;
915  }
916 
917  if (nPrimaryParticleDaughters > 1) {
918  for (auto& daughter : daughters) {
919  if (!daughter->hasStatus(MCParticle::c_PrimaryParticle))
920  continue;
921 
922  if (abs(daughter->getPDG()) == 12 or abs(daughter->getPDG()) == 14 or abs(daughter->getPDG()) == 16) {
923  if (!recursive) {
924  if (nHadronicParticles == 0) {
925  decayType = 1.0;
926  break;
927  } else {
928  decayType = 2.0;
929  break;
930  }
931  } else {
932  decayType = 3.0;
933  break;
934  }
935  }
936 
937  else if (recursive)
938  checkMCParticleDecay(daughter, decayType, recursive);
939  }
940  }
941  }
942 
943  int nRemainingTracksInEvent(const Particle* particle)
944  {
945 
946  StoreArray<Track> tracks;
947  int event_tracks = tracks.getEntries();
948 
949  int par_tracks = 0;
950  const auto& daughters = particle->getFinalStateDaughters();
951  for (const auto& daughter : daughters) {
952  int pdg = abs(daughter->getPDGCode());
953  if (pdg == Const::electron.getPDGCode() or pdg == Const::muon.getPDGCode() or pdg == Const::pion.getPDGCode()
954  or pdg == Const::kaon.getPDGCode() or pdg == Const::proton.getPDGCode())
955  par_tracks++;
956  }
957  return event_tracks - par_tracks;
958  }
959 
960  double trackMatchType(const Particle* particle)
961  {
962  // Particle does not contain a ECL Cluster
963  double result = Const::doubleNaN;
964 
965  const ECLCluster* cluster = particle->getECLCluster();
966  if (cluster) {
967  // No associated track is default
968  result = 0;
969  if (cluster->isTrack()) {
970  // There is a track match
971  result = 1.0;
972  }
973  }
974  return result;
975  }
976 
977 
978  bool False(const Particle*)
979  {
980  return 0;
981  }
982 
983  bool True(const Particle*)
984  {
985  return 1;
986  }
987 
988  double infinity(const Particle*)
989  {
990  double inf = std::numeric_limits<double>::infinity();
991  return inf;
992  }
993 
994  double random(const Particle*)
995  {
996  return gRandom->Uniform(0, 1);
997  }
998 
999  double eventRandom(const Particle*)
1000  {
1001  std::string key = "__eventRandom";
1002  StoreObjPtr<EventExtraInfo> eventExtraInfo;
1003  if (not eventExtraInfo.isValid())
1004  eventExtraInfo.create();
1005  if (eventExtraInfo->hasExtraInfo(key)) {
1006  return eventExtraInfo->getExtraInfo(key);
1007  } else {
1008  double value = gRandom->Uniform(0, 1);
1009  eventExtraInfo->addExtraInfo(key, value);
1010  return value;
1011  }
1012  }
1013 
1014  Manager::FunctionPtr particleDistToClosestExtTrk(const std::vector<std::string>& arguments)
1015  {
1016  if (arguments.size() != 3 && arguments.size() != 4) {
1017  B2ERROR("Wrong number of arguments (3 or 4 required) for meta variable minET2ETDist");
1018  return nullptr;
1019  }
1020  bool useHighestProbMassForExt(true);
1021  if (arguments.size() == 4) {
1022  try {
1023  useHighestProbMassForExt = static_cast<bool>(Belle2::convertString<int>(arguments[3]));
1024  } catch (std::invalid_argument& e) {
1025  B2ERROR("Fourth (optional) argument of minET2ETDist must be an integer flag.");
1026  return nullptr;
1027  }
1028  }
1029 
1030  std::string detName = arguments[0];
1031  std::string detLayer = arguments[1];
1032  std::string referenceListName = arguments[2];
1033  std::string extraSuffix = (useHighestProbMassForExt) ? "__useHighestProbMassForExt" : "";
1034  // Distance to closets neighbour at this detector layer.
1035  std::string extraInfo = "distToClosestTrkAt" + detName + detLayer + "_VS_" + referenceListName + extraSuffix;
1036 
1037  auto func = [ = ](const Particle * part) -> double {
1038  auto dist = (part->hasExtraInfo(extraInfo)) ? part->getExtraInfo(extraInfo) : Const::doubleNaN;
1039  return dist;
1040  };
1041 
1042  return func;
1043  }
1044 
1045 // Track helix extrapolation-based isolation --------------------------------------------------
1046 
1047  Manager::FunctionPtr particleDistToClosestExtTrkVar(const std::vector<std::string>& arguments)
1048  {
1049  if (arguments.size() != 4) {
1050  B2ERROR("Wrong number of arguments (4 required) for meta variable minET2ETDistVar");
1051  return nullptr;
1052  }
1053 
1054  std::string detName = arguments[0];
1055  std::string detLayer = arguments[1];
1056  std::string referenceListName = arguments[2];
1057  std::string variableName = arguments[3];
1058  // Mdst array index of the particle that is closest to the particle in question.
1059  std::string extraInfo = "idxOfClosestPartAt" + detName + detLayer + "In_" + referenceListName;
1060 
1061  auto func = [ = ](const Particle * part) -> double {
1062 
1063  StoreObjPtr<ParticleList> refPartList(referenceListName);
1064  if (!refPartList.isValid())
1065  {
1066  B2FATAL("Invalid Listname " << referenceListName << " given to minET2ETDistVar!");
1067  }
1068 
1069  if (!part->hasExtraInfo(extraInfo))
1070  {
1071  return Const::doubleNaN;
1072  }
1073 
1074  const Variable::Manager::Var* var = Manager::Instance().getVariable(variableName);
1075  auto refPart = refPartList->getParticleWithMdstIdx(part->getExtraInfo(extraInfo));
1076 
1077  return std::get<double>(var->function(refPart));
1078  };
1079 
1080  return func;
1081  }
1082 
1083 
1084  Manager::FunctionPtr particleExtTrkIsoScoreVar(const std::vector<std::string>& arguments)
1085  {
1086 
1087  if (arguments.size() < 3) {
1088  B2ERROR("Wrong number of arguments (at least 3 required) for meta variable minET2ETIsoScore");
1089  return nullptr;
1090  }
1091 
1092  std::string referenceListName = arguments[0];
1093  bool useHighestProbMassForExt;
1094  try {
1095  useHighestProbMassForExt = static_cast<bool>(Belle2::convertString<int>(arguments[1]));
1096  } catch (std::invalid_argument& e) {
1097  B2ERROR("Second argument must be an integer flag.");
1098  return nullptr;
1099  }
1100  std::string extraSuffix = (useHighestProbMassForExt) ? "__useHighestProbMassForExt" : "";
1101 
1102  std::vector<std::string> detectorNames(arguments.begin() + 2, arguments.end());
1103 
1104  std::string detNamesConcat("");
1105  for (auto& detName : detectorNames) {
1106  boost::to_upper(detName);
1107  detNamesConcat += "_" + detName;
1108  }
1109 
1110  std::string extraInfo = "trkIsoScore" + detNamesConcat + "_VS_" + referenceListName + extraSuffix;
1111 
1112  auto func = [ = ](const Particle * part) -> double {
1113 
1114  StoreObjPtr<ParticleList> refPartList(referenceListName);
1115  if (!refPartList.isValid())
1116  {
1117  B2FATAL("Invalid Listname " << referenceListName << " given to minET2ETIsoScore!");
1118  }
1119 
1120  if (!part->hasExtraInfo(extraInfo))
1121  {
1122  return Const::doubleNaN;
1123  }
1124  auto scoreDet = part->getExtraInfo(extraInfo);
1125  if (std::isnan(scoreDet))
1126  {
1127  return Const::doubleNaN;
1128  }
1129 
1130  return scoreDet;
1131 
1132  };
1133 
1134  return func;
1135  }
1136 
1137 
1138  Manager::FunctionPtr particleExtTrkIsoScoreVarAsWeightedAvg(const std::vector<std::string>& arguments)
1139  {
1140 
1141  if (arguments.size() < 3) {
1142  B2ERROR("Wrong number of arguments (at least 3 required) for meta variable minET2ETIsoScoreAsWeightedAvg");
1143  return nullptr;
1144  }
1145 
1146  std::string referenceListName = arguments[0];
1147  bool useHighestProbMassForExt;
1148  try {
1149  useHighestProbMassForExt = static_cast<bool>(Belle2::convertString<int>(arguments[1]));
1150  } catch (std::invalid_argument& e) {
1151  B2ERROR("Second argument must be an integer flag.");
1152  return nullptr;
1153  }
1154  std::string extraSuffix = (useHighestProbMassForExt) ? "__useHighestProbMassForExt" : "";
1155 
1156  std::vector<std::string> detectorNames(arguments.begin() + 2, arguments.end());
1157 
1158  std::string detNamesConcat("");
1159  for (auto& detName : detectorNames) {
1160  boost::to_upper(detName);
1161  detNamesConcat += "_" + detName;
1162  }
1163 
1164  std::string extraInfo = "trkIsoScoreAsWeightedAvg" + detNamesConcat + "_VS_" + referenceListName + extraSuffix;
1165 
1166  auto func = [ = ](const Particle * part) -> double {
1167 
1168  StoreObjPtr<ParticleList> refPartList(referenceListName);
1169  if (!refPartList.isValid())
1170  {
1171  B2FATAL("Invalid Listname " << referenceListName << " given to minET2ETIsoScoreAsWeightedAvg!");
1172  }
1173 
1174  if (!part->hasExtraInfo(extraInfo))
1175  {
1176  return Const::doubleNaN;
1177  }
1178  auto scoreDet = part->getExtraInfo(extraInfo);
1179  if (std::isnan(scoreDet))
1180  {
1181  return Const::doubleNaN;
1182  }
1183 
1184  return scoreDet;
1185 
1186  };
1187 
1188  return func;
1189  }
1190 
1191  VARIABLE_GROUP("Kinematics");
1192  REGISTER_VARIABLE("p", particleP, "momentum magnitude\n\n", "GeV/c");
1193  REGISTER_VARIABLE("E", particleE, "energy\n\n", "GeV");
1194 
1195  REGISTER_VARIABLE("E_uncertainty", particleEUncertainty, R"DOC(
1196  energy uncertainty (:math:`\sqrt{\sigma^2}`)
1197 
1198  )DOC", "GeV");
1199  REGISTER_VARIABLE("ECLClusterE_uncertainty", particleClusterEUncertainty,
1200  "energy uncertainty as given by the underlying ECL cluster\n\n", "GeV");
1201  REGISTER_VARIABLE("px", particlePx, "momentum component x\n\n", "GeV/c");
1202  REGISTER_VARIABLE("py", particlePy, "momentum component y\n\n", "GeV/c");
1203  REGISTER_VARIABLE("pz", particlePz, "momentum component z\n\n", "GeV/c");
1204  REGISTER_VARIABLE("pt", particlePt, "transverse momentum\n\n", "GeV/c");
1205  REGISTER_VARIABLE("xp", particleXp,
1206  "scaled momentum: the momentum of the particle in the CMS as a fraction of its maximum available momentum in the collision");
1207  REGISTER_VARIABLE("pErr", particlePErr, "error of momentum magnitude\n\n", "GeV/c");
1208  REGISTER_VARIABLE("pxErr", particlePxErr, "error of momentum component x\n\n", "GeV/c");
1209  REGISTER_VARIABLE("pyErr", particlePyErr, "error of momentum component y\n\n", "GeV/c");
1210  REGISTER_VARIABLE("pzErr", particlePzErr, "error of momentum component z\n\n", "GeV/c");
1211  REGISTER_VARIABLE("ptErr", particlePtErr, "error of transverse momentum\n\n", "GeV/c");
1212  REGISTER_VARIABLE("momVertCovM(i,j)", covMatrixElement,
1213  "returns the (i,j)-th element of the MomentumVertex Covariance Matrix (7x7).\n"
1214  "Order of elements in the covariance matrix is: px, py, pz, E, x, y, z.\n\n", "GeV/c, GeV/c, GeV/c, GeV, cm, cm, cm");
1215  REGISTER_VARIABLE("momDevChi2", momentumDeviationChi2, R"DOC(
1216 momentum deviation :math:`\chi^2` value calculated as :math:`\chi^2 = \sum_i (p_i - mc(p_i))^2/\sigma(p_i)^2`,
1217 where :math:`\sum` runs over i = px, py, pz and :math:`mc(p_i)` is the mc truth value and :math:`\sigma(p_i)` is the estimated error of i-th component of momentum vector
1218 )DOC");
1219  REGISTER_VARIABLE("theta", particleTheta, "polar angle\n\n", "rad");
1220  REGISTER_VARIABLE("thetaErr", particleThetaErr, "error of polar angle\n\n", "rad");
1221  REGISTER_VARIABLE("cosTheta", particleCosTheta, "momentum cosine of polar angle");
1222  REGISTER_VARIABLE("cosThetaErr", particleCosThetaErr, "error of momentum cosine of polar angle");
1223  REGISTER_VARIABLE("phi", particlePhi, "momentum azimuthal angle\n\n", "rad");
1224  REGISTER_VARIABLE("phiErr", particlePhiErr, "error of momentum azimuthal angle\n\n", "rad");
1225  REGISTER_VARIABLE("PDG", particlePDGCode, "PDG code");
1226  REGISTER_VARIABLE("cosAngleBetweenMomentumAndVertexVectorInXYPlane",
1227  cosAngleBetweenMomentumAndVertexVectorInXYPlane,
1228  "cosine of the angle between momentum and vertex vector (vector connecting ip and fitted vertex) of this particle in xy-plane");
1229  REGISTER_VARIABLE("cosAngleBetweenMomentumAndVertexVector",
1230  cosAngleBetweenMomentumAndVertexVector,
1231  "cosine of the angle between momentum and vertex vector (vector connecting ip and fitted vertex) of this particle");
1232  REGISTER_VARIABLE("cosThetaBetweenParticleAndNominalB",
1233  cosThetaBetweenParticleAndNominalB,
1234  "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.");
1235  REGISTER_VARIABLE("cosToThrustOfEvent", cosToThrustOfEvent,
1236  "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");
1237 
1238  REGISTER_VARIABLE("ImpactXY" , ImpactXY , "The impact parameter of the given particle in the xy plane\n\n", "cm");
1239 
1240  REGISTER_VARIABLE("M", particleMass, R"DOC(
1241 The particle's mass.
1242 
1243 Note that this is context-dependent variable and can take different values depending on the situation. This should be the "best"
1244 value possible with the information provided.
1245 
1246 - If this particle is track- or cluster- based, then this is the value of the mass hypothesis.
1247 - If this particle is an MC particle then this is the mass of that particle.
1248 - If this particle is composite, then *initially* this takes the value of the invariant mass of the daughters.
1249 - If this particle is composite and a *mass or vertex fit* has been performed then this may be updated by the fit.
1250 
1251 * You will see a difference between this mass and the :b2:var:`InvM`.
1252 
1253 )DOC", "GeV/:math:`\\text{c}^2`");
1254  REGISTER_VARIABLE("dM", particleDMass, "mass minus nominal mass\n\n", "GeV/:math:`\\text{c}^2`");
1255  REGISTER_VARIABLE("Q", particleQ, "energy released in decay\n\n", "GeV");
1256  REGISTER_VARIABLE("dQ", particleDQ, ":b2:var:`Q` minus nominal energy released in decay\n\n", "GeV");
1257  REGISTER_VARIABLE("Mbc", particleMbc, "beam constrained mass\n\n", "GeV/:math:`\\text{c}^2`");
1258  REGISTER_VARIABLE("deltaE", particleDeltaE, "difference between :b2:var:`E` and half the center of mass energy\n\n", "GeV");
1259  REGISTER_VARIABLE("M2", particleMassSquared, "The particle's mass squared.\n\n", ":math:`[\\text{GeV}/\\text{c}^2]^2`");
1260 
1261  REGISTER_VARIABLE("InvM", particleInvariantMassFromDaughtersDisplaced,
1262  "invariant mass (determined from particle's daughter 4-momentum vectors). If this particle is V0 or decays at rho > 5 mm, its daughter 4-momentum vectors at fitted vertex are taken.\n"
1263  "If this particle has no daughters, defaults to :b2:var:`M`.\n\n", "GeV/:math:`\\text{c}^2`");
1264  REGISTER_VARIABLE("InvMLambda", particleInvariantMassLambda,
1265  "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"
1266  "If the particle has not 2 daughters, it returns just the mass value.\n\n", "GeV/:math:`\\text{c}^2`");
1267 
1268  REGISTER_VARIABLE("ErrM", particleInvariantMassError,
1269  "uncertainty of invariant mass\n\n", "GeV/:math:`\\text{c}^2`");
1270  REGISTER_VARIABLE("SigM", particleInvariantMassSignificance,
1271  "signed deviation of particle's invariant mass from its nominal mass in units of the uncertainty on the invariant mass (:b2:var:`dM`/:b2:var:`ErrM`)");
1272 
1273  REGISTER_VARIABLE("pxRecoil", recoilPx,
1274  "component x of 3-momentum recoiling against given Particle\n\n", "GeV/c");
1275  REGISTER_VARIABLE("pyRecoil", recoilPy,
1276  "component y of 3-momentum recoiling against given Particle\n\n", "GeV/c");
1277  REGISTER_VARIABLE("pzRecoil", recoilPz,
1278  "component z of 3-momentum recoiling against given Particle\n\n", "GeV/c");
1279 
1280  REGISTER_VARIABLE("pRecoil", recoilMomentum,
1281  "magnitude of 3 - momentum recoiling against given Particle\n\n", "GeV/c");
1282  REGISTER_VARIABLE("pRecoilTheta", recoilMomentumTheta,
1283  "Polar angle of a particle's missing momentum\n\n", "rad");
1284  REGISTER_VARIABLE("pRecoilPhi", recoilMomentumPhi,
1285  "Azimuthal angle of a particle's missing momentum\n\n", "rad");
1286  REGISTER_VARIABLE("eRecoil", recoilEnergy,
1287  "energy recoiling against given Particle\n\n", "GeV");
1288  REGISTER_VARIABLE("mRecoil", recoilMass,
1289  "Invariant mass of the system recoiling against given Particle\n\n", "GeV/:math:`\\text{c}^2`");
1290  REGISTER_VARIABLE("m2Recoil", recoilMassSquared,
1291  "invariant mass squared of the system recoiling against given Particle\n\n", ":math:`[\\text{GeV}/\\text{c}^2]^2`");
1292  REGISTER_VARIABLE("m2RecoilSignalSide", m2RecoilSignalSide, R"DOC(
1293  Squared recoil mass of the signal side which is calculated in the CMS frame under the assumption that the
1294  signal and tag side are produced back to back and the tag side energy equals the beam energy. The variable
1295  must be applied to the Upsilon and the tag side must be the first, the signal side the second daughter
1296 
1297  )DOC", ":math:`[\\text{GeV}/\\text{c}^2]^2`");
1298 
1299  REGISTER_VARIABLE("b2bTheta", b2bTheta,
1300  "Polar angle in the lab system that is back-to-back to the particle in the CMS. Useful for low multiplicity studies.\n\n", "rad");
1301  REGISTER_VARIABLE("b2bPhi", b2bPhi,
1302  "Azimuthal angle in the lab system that is back-to-back to the particle in the CMS. Useful for low multiplicity studies.\n\n",
1303  "rad");
1304  REGISTER_VARIABLE("b2bClusterTheta", b2bClusterTheta,
1305  "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.\n\n",
1306  "rad");
1307  REGISTER_VARIABLE("b2bClusterPhi", b2bClusterPhi,
1308  "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.\n\n",
1309  "rad");
1310  REGISTER_VARIABLE("ArmenterosLongitudinalMomentumAsymmetry", ArmenterosLongitudinalMomentumAsymmetry,
1311  "Longitudinal momentum asymmetry of V0's daughters.\n"
1312  "The mother (V0) is required to have exactly two daughters");
1313  REGISTER_VARIABLE("ArmenterosDaughter1Qt", ArmenterosDaughter1Qt, R"DOC(
1314  Transverse momentum of the first daughter with respect to the V0 mother.
1315  The mother is required to have exactly two daughters
1316 
1317  )DOC", "GeV/c");
1318  REGISTER_VARIABLE("ArmenterosDaughter2Qt", ArmenterosDaughter2Qt, R"DOC(
1319  Transverse momentum of the second daughter with respect to the V0 mother.
1320  The mother is required to have exactly two daughters
1321 
1322  )DOC", "GeV/c");
1323 
1324  VARIABLE_GROUP("Miscellaneous");
1325  REGISTER_VARIABLE("nRemainingTracksInEvent", nRemainingTracksInEvent,
1326  "Number of tracks in the event - Number of tracks( = charged FSPs) of particle.");
1327  REGISTER_VARIABLE("trackMatchType", trackMatchType, R"DOC(
1328 
1329 * -1 particle has no ECL cluster
1330 * 0 particle has no associated track
1331 * 1 there is a matched track called connected - region(CR) track match
1332 
1333  )DOC");
1334  MAKE_DEPRECATED("trackMatchType", true, "light-2012-minos", R"DOC(
1335  Use better variables like `trackNECLClusters`, `clusterTrackMatch`, and `nECLClusterTrackMatches`.)DOC");
1336 
1337  REGISTER_VARIABLE("decayTypeRecoil", recoilMCDecayType,
1338  "type of the particle decay(no related mcparticle = -1, hadronic = 0, direct leptonic = 1, direct semileptonic = 2,"
1339  "lower level leptonic = 3.");
1340 
1341  REGISTER_VARIABLE("printParticle", printParticle,
1342  "For debugging, print Particle and daughter PDG codes, plus MC match. Returns 0.");
1343  REGISTER_VARIABLE("mcMomTransfer2", particleMCMomentumTransfer2,
1344  "Return the true momentum transfer to lepton pair in a B(semi -) leptonic B meson decay.\n\n", "GeV/c");
1345  REGISTER_VARIABLE("False", False,
1346  "returns always 0, used for testing and debugging.");
1347  REGISTER_VARIABLE("True", True,
1348  "returns always 1, used for testing and debugging.");
1349  REGISTER_VARIABLE("infinity", infinity,
1350  "returns std::numeric_limits<double>::infinity()");
1351  REGISTER_VARIABLE("random", random,
1352  "return a random number between 0 and 1 for each candidate. Can be used, e.g. for picking a random"
1353  "candidate in the best candidate selection.");
1354  REGISTER_VARIABLE("eventRandom", eventRandom,
1355  "[Eventbased] Returns a random number between 0 and 1 for this event. Can be used, e.g. for applying an event prescale.");
1356  REGISTER_METAVARIABLE("minET2ETDist(detName, detLayer, referenceListName, useHighestProbMassForExt=1)", particleDistToClosestExtTrk,
1357  R"DOC(Returns the distance :math:`d_{\mathrm{i}}` in [cm] between the particle and the nearest particle in the reference list at the given detector :math:`i`-th layer surface.
1358 The definition is based on the track helices extrapolation.
1359 
1360 * The first argument is the name of the detector to consider.
1361 * The second argument is the detector layer on whose surface we search for the nearest neighbour.
1362 * The third argument is the reference particle list name used to search for the nearest neighbour.
1363 * The fourth (optional) argument is an integer ("boolean") flag: if 1 (the default, if nothing is set), it is assumed the extrapolation was done with the most probable mass hypothesis for the track fit;
1364  if 0, it is assumed the mass hypothesis matching the particle lists' PDG was used.
1365 
1366 .. note::
1367  This variable requires to run the ``TrackIsoCalculator`` module first.
1368  Note that the choice of input parameters of this metafunction must correspond to the settings used to configure the module!
1369 )DOC",
1370  Manager::VariableDataType::c_double);
1371 
1372  REGISTER_METAVARIABLE("minET2ETDistVar(detName, detLayer, referenceListName, variableName)", particleDistToClosestExtTrkVar,
1373  R"DOC(Returns the value of the variable for the nearest neighbour to this particle as taken from the reference list at the given detector :math:`i`-th layer surface
1374 , according to the distance definition of `minET2ETDist`.
1375 
1376 * The first argument is the name of the detector to consider.
1377 * The second argument is the detector layer on whose surface we search for the nearest neighbour.
1378 * The third argument is the reference particle list name used to search for the nearest neighbour.
1379 * The fourth argument is a variable name, e.g. `nCDCHits`.
1380 
1381 .. note::
1382  This variable requires to run the ``TrackIsoCalculator`` module first.
1383  Note that the choice of input parameters of this metafunction must correspond to the settings used to configure the module!
1384 )DOC",
1385  Manager::VariableDataType::c_double);
1386 
1387  REGISTER_METAVARIABLE("minET2ETIsoScore(referenceListName, useHighestProbMassForExt, detectorList)", particleExtTrkIsoScoreVar,
1388  R"DOC(Returns a particle's isolation score :math:`s` defined as:
1389 
1390 .. math::
1391  :nowrap:
1392 
1393  \begin{split}
1394 
1395  s &= \sum_{\mathrm{det}} 1 - \left(-w_{\mathrm{det}} \cdot \frac{\sum_{i}^{N_{\mathrm{det}}^{\mathrm{layers}}} H(i)}{N_{\mathrm{det}}^{\mathrm{layers}}}\right), \\
1396 
1397  H(i) &=
1398  \begin{cases}
1399  0 & d_{\mathrm{i}} > D_{\mathrm{det}}^{\mathrm{thresh}} \\
1400  1 & d_{\mathrm{i}} <= D_{\mathrm{det}}^{\mathrm{thresh}}, \\
1401  \end{cases}
1402 
1403  \end{split}
1404 
1405 where :math:`d_{\mathrm{i}}` is the distance to the closest neighbour at the :math:`i`-th layer of the given detector (c.f., `minET2ETDist`), :math:`N_{\mathrm{det}}^{\mathrm{layers}}` is the
1406 number of layers of the detector, :math:`D_{\mathrm{det}}^{\mathrm{thresh}}` is a threshold length related to the detector's granularity defined in the ``TrackIsoCalculator`` module,
1407 and :math:`w_{\mathrm{det}}` are (negative) weights associated to the detector's impact on PID for this particle type, read from a CDB payload.
1408 
1409 The score is normalised in [0, 1], where values closer to 1 indicates a well-isolated particle.
1410 
1411 * The first argument is the reference particle list name used to search for the nearest neighbour.
1412 * The second argument is an integer ("boolean") flag: if 1, it is assumed the extrapolation was done with the most probable mass hypothesis for the track fit;
1413  if 0, it is assumed the mass hypothesis matching the particle lists' PDG was used.
1414 * The remaining arguments are a comma-separated list of detector names, which must correspond to the one given to the `TrackIsoCalculator` module.
1415 
1416 .. note::
1417  The PID detector weights :math:`w_{\mathrm{det}}` are non-trivial only if ``excludePIDDetWeights=false`` in the ``TrackIsoCalculator`` module configuration.
1418  Otherwise :math:`w_{\mathrm{det}} = -1`.
1419 
1420 .. note::
1421  This variable requires to run the `TrackIsoCalculator` module first.
1422  Note that the choice of input parameters of this metafunction must correspond to the settings used to configure the module!
1423 
1424 )DOC",
1425  Manager::VariableDataType::c_double);
1426 
1427  REGISTER_METAVARIABLE("minET2ETIsoScoreAsWeightedAvg(referenceListName, useHighestProbMassForExt, detectorList)", particleExtTrkIsoScoreVarAsWeightedAvg,
1428  R"DOC(Returns a particle's isolation score :math:`s` based on the weighted average:
1429 
1430 .. math::
1431 
1432  s = 1 - \frac{\sum_{\mathrm{det}} \sum_{i}^{N_{\mathrm{det}}^{\mathrm{layers}}} w_{\mathrm{det}} \cdot \frac{D_{\mathrm{det}}^{\mathrm{thresh}}}{d_{\mathrm{i}}} }{ \sum_{\mathrm{det}} w_{\mathrm{det}}},
1433 
1434 where :math:`d_{\mathrm{i}}` is the distance to the closest neighbour at the :math:`i`-th layer of the given detector (c.f., `minET2ETDist`), :math:`N_{\mathrm{det}}^{\mathrm{layers}}` is the
1435 number of layers of the detector, :math:`D_{\mathrm{det}}^{\mathrm{thresh}}` is a threshold length related to the detector's granularity defined in the ``TrackIsoCalculator`` module,
1436 and :math:`w_{\mathrm{det}}` are (negative) weights associated to the detector's impact on PID for this particle type, read from a CDB payload.
1437 
1438 The score is normalised in [0, 1], where values closer to 1 indicates a well-isolated particle.
1439 
1440 * The first argument is the reference particle list name used to search for the nearest neighbour.
1441 * The second argument is an integer ("boolean") flag: if 1, it is assumed the extrapolation was done with the most probable mass hypothesis for the track fit;
1442  if 0, it is assumed the mass hypothesis matching the particle lists' PDG was used.
1443 * The remaining arguments are a comma-separated list of detector names, which must correspond to the one given to the `TrackIsoCalculator` module.
1444 
1445 .. note::
1446  The PID detector weights :math:`w_{\mathrm{det}}` are non-trivial only if ``excludePIDDetWeights=false`` in the ``TrackIsoCalculator`` module configuration.
1447  Otherwise :math:`\lvert w_{\mathrm{det}} \rvert = 1`.
1448 
1449 .. note::
1450  This variable requires to run the `TrackIsoCalculator` module first.
1451  Note that the choice of input parameters of this metafunction must correspond to the settings used to configure the module!
1452 
1453 )DOC",
1454  Manager::VariableDataType::c_double);
1455 
1456  }
1458 }
E
R E
internal precision of FFTW codelets
Definition: DiscreteCosineTransform_31points.cc:36
Belle2::B2Vector3::Perp
DataType Perp() const
The transverse component (R in cylindrical coordinate system).
Definition: B2Vector3.h:200
LogVar
Class to store variables with their name which were sent to the logging service.
Definition: LogVariableStream.h:24
MAKE_DEPRECATED
#define MAKE_DEPRECATED(name, make_fatal, version, description)
Registers a variable as deprecated.
Definition: Manager.h:443
Belle2::B2Vector3D
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17