Belle II Software development
BelleVariables.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/BelleVariables.h>
11
12// include VariableManager
13#include <analysis/VariableManager/Manager.h>
14
15#include <analysis/dataobjects/Particle.h>
16#include <analysis/variables/Variables.h>
17#include <analysis/variables/VertexVariables.h>
18#include <analysis/variables/ECLVariables.h>
19#include <analysis/variables/V0DaughterTrackVariables.h>
20#include <mdst/dataobjects/Track.h>
21
22#include <framework/logging/Logger.h>
23#include <framework/gearbox/Const.h>
24
25#include <b2bii/dataobjects/BelleTrkExtra.h>
26
27#include <cmath>
28#include <TVectorF.h>
29
30namespace Belle2 {
35 namespace Variable {
36 bool goodBelleKshort(const Particle* KS)
37 {
38 // check input
39 if (KS->getNDaughters() != 2) {
40 B2WARNING("goodBelleKshort is only defined for a particle with two daughters");
41 return false;
42 }
43 const Particle* d0 = KS->getDaughter(0);
44 const Particle* d1 = KS->getDaughter(1);
45 if ((d0->getCharge() == 0) || (d1->getCharge() == 0)) {
46 B2WARNING("goodBelleKshort is only defined for a particle with charged daughters");
47 return false;
48 }
49 if (std::abs(KS->getPDGCode()) != Const::Kshort.getPDGCode())
50 B2WARNING("goodBelleKshort is being applied to a candidate with PDG " << KS->getPDGCode());
51
52 // If goodKs exists, return the value
53 if (KS->hasExtraInfo("goodKs")) {
54 return bool(KS->getExtraInfo("goodKs"));
55 }
56
57 // Belle selection
58 double p = particleP(KS);
59 double fl = particleDRho(KS);
60 double dphi = acos(((particleDX(KS) * particlePx(KS)) + (particleDY(KS) * particlePy(KS))) / (fl * sqrt(particlePx(KS) * particlePx(
61 KS) + particlePy(KS) * particlePy(KS))));
62 // particleDRho returns track d0 relative to IP for tracks
63 double dr = std::min(std::abs(particleDRho(d0)), std::abs(particleDRho(d1)));
64 double zdist = v0DaughterZ0Diff(KS);
65
66 bool low = p < 0.5 && std::abs(zdist) < 0.8 && dr > 0.05 && dphi < 0.3;
67 bool mid = p < 1.5 && p > 0.5 && std::abs(zdist) < 1.8 && dr > 0.03 && dphi < 0.1 && fl > .08;
68 bool high = p > 1.5 && std::abs(zdist) < 2.4 && dr > 0.02 && dphi < 0.03 && fl > .22;
69
70 if (low || mid || high) {
71 return true;
72 } else
73 return false;
74 }
75
76
77 double goodBelleLambda(const Particle* Lambda)
78 {
79 if (Lambda->getNDaughters() != 2) {
80 B2WARNING("goodBelleLambda is only defined for a particle with two daughters");
81 return 0.;
82 }
83 const Particle* d0 = Lambda->getDaughter(0);
84 const Particle* d1 = Lambda->getDaughter(1);
85 if ((d0->getCharge() == 0) || (d1->getCharge() == 0)) {
86 B2WARNING("goodBelleLambda is only defined for a particle with charged daughters");
87 return 0.;
88 }
89 if (std::abs(Lambda->getPDGCode()) != Const::Lambda.getPDGCode()) {
90 B2WARNING("goodBelleLambda is being applied to a candidate with PDG " << Lambda->getPDGCode());
91 }
92
93 if (Lambda->hasExtraInfo("goodLambda"))
94 return Lambda->getExtraInfo("goodLambda");
95
96 double p = particleP(Lambda);
97 double dr = std::min(std::abs(particleDRho(d0)), std::abs(particleDRho(d1)));
98 double zdist = v0DaughterZ0Diff(Lambda);
99 double dphi = acos(cosAngleBetweenMomentumAndVertexVectorInXYPlane(Lambda));
100 // Flight distance of Lambda0 in xy plane
101 double fl = particleDRho(Lambda);
102
103 // goodBelleLambda == 1 (optimized for proton PID > 0.6)
104 bool high1 = p >= 1.5 && std::abs(zdist) < 12.9 && dr > 0.008 && dphi < 0.09 && fl > 0.22;
105 bool mid1 = p >= 0.5 && p < 1.5 && std::abs(zdist) < 9.8 && dr > 0.01 && dphi < 0.18 && fl > 0.16;
106 bool low1 = p < 0.5 && std::abs(zdist) < 2.4 && dr > 0.027 && dphi < 1.2 && fl > 0.11;
107
108 // goodBelleLambda == 2 (optimized without PID selection)
109 bool high2 = p >= 1.5 && std::abs(zdist) < 7.7 && dr > 0.018 && dphi < 0.07 && fl > 0.35;
110 bool mid2 = p >= 0.5 && p < 1.5 && std::abs(zdist) < 2.1 && dr > 0.033 && dphi < 0.10 && fl > 0.24;
111 bool low2 = p < 0.5 && std::abs(zdist) < 1.9 && dr > 0.059 && dphi < 0.6 && fl > 0.17;
112
113 if (low2 || mid2 || high2) {
114 return 2.0;
115 } else if (low1 || mid1 || high1) {
116 return 1.0;
117 } else {
118 return 0.0;
119 }
120 }
121
122
123 bool isGoodBelleGamma(int region, double energy)
124 {
125 bool goodGammaRegion1, goodGammaRegion2, goodGammaRegion3;
126 goodGammaRegion1 = region == 1 && energy > 0.100;
127 goodGammaRegion2 = region == 2 && energy > 0.050;
128 goodGammaRegion3 = region == 3 && energy > 0.150;
129
130 return goodGammaRegion1 || goodGammaRegion2 || goodGammaRegion3;
131 }
132
133 bool goodBelleGamma(const Particle* particle)
134 {
135 double energy = eclClusterE(particle);
136 int region = eclClusterDetectionRegion(particle);
137
138 return isGoodBelleGamma(region, energy);
139 }
140
141 BelleTrkExtra* getBelleTrkExtraInfoFromParticle(Particle const* particle)
142 {
143 const Track* track = particle->getTrack();
144 if (!track) {
145 return nullptr;
146 }
147 auto belleTrkExtra = track->getRelatedTo<BelleTrkExtra>();
148 if (!belleTrkExtra) {
149 return nullptr;
150 }
151 return belleTrkExtra;
152 }
153
154 double BelleFirstCDCHitX(const Particle* particle)
155 {
156 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
157 if (!belleTrkExtra) {
158 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
159 return Const::doubleNaN;
160 }
161 return belleTrkExtra->getTrackFirstX();
162 }
163
164 double BelleFirstCDCHitY(const Particle* particle)
165 {
166 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
167 if (!belleTrkExtra) {
168 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
169 return Const::doubleNaN;
170 }
171 return belleTrkExtra->getTrackFirstY();
172 }
173
174 double BelleFirstCDCHitZ(const Particle* particle)
175 {
176 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
177 if (!belleTrkExtra) {
178 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
179 return Const::doubleNaN;
180 }
181 return belleTrkExtra->getTrackFirstZ();
182 }
183
184 double BelleLastCDCHitX(const Particle* particle)
185 {
186 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
187 if (!belleTrkExtra) {
188 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
189 return Const::doubleNaN;
190 }
191 return belleTrkExtra->getTrackLastX();
192 }
193
194 double BelleLastCDCHitY(const Particle* particle)
195 {
196 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
197 if (!belleTrkExtra) {
198 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
199 return Const::doubleNaN;
200 }
201 return belleTrkExtra->getTrackLastY();
202 }
203
204 double BelleLastCDCHitZ(const Particle* particle)
205 {
206 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
207 if (!belleTrkExtra) {
208 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
209 return Const::doubleNaN;
210 }
211 return belleTrkExtra->getTrackLastZ();
212 }
213
214 double BellePi0InvariantMassSignificance(const Particle* particle)
215 {
216 TMatrixFSym covarianceMatrix(Particle::c_DimMomentum);
217 for (auto daughter : particle->getDaughters()) {
218 covarianceMatrix += daughter->getMomentumErrorMatrix();
219 }
220
221 TVectorF jacobian(Particle::c_DimMomentum);
222 jacobian[0] = -1.0 * particle->getPx() / particle->getMass();
223 jacobian[1] = -1.0 * particle->getPy() / particle->getMass();
224 jacobian[2] = -1.0 * particle->getPz() / particle->getMass();
225 jacobian[3] = 1.0 * particle->getEnergy() / particle->getMass();
226
227 double massErrSquared = jacobian * (covarianceMatrix * jacobian);
228
229 if (massErrSquared < 0.0)
230 return Const::doubleNaN;
231
232 double invMass = particleInvariantMassFromDaughters(particle);
233 double nomMass = particle->getPDGMass();
234
235 return (invMass - nomMass) / sqrt(massErrSquared);
236 }
237
238 double BelleTof(const Particle* particle)
239 {
240 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
241 if (!belleTrkExtra) {
242 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
243 return Const::doubleNaN;
244 }
245 return belleTrkExtra->getTof();
246 }
247
248 double BelleTofQuality(const Particle* particle)
249 {
250 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
251 if (!belleTrkExtra) {
252 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
253 return Const::doubleNaN;
254 }
255 return belleTrkExtra->getTofQuality();
256 }
257
258 double BelleTofSigma(const Particle* particle)
259 {
260 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
261 if (!belleTrkExtra) {
262 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
263 return Const::doubleNaN;
264 }
265 return belleTrkExtra->getTofSigma();
266 }
267
268 double BellePathLength(const Particle* particle)
269 {
270 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
271 if (!belleTrkExtra) {
272 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
273 return Const::doubleNaN;
274 }
275 return belleTrkExtra->getPathLength();
276 }
277
278 double BelleTofMass(const Particle* particle)
279 {
280 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
281 if (!belleTrkExtra) {
282 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
283 return Const::doubleNaN;
284 }
285 double time = belleTrkExtra->getTof();
286 double length = belleTrkExtra->getPathLength();
287 double p = particle->getP(); //3-momentum
288 double tofbeta = length / time / Const::speedOfLight;
289 double tofmass = p * sqrt(1. / (tofbeta * tofbeta) - 1.); //(GeV)
290
291 return tofmass;
292 }
293
294 double BelledEdx(const Particle* particle)
295 {
296 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
297 if (!belleTrkExtra) {
298 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
299 return Const::doubleNaN;
300 }
301 return belleTrkExtra->getdEdx();
302 }
303
304 double BelledEdxQuality(const Particle* particle)
305 {
306 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
307 if (!belleTrkExtra) {
308 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
309 return Const::doubleNaN;
310 }
311 return belleTrkExtra->getdEdxQuality();
312 }
313
314 double BelleACCnPe(const Particle* particle)
315 {
316 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
317 if (!belleTrkExtra) {
318 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
319 return Const::doubleNaN;
320 }
321 return belleTrkExtra->getACCPe();
322 }
323
324 double BelleACCQuality(const Particle* particle)
325 {
326 auto belleTrkExtra = getBelleTrkExtraInfoFromParticle(particle);
327 if (!belleTrkExtra) {
328 B2WARNING("Cannot find BelleTrkExtra, did you forget to enable BelleTrkExtra during the conversion?");
329 return Const::doubleNaN;
330 }
331 return belleTrkExtra->getACCQuality();
332 }
333
334
335
336 VARIABLE_GROUP("Belle Variables");
337
338 REGISTER_VARIABLE("goodBelleKshort", goodBelleKshort, R"DOC(
339[Legacy] GoodKs Returns true if a :math:`K_{S}^0\to\pi\pi` candidate passes the Belle algorithm:
340a momentum-binned selection including requirements on impact parameter of, and
341angle between the daughter pions as well as separation from the vertex and
342flight distance in the transverse plane.
343
344.. seealso:: `BELLE2-NOTE-PH-2018-017 <https://docs.belle2.org/record/957>`_
345)DOC");
346
347 REGISTER_VARIABLE("goodBelleLambda", goodBelleLambda, R"DOC(
348[Legacy] Returns 2.0, 1.0, 0.0 as an indication of goodness of :math:`\Lambda^0` candidates,
349based on:
350
351* The distance of the two daughter tracks at their interception at z axis,
352* the minimum distance of the daughter tracks and the IP in xy plane,
353* the difference of the azimuthal angle of the vertex vector and the momentum vector,
354* and the flight distance of the Lambda0 candidates in xy plane.
355
356It reproduces the ``goodLambda()`` function in Belle.
357
358``goodBelleLambda`` selection 1 (selected with: ``goodBelleLambda>0``) maximizes the signal significance after applying
359``atcPIDBelle(4,2) > 0.6``, while ``goodBelleLambda`` selection 2 (``goodBelleLambda>1``) is tighter and maximizes the signal
360significance of a :math:`\Lambda^0` sample without any proton PID cut. However, it might still be beneficial to apply a proton PID
361cut on top of it. Which combination of proton PID cut and ``goodBelleLambda`` selection scenario is ideal, is probably
362analysis-dependent.
363
364.. warning:: ``goodBelleLambda`` is not optimized or tested on Belle II data.
365
366See Also:
367 * `BN-684`_ Lambda selection at Belle. K F Chen et al.
368 * The ``FindLambda`` class can be found at ``/belle_legacy/findLambda/findLambda.h``
369
370.. _BN-684: https://belle.kek.jp/secured/belle_note/gn684/bn684.ps.gz
371
372)DOC");
373
374 REGISTER_VARIABLE("goodBelleGamma", goodBelleGamma, R"DOC(
375[Legacy] Returns 1.0 if the photon candidate passes the simple region dependent
376energy selection for Belle data and MC.
377
378.. math::
379
380 E > 50 \textrm{ MeV; barrel}\\
381 E > 100 \textrm{ MeV; forward endcap}\\
382 E > 150 \textrm{ MeV; backward endcap}
383)DOC");
384
385 REGISTER_VARIABLE("BelleFirstCDCHitX", BelleFirstCDCHitX, R"DOC(
386[Legacy] Returns x component of starting point of the track near the 1st SVD or CDC hit for SVD1 data (exp. 7 - 27) or the 1st CDC hit for SVD2 data (from exp. 31). (Belle only, originally stored in mdst_trk_fit.)
387
388)DOC","cm");
389
390 REGISTER_VARIABLE("BelleFirstCDCHitY", BelleFirstCDCHitY, R"DOC(
391[Legacy] Returns y component of starting point of the track near the 1st SVD or CDC hit for SVD1 data (exp. 7 - 27) or the 1st CDC hit for SVD2 data (from exp. 31). (Belle only, originally stored in mdst_trk_fit.)
392
393)DOC","cm");
394
395 REGISTER_VARIABLE("BelleFirstCDCHitZ", BelleFirstCDCHitZ, R"DOC(
396[Legacy] Returns z component of starting point of the track near the 1st SVD or CDC hit for SVD1 data (exp. 7 - 27) or the 1st CDC hit for SVD2 data (from exp. 31). (Belle only, originally stored in mdst_trk_fit.)
397
398)DOC","cm");
399
400 REGISTER_VARIABLE("BelleLastCDCHitX", BelleLastCDCHitX, R"DOC(
401[Legacy] Returns x component of end point of the track near the last CDC hit. (Belle only, originally stored in mdst_trk_fit.)
402
403)DOC","cm");
404
405 REGISTER_VARIABLE("BelleLastCDCHitY", BelleLastCDCHitY, R"DOC(
406[Legacy] Returns y component of end point of the track near the last CDC hit. (Belle only, originally stored in mdst_trk_fit.)
407
408)DOC","cm");
409
410 REGISTER_VARIABLE("BelleLastCDCHitZ", BelleLastCDCHitZ, R"DOC(
411[Legacy] Returns z component of end point of the track near the last CDC hit. (Belle only, originally stored in mdst_trk_fit.)
412
413)DOC","cm");
414
415 REGISTER_VARIABLE("BellePi0SigM", BellePi0InvariantMassSignificance, R"DOC(
416[Legacy] Returns the significance of the :math:`\pi^0` mass used in the FEI for B2BII.
417The significance is calculated as the difference between the reconstructed and the nominal mass divided by the mass uncertainty:
418
419.. math::
420 \frac{m_{\gamma\gamma} - m_{\pi^0}^\textrm{PDG}}{\sigma_{m_{\gamma\gamma}}}
421
422Since the :math:`\pi^0`'s covariance matrix for B2BII is empty, the latter is calculated using the photon daughters' covariance matrices.
423 )DOC");
424
425 REGISTER_VARIABLE("BelleTof", BelleTof, R"DOC(
426[Legacy] Returns the time of flight of a track. (Belle only).
427
428)DOC", "ns");
429
430 REGISTER_VARIABLE("BelleTofQuality", BelleTofQuality, R"DOC(
431[Legacy] Returns the quality flag of the time of flight of a track. Original description from the panther table: 0 if consistency between z of hit by charge Q and corrected times, 1 if zhit from Q NOT consistent with zhit from and correct times. (Belle only).
432)DOC");
433
434 REGISTER_VARIABLE("BelleTofSigma", BelleTofSigma, R"DOC(
435[Legacy] Returns the expected resolution on the time of flight of a track. (Belle only).
436
437)DOC", "ns");
438
439 REGISTER_VARIABLE("BellePathLength", BellePathLength, R"DOC(
440[Legacy] Returns the track path length. This is defined from the closest point to the z-axis up to TOF counter. (Belle only).
441
442)DOC", "cm");
443
444 REGISTER_VARIABLE("BelleTofMass", BelleTofMass, R"DOC(
445[Legacy] Returns the TOF mass calculated from the time of flight and path length. (Belle only).
446)DOC", "GeV/:math:`\\text{c}^2`");
447
448 REGISTER_VARIABLE("BelledEdx", BelledEdx, R"DOC(
449[Legacy] Returns the dE/dx measured in the CDC. (Belle only).
450
451)DOC", "keV/cm");
452
453 REGISTER_VARIABLE("BelledEdxQuality", BelledEdxQuality, R"DOC(
454[Legacy] Returns the quality flag of the dE/dx measured in the CDC. Sadly no information about the code meaning is given in the original panther tables. (Belle only).
455)DOC");
456
457 REGISTER_VARIABLE("BelleACCnPe", BelleACCnPe, R"DOC(
458[Legacy] Returns the number of photo-electrons associated to the track in the ACC. (Belle only).
459)DOC");
460
461 REGISTER_VARIABLE("BelleACCQuality", BelleACCQuality, R"DOC(
462[Legacy] Returns the ACC quality flag. Original definition in the panther tables: if 0 normal, if 1 the track is out of ACC acceptance. (Belle only).
463)DOC");
464
465
466 // this is defined in ECLVariables.{h,cc}
467 REGISTER_VARIABLE("clusterBelleQuality", eclClusterDeltaL, R"DOC(
468[Legacy] Returns ECL cluster's quality indicating a good cluster in GSIM (stored in deltaL of ECL cluster object).
469Belle analysis typically used clusters with quality == 0 in their :math:`E_{\text{extra ECL}}` (Belle only).
470)DOC");
471 }
473}
int getPDGCode() const
PDG code.
Definition Const.h:473
static const ParticleType Lambda
Lambda particle.
Definition Const.h:679
static const double speedOfLight
[cm/ns]
Definition Const.h:695
static const ParticleType Kshort
K^0_S particle.
Definition Const.h:677
static const double doubleNaN
quiet_NaN
Definition Const.h:703
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Abstract base class for different kinds of events.