Belle II Software development
SpecificKinematicVariables.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/SpecificKinematicVariables.h>
11
12// include VariableManager
13#include <analysis/VariableManager/Manager.h>
14
15#include <analysis/utility/PCmsLabTransform.h>
16#include <analysis/variables/Variables.h>
17
18// dataobjects
19#include <analysis/dataobjects/Particle.h>
20
21// framework aux
22#include <framework/logging/Logger.h>
23
24#include <TRandom.h>
25#include <TMath.h>
26#include <Math/Vector3D.h>
27#include <Math/Vector4D.h>
28using namespace ROOT::Math;
29
30#include <iostream>
31using namespace std;
32
33namespace Belle2 {
38 namespace Variable {
39
40
41 double REC_q2BhSimple(const Particle* particle)
42 {
43 // calculates q^2 = (p_B - p_h) in decays of B -> h_1 .. h_n ell nu_ell,
44 // where p_h = Sum_i^n p_h_i is the 4-momentum of hadrons in the final
45 // state. The calculation is performed in the CMS system, where B-meson
46 // is assumed to be at rest p_B = (m_B, 0).
47
48 PxPyPzEVector hadron4vec;
49
50 unsigned n = particle->getNDaughters();
51
52 if (n < 1)
53 return Const::doubleNaN;
54
55 for (unsigned i = 0; i < n; i++) {
56 int absPDG = abs(particle->getDaughter(i)->getPDGCode());
57 if (absPDG == Const::electron.getPDGCode() || absPDG == Const::muon.getPDGCode() || absPDG == 15)
58 continue;
59
60 hadron4vec += particle->getDaughter(i)->get4Vector();
61 }
62
63 // boost to CMS
64 PCmsLabTransform T;
65 PxPyPzEVector phCMS = T.rotateLabToCms() * hadron4vec;
66 PxPyPzMVector pBCMS(0.0, 0.0, 0.0, particle->getPDGMass());
67 return (pBCMS - phCMS).mag2();
68 }
69
70 double REC_q2Bh(const Particle* particle)
71 {
72 // calculates q^2 = (p_B - p_h) in decays of B -> h_1 .. h_n ell nu_ell,
73 // where p_h = Sum_i^n p_h_i is the 4-momentum of hadrons in the final
74 // state. The calculation is performed in the CMS system,
75 // with a weighter average in a cone around the true B direction
76
77 PxPyPzEVector hadron4vec;
78
79 unsigned n = particle->getNDaughters();
80
81 if (n < 1)
82 return Const::doubleNaN;
83
84 for (unsigned i = 0; i < n; i++) {
85 int absPDG = abs(particle->getDaughter(i)->getPDGCode());
86 if (absPDG == Const::electron.getPDGCode() || absPDG == Const::muon.getPDGCode() || absPDG == 15)
87 continue;
88
89 hadron4vec += particle->getDaughter(i)->get4Vector();
90 }
91
92 // boost to CMS
93 PCmsLabTransform T;
94 PxPyPzEVector had_cm = T.rotateLabToCms() * hadron4vec;
95 PxPyPzEVector Y_cm = T.rotateLabToCms() * particle->get4Vector();
96
97
98 double bmass = particle->getPDGMass();
99 // B theta angle
100 double cos_cone_angle = Variable::cosThetaBetweenParticleAndNominalB(particle);
101 if (abs(cos_cone_angle) > 1) {
102 //makes no sense in this case, return simple value
103 double q2 = Variable::REC_q2BhSimple(particle);
104 if (q2 < 0) {
105 return 0;
106 }
107 return q2;
108 }
109 double thetaBY = TMath::ACos(cos_cone_angle);
110 const double E_B = T.getCMSEnergy() / 2.0;
111 const double p_B = sqrt(E_B * E_B - bmass * bmass);
112
113 double phi_start = gRandom->Uniform(0, TMath::Pi() / 2);
114
115 double q2 = 0;
116 double denom = 0;
117
118 XYZVector zHatY(Y_cm.Vect().Unit());
119 XYZVector yHatY((had_cm.Vect().Cross(zHatY)).Unit());
120 XYZVector xHatY(yHatY.Cross(zHatY));
121
122 double phi = phi_start;
123 double dphi = TMath::Pi() / 2;
124 PxPyPzEVector static_B(0, 0, 0, bmass);
125
126 for (int around_the_cone = 0; around_the_cone < 4; around_the_cone++) {
127 //Define the momentum of B in the Y rest frame using the angles thetaBY and the
128 //current phi in the loop.
129
130 //The purpose of these lines is to calculate the B momentum in the BY cone in the coordinate system previously developed.
131 double B0_px_Y_frame = p_B * TMath::Sin(thetaBY) * TMath::Cos(phi);
132 double B0_py_Y_frame = p_B * TMath::Sin(thetaBY) * TMath::Sin(phi);
133 double B0_pz_Y_frame = p_B * TMath::Cos(thetaBY);
134
135 //The 3 components of the guess for the B0 p direction at the current phi are
136 //calculated by scaling the 3 basis vectors computed before by the corresponding B
137 //momentum in that direction.
138
139 XYZVector B0_px_Dframe(xHatY);
140 B0_px_Dframe *= B0_px_Y_frame;
141 XYZVector B0_py_Dframe(yHatY);
142 B0_py_Dframe *= B0_py_Y_frame;
143 XYZVector B0_pz_Dframe(zHatY);
144 B0_pz_Dframe *= B0_pz_Y_frame;
145 //Construct the B0 p 3-vector with the current phi by summing the 3 components.
146
147 XYZVector B0_p3_Dframe(B0_px_Dframe + B0_py_Dframe + B0_pz_Dframe);
148 PxPyPzEVector B0_p4_Dframe(B0_p3_Dframe.X(), B0_p3_Dframe.Y(), B0_p3_Dframe.Z(), E_B);
149
150 //This is the polar angle of B0.
151
152 double cosThetaB = cos(B0_p3_Dframe.Theta());
153 double sinThetaB2 = (1 - cosThetaB * cosThetaB);
154
155 //The weight is given by the sin squared of such angle.
156 double wt = sinThetaB2;
157
158 //Boost the hadronic daughter to the computed B0 rest frame.
159 // In that frame, q2 can simply be calculated by subtracting the hadron 4-momentum from the momentum of a static B.
160 PxPyPzEVector had_B0(had_cm);
161 had_B0 = Boost(B0_p4_Dframe.BoostToCM()) * had_B0;
162 q2 += wt * ((static_B - had_B0).M2());
163 denom += wt;
164 phi += dphi;
165 }
166
167 q2 /= denom;
168 if (q2 < 0) {
169 q2 = 0.0;
170 }
171 return q2;
172 }
173
174 double REC_MissM2(const Particle* particle)
175 {
176 PCmsLabTransform T;
177 PxPyPzEVector rec4vecLAB = particle->get4Vector();
178 PxPyPzEVector rec4vec = T.rotateLabToCms() * rec4vecLAB;
179
180 double E_beam_cms = T.getCMSEnergy() / 2.0;
181
182 PxPyPzEVector miss4vec(-rec4vec.px(), -rec4vec.py(), -rec4vec.pz(), E_beam_cms - rec4vec.E());
183
184 return miss4vec.mag2();
185 }
186
187
188 VARIABLE_GROUP("Specific kinematic variables");
189
190 REGISTER_VARIABLE("recQ2BhSimple", REC_q2BhSimple,
191 "Returns the momentum transfer squared, :math:`q^2`, calculated in CMS as :math:`q^2 = (p_B - p_h)^2`, \n"
192 "where p_h is the CMS momentum of all hadrons in the decay :math:`B \\to H_1 ... H_n \\ell \\nu_\\ell`.\n"
193 "The B meson momentum in CMS is assumed to be 0.\n\n", ":math:`[\\text{GeV}/\\text{c}]^2`");
194
195 REGISTER_VARIABLE("recQ2Bh", REC_q2Bh,
196 "Returns the momentum transfer squared, :math:`q^2`, calculated in CMS as :math:`q^2 = (p_B - p_h)^2`, \n"
197 "where p_h is the CMS momentum of all hadrons in the decay :math:`B \\to H_1\\dots H_n \\ell \\nu_\\ell`.\n"
198 "This calculation uses a weighted average of the B meson around the reco B cone. \n"
199 "Based on diamond frame calculation of :math:`q^2` following the idea presented in https://www.osti.gov/biblio/1442697 \n"
200 "It will switch to use of :b2:var:`recQ2BhSimple` if absolute of :b2:var:`cosThetaBetweenParticleAndNominalB` > 1.\n\n",
201 ":math:`[\\text{GeV}/\\text{c}]^2`");
202
203 REGISTER_VARIABLE("recMissM2", REC_MissM2, R"DOC(
204 Returns the invariant mass squared of the missing momentum calculated assumings the
205 reco B is at rest and calculating the neutrino (missing) momentum from :math:`p_\nu = p_B - p_{\rm had} - p_{\rm lep}`
206
207 )DOC", ":math:`[\\text{GeV}/\\text{c}^2]^2`");
208
209
210
211 }
213}
static const ChargedStable muon
muon particle
Definition: Const.h:660
static const double doubleNaN
quiet_NaN
Definition: Const.h:703
static const ChargedStable electron
electron particle
Definition: Const.h:659
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Abstract base class for different kinds of events.
STL namespace.