Belle II Software development
KLMClusterVariables.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/KLMClusterVariables.h>
11
12// include VariableManager
13#include <analysis/VariableManager/Manager.h>
14
15/* Analysis headers. */
16#include <analysis/dataobjects/Particle.h>
17#include <analysis/utility/PCmsLabTransform.h>
18
19/* Basf2 headers. */
20#include <framework/datastore/StoreArray.h>
21#include <mdst/dataobjects/ECLCluster.h>
22#include <mdst/dataobjects/KlId.h>
23#include <mdst/dataobjects/KLMCluster.h>
24#include <mdst/dataobjects/Track.h>
25#include <mdst/dataobjects/TrackFitResult.h>
26
27/* ROOT headers. */
28#include <Math/Vector3D.h>
29#include <Math/Vector4D.h>
30#include <Math/VectorUtil.h>
31
32using namespace std;
33
34namespace Belle2::Variable {
35
36 double klmClusterKlId(const Particle* particle)
37 {
38 const KLMCluster* cluster = particle->getKLMCluster();
39 if (!cluster) {
40 return Const::doubleNaN;
41 }
42 const KlId* klid = cluster->getRelatedTo<KlId>();
43 if (!klid) {
44 return Const::doubleNaN;
45 }
46 return klid->getKlId();
47 }
48
49 int klmClusterBelleTrackFlag(const Particle* particle)
50 {
51 const float angle = 0.24;
52 const KLMCluster* cluster = particle->getKLMCluster();
53 if (!cluster) {
54 return 0;
55 }
56 const ROOT::Math::XYZVector& pos = cluster->getClusterPosition();
57 StoreArray<TrackFitResult> tracks;
58 for (const TrackFitResult& track : tracks) {
59 const ROOT::Math::XYZVector& trackPos = track.getPosition();
60 if (ROOT::Math::VectorUtil::Angle(trackPos, pos) < angle) {
61 return 1;
62 }
63 }
64 return 0;
65 }
66
67 int klmClusterBelleECLFlag(const Particle* particle)
68 {
69 const float angle = 0.24;
70 const KLMCluster* klmCluster = particle->getKLMCluster();
71 if (!klmCluster) {
72 return 0;
73 }
74 const ROOT::Math::XYZVector& klmClusterPos = klmCluster->getClusterPosition();
75 StoreArray<ECLCluster> eclClusters;
76 for (const ECLCluster& eclCluster : eclClusters) {
77 const ROOT::Math::XYZVector& eclClusterPos = eclCluster.getClusterPosition();
78 if (ROOT::Math::VectorUtil::Angle(eclClusterPos, klmClusterPos) < angle) {
79 return 1;
80 }
81 }
82 return 0;
83 }
84
85 double klmClusterTiming(const Particle* particle)
86 {
87 const KLMCluster* cluster = particle->getKLMCluster();
88 if (!cluster) {
89 return Const::doubleNaN;
90 }
91 return cluster->getTime();
92 }
93
94
95 double klmClusterPositionX(const Particle* particle)
96 {
97 const KLMCluster* cluster = particle->getKLMCluster();
98 if (!cluster) {
99 return Const::doubleNaN;
100 }
101 return cluster->getClusterPosition().X();
102 }
103
104
105 double klmClusterPositionY(const Particle* particle)
106 {
107 const KLMCluster* cluster = particle->getKLMCluster();
108 if (!cluster) {
109 return Const::doubleNaN;
110 }
111 return cluster->getClusterPosition().Y();
112 }
113
114
115 double klmClusterPositionZ(const Particle* particle)
116 {
117 const KLMCluster* cluster = particle->getKLMCluster();
118 if (!cluster) {
119 return Const::doubleNaN;
120 }
121 return cluster->getClusterPosition().Z();
122 }
123
124
125 double klmClusterInnermostLayer(const Particle* particle)
126 {
127 const KLMCluster* cluster = particle->getKLMCluster();
128 if (!cluster) {
129 return Const::doubleNaN;
130 }
131 return cluster->getInnermostLayer();
132 }
133
134
135 double klmClusterLayers(const Particle* particle)
136 {
137 const KLMCluster* cluster = particle->getKLMCluster();
138 if (!cluster) {
139 return Const::doubleNaN;
140 }
141 return cluster->getLayers();
142 }
143
144 double klmClusterEnergy(const Particle* particle)
145 {
146 const KLMCluster* cluster = particle->getKLMCluster();
147 if (!cluster) {
148 return Const::doubleNaN;
149 }
150 return cluster->getEnergy();
151 }
152
153 double klmClusterMomentum(const Particle* particle)
154 {
155 const KLMCluster* cluster = particle->getKLMCluster();
156 if (!cluster) {
157 return Const::doubleNaN;
158 }
159 return cluster->getMomentumMag();
160 }
161
162 double klmClusterIsBKLM(const Particle* particle)
163 {
164 const KLMCluster* cluster = particle->getKLMCluster();
165 if (!cluster) {
166 return Const::doubleNaN;
167 }
168 float clusterZ = cluster->getClusterPosition().Z();
169 if ((clusterZ > -180) && (clusterZ < 275)) {
170 return 1;
171 }
172 return 0;
173 }
174
175 double klmClusterIsEKLM(const Particle* particle)
176 {
177 const KLMCluster* cluster = particle->getKLMCluster();
178 if (!cluster) {
179 return Const::doubleNaN;
180 }
181 float clusterZ = cluster->getClusterPosition().Z();
182 if ((clusterZ < -180) || (clusterZ > 275)) {
183 return 1;
184 }
185 return 0;
186 }
187
188 double klmClusterIsForwardEKLM(const Particle* particle)
189 {
190 const KLMCluster* cluster = particle->getKLMCluster();
191 if (!cluster) {
192 return Const::doubleNaN;
193 }
194 float clusterZ = cluster->getClusterPosition().Z();
195 if (clusterZ > 275) {
196 return 1;
197 }
198 return 0;
199 }
200
201 double klmClusterIsBackwardEKLM(const Particle* particle)
202 {
203 const KLMCluster* cluster = particle->getKLMCluster();
204 if (!cluster) {
205 return Const::doubleNaN;
206 }
207 float clusterZ = cluster->getClusterPosition().Z();
208 if (clusterZ < -180) {
209 return 1;
210 }
211 return 0;
212 }
213
214 double klmClusterTheta(const Particle* particle)
215 {
216 const KLMCluster* cluster = particle->getKLMCluster();
217 if (!cluster) {
218 return Const::doubleNaN;
219 }
220 return cluster->getClusterPosition().Theta();
221 }
222
223 double klmClusterPhi(const Particle* particle)
224 {
225 const KLMCluster* cluster = particle->getKLMCluster();
226 if (!cluster) {
227 return Const::doubleNaN;
228 }
229 return cluster->getClusterPosition().Phi();
230 }
231
232 double maximumKLMAngleCMS(const Particle* particle)
233 {
234 // check there actually are KLM clusters in the event
235 StoreArray<KLMCluster> clusters;
236 if (clusters.getEntries() == 0) return Const::doubleNaN;
237
238 // get the input particle's vector momentum in the CMS frame
239 PCmsLabTransform T;
240 const ROOT::Math::PxPyPzEVector pCms = T.rotateLabToCms() * particle->get4Vector();
241
242 // find the KLM cluster with the largest angle
243 double maxAngle = 0.0;
244 for (int iKLM = 0; iKLM < clusters.getEntries(); iKLM++) {
245 const ROOT::Math::PxPyPzEVector clusterMomentumCms = T.rotateLabToCms() * clusters[iKLM]->getMomentum();
246 double angle = ROOT::Math::VectorUtil::Angle(pCms, clusterMomentumCms);
247 if (angle > maxAngle) maxAngle = angle;
248 }
249 return maxAngle;
250 }
251
252 double nKLMClusterTrackMatches(const Particle* particle)
253 {
254 const KLMCluster* cluster = particle->getKLMCluster();
255 if (!cluster)
256 return Const::doubleNaN;
257 size_t out = cluster->getRelationsFrom<Track>().size();
258 return double(out);
259 }
260
261 double nMatchedKLMClusters(const Particle* particle)
262 {
263 Particle::EParticleSourceObject particleSource = particle->getParticleSource();
264 if (particleSource == Particle::EParticleSourceObject::c_Track) {
265 return particle->getTrack()->getRelationsTo<KLMCluster>().size();
266 } else if (particleSource == Particle::EParticleSourceObject::c_ECLCluster) {
267 return particle->getECLCluster()->getRelationsTo<KLMCluster>().size();
268 } else {
269 return Const::doubleNaN;
270 }
271 }
272
273 double nKLMClusterECLClusterMatches(const Particle* particle)
274 {
275 const KLMCluster* cluster = particle->getKLMCluster();
276 if (!cluster)
277 return Const::doubleNaN;
278 size_t out = cluster->getRelationsFrom<ECLCluster>().size();
279 return double(out);
280 }
281
282 double klmClusterTrackDistance(const Particle* particle)
283 {
284 const KLMCluster* cluster = particle->getKLMCluster();
285 if (!cluster)
286 return Const::doubleNaN;
287 return cluster->getClusterTrackSeparation();
288 }
289
290 double klmClusterTrackSeparationAngle(const Particle* particle)
291 {
292 const KLMCluster* cluster = particle->getKLMCluster();
293 if (!cluster)
294 return Const::doubleNaN;
295 return cluster->getClusterTrackSeparationAngle();
296 }
297
298 double klmClusterTrackRotationAngle(const Particle* particle)
299 {
300 const KLMCluster* cluster = particle->getKLMCluster();
301 if (!cluster)
302 return Const::doubleNaN;
303 return cluster->getClusterTrackRotationAngle();
304 }
305
306 double klmClusterShapeStdDev1(const Particle* particle)
307 {
308 const KLMCluster* cluster = particle->getKLMCluster();
309 if (!cluster)
310 return Const::doubleNaN;
311 return cluster->getShapeStdDev1();
312 }
313
314 double klmClusterShapeStdDev2(const Particle* particle)
315 {
316 const KLMCluster* cluster = particle->getKLMCluster();
317 if (!cluster)
318 return Const::doubleNaN;
319 return cluster->getShapeStdDev2();
320 }
321
322 double klmClusterShapeStdDev3(const Particle* particle)
323 {
324 const KLMCluster* cluster = particle->getKLMCluster();
325 if (!cluster)
326 return Const::doubleNaN;
327 return cluster->getShapeStdDev3();
328 }
329
330 VARIABLE_GROUP("KLM Cluster and KlongID");
331
332 REGISTER_VARIABLE("klmClusterKlId", klmClusterKlId,
333 "Returns the KlId classifier output associated to the KLMCluster.");
334 REGISTER_VARIABLE("klmClusterBelleTrackFlag", klmClusterBelleTrackFlag,
335 "Returns the Belle-style Track flag.");
336 REGISTER_VARIABLE("klmClusterBelleECLFlag", klmClusterBelleECLFlag,
337 "Returns the Belle-style ECL flag.");
338 REGISTER_VARIABLE("klmClusterTiming", klmClusterTiming, R"DOC(
339Returns the timing information of the associated KLMCluster.
340
341)DOC","ns");
342 REGISTER_VARIABLE("klmClusterPositionX", klmClusterPositionX, R"DOC(
343Returns the :math:`x` position of the associated KLMCluster.
344
345)DOC","cm");
346 REGISTER_VARIABLE("klmClusterPositionY", klmClusterPositionY, R"DOC(
347Returns the :math:`y` position of the associated KLMCluster.
348
349)DOC","cm");
350 REGISTER_VARIABLE("klmClusterPositionZ", klmClusterPositionZ, R"DOC(
351Returns the :math:`z` position of the associated KLMCluster.
352
353)DOC","cm");
354 REGISTER_VARIABLE("klmClusterInnermostLayer", klmClusterInnermostLayer,
355 "Returns the number of the innermost KLM layer with a 2-dimensional hit of the associated KLMCluster.");
356 REGISTER_VARIABLE("klmClusterLayers", klmClusterLayers,
357 "Returns the number of KLM layers with 2-dimensional hits of the associated KLMCluster.");
358 REGISTER_VARIABLE("klmClusterEnergy", klmClusterEnergy, R"DOC(
359Returns the energy of the associated KLMCluster.
360
361.. warning::
362 This variable returns an approximation of the energy: it uses :b2:var:`klmClusterMomentum` as momentum and the hypothesis that the KLMCluster is originated by a :math:`K_{L}^0`
363 (:math:`E_{\text{KLM}} = \sqrt{M_{K^0_L}^2 + p_{\text{KLM}}^2}`, where :math:`E_{\text{KLM}}` is this variable, :math:`M_{K^0_L}` is the :math:`K^0_L` mass and :math:`p_{\text{KLM}}` is :b2:var:`klmClusterMomentum`).
364 It should be used with caution, and may not be physically meaningful, especially for :math:`n` candidates.
365
366)DOC","GeV");
367 REGISTER_VARIABLE("klmClusterMomentum", klmClusterMomentum, R"DOC(
368Returns the momentum magnitude of the associated KLMCluster.
369
370.. warning::
371 This variable returns an approximation of the momentum, since it is proportional to :b2:var:`klmClusterLayers`
372 (:math:`p_{\text{KLM}} = 0.215 \cdot N_{\text{layers}}`, where :math:`p_{\text{KLM}}` is this variable and :math:`N_{\text{layers}}` is :b2:var:`klmClusterLayers`).
373 It should be used with caution, and may not be physically meaningful.
374
375)DOC","GeV/c");
376 REGISTER_VARIABLE("klmClusterIsBKLM", klmClusterIsBKLM,
377 "Returns 1 if the associated KLMCluster is in barrel KLM.");
378 REGISTER_VARIABLE("klmClusterIsEKLM", klmClusterIsEKLM,
379 "Returns 1 if the associated KLMCluster is in endcap KLM.");
380 REGISTER_VARIABLE("klmClusterIsForwardEKLM", klmClusterIsForwardEKLM,
381 "Returns 1 if the associated KLMCluster is in forward endcap KLM.");
382 REGISTER_VARIABLE("klmClusterIsBackwardEKLM", klmClusterIsBackwardEKLM,
383 "Returns 1 if the associated KLMCluster is in backward endcap KLM.");
384 REGISTER_VARIABLE("klmClusterTheta", klmClusterTheta, R"DOC(
385Returns the polar (:math:`\theta`) angle of the associated KLMCluster.
386
387)DOC","rad");
388 REGISTER_VARIABLE("klmClusterPhi", klmClusterPhi, R"DOC(
389Returns the azimuthal (:math:`\phi`) angle of the associated KLMCluster.
390
391)DOC","rad");
392 REGISTER_VARIABLE("maximumKLMAngleCMS", maximumKLMAngleCMS ,
393 "Returns the maximum angle in the CMS frame between the Particle and all KLMClusters in the event.\n\n","rad");
394 REGISTER_VARIABLE("nKLMClusterTrackMatches", nKLMClusterTrackMatches, R"DOC(
395Returns the number of Tracks matched to the KLMCluster associated to this Particle. This variable can return a number greater than 0 for :math:`K_{L}^0` or :math:`n` candidates originating from KLMClusters and returns NaN for Particles with no KLMClusters associated.
396)DOC");
397 REGISTER_VARIABLE("nMatchedKLMClusters", nMatchedKLMClusters, R"DOC(
398 Returns the number of KLMClusters matched to the particle. It only works for
399 Particles created either from Tracks or from ECLCluster, while it returns NaN
400 for :math:`K_{L}^0` or :math:`n` candidates originating from KLMClusters.
401 )DOC");
402 REGISTER_VARIABLE("nKLMClusterECLClusterMatches", nKLMClusterECLClusterMatches, R"DOC(
403 Returns the number of ECLClusters matched to the KLMCluster associated to this Particle.
404 )DOC");
405 REGISTER_VARIABLE("klmClusterTrackDistance", klmClusterTrackDistance,
406 "Returns the distance between KLMCluster associated to this Particle and the closest track. This variable returns NaN if there is no Track-to-KLMCluster relationship.\n\n",
407 "cm");
408 REGISTER_VARIABLE("klmClusterTrackRotationAngle", klmClusterTrackRotationAngle,
409 "Returns the angle between the direction at the IP and at the POCA to the KLMCluster associated to this Particle for the closest track. This variable returns NaN if there is no Track-to-KLMCluster relationship.\n\n",
410 "rad");
411 REGISTER_VARIABLE("klmClusterTrackSeparationAngle", klmClusterTrackSeparationAngle,
412 "Returns the angle between the KLMCluster associated to this Particle and the closest track. This variable returns NaN if there is no Track-to-KLMCluster relationship.\n\n",
413 "rad");
414
415 REGISTER_VARIABLE("klmClusterShapeStdDev1", klmClusterShapeStdDev1,
416 "Returns the std deviation of the 1st axis from a PCA of the KLMCluster associated to this Particle. This variable returns 0 if this KLMCluster contains only one KLMHit2d cluster.\n\n",
417 "cm");
418 REGISTER_VARIABLE("klmClusterShapeStdDev2", klmClusterShapeStdDev2,
419 "Returns the std deviation of the 2nd axis from a PCA of the KLMCluster associated to this Particle. This variable returns 0 if this KLMCluster contains only one KLMHit2d cluster.\n\n",
420 "cm");
421 REGISTER_VARIABLE("klmClusterShapeStdDev3", klmClusterShapeStdDev3,
422 "Returns the std deviation of the 3rd axis from a PCA of the KLMCluster associated to this Particle. This variable returns 0 if this KLMCluster contains only one KLMHit2d cluster.\n\n",
423 "cm");
424}
Belle2::Const::doubleNaN
static const double doubleNaN
quiet_NaN
Definition Const.h:703
Belle2::Particle::EParticleSourceObject
EParticleSourceObject
particle source enumerators
Definition Particle.h:83
std
STL namespace.