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