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