Belle II Software light-2505-deimos
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 Belle2::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 auto trackWithWeight = cluster->getRelatedFromWithWeight<Track>();
288 if (!trackWithWeight.first)
289 return Const::doubleNaN;
290 return 1. / trackWithWeight.second;
291 }
292
293 VARIABLE_GROUP("KLM Cluster and KlongID");
294
295 REGISTER_VARIABLE("klmClusterKlId", klmClusterKlId,
296 "Returns the KlId classifier output associated to the KLMCluster.");
297 REGISTER_VARIABLE("klmClusterBelleTrackFlag", klmClusterBelleTrackFlag,
298 "Returns the Belle-style Track flag.");
299 REGISTER_VARIABLE("klmClusterBelleECLFlag", klmClusterBelleECLFlag,
300 "Returns the Belle-style ECL flag.");
301 REGISTER_VARIABLE("klmClusterTiming", klmClusterTiming, R"DOC(
302Returns the timing information of the associated KLMCluster.
303
304)DOC","ns");
305 REGISTER_VARIABLE("klmClusterPositionX", klmClusterPositionX, R"DOC(
306Returns the :math:`x` position of the associated KLMCluster.
307
308)DOC","cm");
309 REGISTER_VARIABLE("klmClusterPositionY", klmClusterPositionY, R"DOC(
310Returns the :math:`y` position of the associated KLMCluster.
311
312)DOC","cm");
313 REGISTER_VARIABLE("klmClusterPositionZ", klmClusterPositionZ, R"DOC(
314Returns the :math:`z` position of the associated KLMCluster.
315
316)DOC","cm");
317 REGISTER_VARIABLE("klmClusterInnermostLayer", klmClusterInnermostLayer,
318 "Returns the number of the innermost KLM layer with a 2-dimensional hit of the associated KLMCluster.");
319 REGISTER_VARIABLE("klmClusterLayers", klmClusterLayers,
320 "Returns the number of KLM layers with 2-dimensional hits of the associated KLMCluster.");
321 REGISTER_VARIABLE("klmClusterEnergy", klmClusterEnergy, R"DOC(
322Returns the energy of the associated KLMCluster.
323
324.. warning::
325 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`
326 (: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`).
327 It should be used with caution, and may not be physically meaningful, especially for :math:`n` candidates.
328
329)DOC","GeV");
330 REGISTER_VARIABLE("klmClusterMomentum", klmClusterMomentum, R"DOC(
331Returns the momentum magnitude of the associated KLMCluster.
332
333.. warning::
334 This variable returns an approximation of the momentum, since it is proportional to :b2:var:`klmClusterLayers`
335 (: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`).
336 It should be used with caution, and may not be physically meaningful.
337
338)DOC","GeV/c");
339 REGISTER_VARIABLE("klmClusterIsBKLM", klmClusterIsBKLM,
340 "Returns 1 if the associated KLMCluster is in barrel KLM.");
341 REGISTER_VARIABLE("klmClusterIsEKLM", klmClusterIsEKLM,
342 "Returns 1 if the associated KLMCluster is in endcap KLM.");
343 REGISTER_VARIABLE("klmClusterIsForwardEKLM", klmClusterIsForwardEKLM,
344 "Returns 1 if the associated KLMCluster is in forward endcap KLM.");
345 REGISTER_VARIABLE("klmClusterIsBackwardEKLM", klmClusterIsBackwardEKLM,
346 "Returns 1 if the associated KLMCluster is in backward endcap KLM.");
347 REGISTER_VARIABLE("klmClusterTheta", klmClusterTheta, R"DOC(
348Returns the polar (:math:`\theta`) angle of the associated KLMCluster.
349
350)DOC","rad");
351 REGISTER_VARIABLE("klmClusterPhi", klmClusterPhi, R"DOC(
352Returns the azimuthal (:math:`\phi`) angle of the associated KLMCluster.
353
354)DOC","rad");
355 REGISTER_VARIABLE("maximumKLMAngleCMS", maximumKLMAngleCMS ,
356 "Returns the maximum angle in the CMS frame between the Particle and all KLMClusters in the event.\n\n","rad");
357 REGISTER_VARIABLE("nKLMClusterTrackMatches", nKLMClusterTrackMatches, R"DOC(
358Returns 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.
359)DOC");
360 REGISTER_VARIABLE("nMatchedKLMClusters", nMatchedKLMClusters, R"DOC(
361 Returns the number of KLMClusters matched to the particle. It only works for
362 Particles created either from Tracks or from ECLCluster, while it returns NaN
363 for :math:`K_{L}^0` or :math:`n` candidates originating from KLMClusters.
364 )DOC");
365 REGISTER_VARIABLE("nKLMClusterECLClusterMatches", nKLMClusterECLClusterMatches, R"DOC(
366 Returns the number of ECLClusters matched to the KLMCluster associated to this Particle.
367 )DOC");
368 REGISTER_VARIABLE("klmClusterTrackDistance", klmClusterTrackDistance,
369 "Returns the distance between the Track and the KLMCluster associated to this Particle. This variable returns NaN if there is no Track-to-KLMCluster relationship.\n\n",
370 "cm");
371
372}
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.