Belle II Software  release-05-01-25
KLMClusterVariables.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2015-2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Anze Zupanc, Torben Ferber, Giacomo De Pietro *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 /* Own header. */
12 #include <analysis/variables/KLMClusterVariables.h>
13 
14 /* Analysis headers. */
15 #include <analysis/utility/PCmsLabTransform.h>
16 #include <analysis/VariableManager/Manager.h>
17 
18 /* Belle 2 headers. */
19 #include <framework/datastore/StoreArray.h>
20 #include <mdst/dataobjects/ECLCluster.h>
21 #include <mdst/dataobjects/KlId.h>
22 #include <mdst/dataobjects/KLMCluster.h>
23 #include <mdst/dataobjects/Track.h>
24 #include <mdst/dataobjects/TrackFitResult.h>
25 
26 using namespace std;
27 
28 namespace Belle2::Variable {
29 
30  double klmClusterKlId(const Particle* particle)
31  {
32  const KLMCluster* cluster = particle->getKLMCluster();
33  if (!cluster) {
34  return std::numeric_limits<double>::quiet_NaN();
35  }
36  const KlId* klid = cluster->getRelatedTo<KlId>();
37  if (!klid) {
38  return std::numeric_limits<double>::quiet_NaN();
39  }
40  return klid->getKlId();
41  }
42 
43  int klmClusterBelleTrackFlag(const Particle* particle)
44  {
45  const float angle = 0.24;
46  const KLMCluster* cluster = particle->getKLMCluster();
47  if (!cluster) {
48  return std::numeric_limits<int>::quiet_NaN();
49  }
50  const TVector3& pos = cluster->getClusterPosition();
51  StoreArray<TrackFitResult> tracks;
52  for (const TrackFitResult& track : tracks) {
53  const TVector3& trackPos = track.getPosition();
54  if (trackPos.Angle(pos) < angle) {
55  return 1;
56  }
57  }
58  return 0;
59  }
60 
61  int klmClusterBelleECLFlag(const Particle* particle)
62  {
63  const float angle = 0.24;
64  const KLMCluster* klmCluster = particle->getKLMCluster();
65  if (!klmCluster) {
66  return std::numeric_limits<int>::quiet_NaN();
67  }
68  const TVector3& klmClusterPos = klmCluster->getClusterPosition();
69  StoreArray<ECLCluster> eclClusters;
70  for (const ECLCluster& eclCluster : eclClusters) {
71  const TVector3& eclClusterPos = eclCluster.getClusterPosition();
72  if (eclClusterPos.Angle(klmClusterPos) < angle) {
73  return 1;
74  }
75  }
76  return 0;
77  }
78 
79  double klmClusterTiming(const Particle* particle)
80  {
81  const KLMCluster* cluster = particle->getKLMCluster();
82  if (!cluster) {
83  return std::numeric_limits<double>::quiet_NaN();
84  }
85  return cluster->getTime();
86  }
87 
88 
89  double klmClusterPositionX(const Particle* particle)
90  {
91  const KLMCluster* cluster = particle->getKLMCluster();
92  if (!cluster) {
93  return std::numeric_limits<double>::quiet_NaN();
94  }
95  return cluster->getClusterPosition().x();
96  }
97 
98 
99  double klmClusterPositionY(const Particle* particle)
100  {
101  const KLMCluster* cluster = particle->getKLMCluster();
102  if (!cluster) {
103  return std::numeric_limits<double>::quiet_NaN();
104  }
105  return cluster->getClusterPosition().y();
106  }
107 
108 
109  double klmClusterPositionZ(const Particle* particle)
110  {
111  const KLMCluster* cluster = particle->getKLMCluster();
112  if (!cluster) {
113  return std::numeric_limits<double>::quiet_NaN();
114  }
115  return cluster->getClusterPosition().z();
116  }
117 
118 
119  double klmClusterInnermostLayer(const Particle* particle)
120  {
121  const KLMCluster* cluster = particle->getKLMCluster();
122  if (!cluster) {
123  return std::numeric_limits<double>::quiet_NaN();
124  }
125  return cluster->getInnermostLayer();
126  }
127 
128 
129  double klmClusterLayers(const Particle* particle)
130  {
131  const KLMCluster* cluster = particle->getKLMCluster();
132  if (!cluster) {
133  return std::numeric_limits<double>::quiet_NaN();
134  }
135  return cluster->getLayers();
136  }
137 
138  double klmClusterEnergy(const Particle* particle)
139  {
140  const KLMCluster* cluster = particle->getKLMCluster();
141  if (!cluster) {
142  return std::numeric_limits<double>::quiet_NaN();
143  }
144  return cluster->getEnergy();
145  }
146 
147  double klmClusterMomentum(const Particle* particle)
148  {
149  const KLMCluster* cluster = particle->getKLMCluster();
150  if (!cluster) {
151  return std::numeric_limits<double>::quiet_NaN();
152  }
153  return cluster->getMomentumMag();
154  }
155 
156  double klmClusterIsBKLM(const Particle* particle)
157  {
158  const KLMCluster* cluster = particle->getKLMCluster();
159  if (!cluster) {
160  return std::numeric_limits<double>::quiet_NaN();
161  }
162  float clusterZ = cluster->getClusterPosition().Z();
163  if ((clusterZ > -180) && (clusterZ < 275)) {
164  return 1;
165  }
166  return 0;
167  }
168 
169  double klmClusterIsEKLM(const Particle* particle)
170  {
171  const KLMCluster* cluster = particle->getKLMCluster();
172  if (!cluster) {
173  return std::numeric_limits<double>::quiet_NaN();
174  }
175  float clusterZ = cluster->getClusterPosition().Z();
176  if ((clusterZ < -180) || (clusterZ > 275)) {
177  return 1;
178  }
179  return 0;
180  }
181 
182  double klmClusterIsForwardEKLM(const Particle* particle)
183  {
184  const KLMCluster* cluster = particle->getKLMCluster();
185  if (!cluster) {
186  return std::numeric_limits<double>::quiet_NaN();
187  }
188  float clusterZ = cluster->getClusterPosition().Z();
189  if (clusterZ > 275) {
190  return 1;
191  }
192  return 0;
193  }
194 
195  double klmClusterIsBackwardEKLM(const Particle* particle)
196  {
197  const KLMCluster* cluster = particle->getKLMCluster();
198  if (!cluster) {
199  return std::numeric_limits<double>::quiet_NaN();
200  }
201  float clusterZ = cluster->getClusterPosition().Z();
202  if (clusterZ < -180) {
203  return 1;
204  }
205  return 0;
206  }
207 
208  double klmClusterTheta(const Particle* particle)
209  {
210  const KLMCluster* cluster = particle->getKLMCluster();
211  if (!cluster) {
212  return std::numeric_limits<double>::quiet_NaN();
213  }
214  return cluster->getClusterPosition().Theta();
215  }
216 
217  double klmClusterPhi(const Particle* particle)
218  {
219  const KLMCluster* cluster = particle->getKLMCluster();
220  if (!cluster) {
221  return std::numeric_limits<double>::quiet_NaN();
222  }
223  return cluster->getClusterPosition().Phi();
224  }
225 
226  double maximumKLMAngleCMS(const Particle* particle)
227  {
228  // check there actually are KLM clusters in the event
229  StoreArray<KLMCluster> clusters;
230  if (clusters.getEntries() == 0) return std::numeric_limits<double>::quiet_NaN();
231 
232  // get the input particle's vector momentum in the CMS frame
233  PCmsLabTransform T;
234  const TVector3 pCms = (T.rotateLabToCms() * particle->get4Vector()).Vect();
235 
236  // find the KLM cluster with the largest angle
237  double maxAngle = 0.0;
238  for (int iKLM = 0; iKLM < clusters.getEntries(); iKLM++) {
239  const TVector3 clusterMomentumCms = (T.rotateLabToCms() * clusters[iKLM]->getMomentum()).Vect();
240  double angle = pCms.Angle(clusterMomentumCms);
241  if (angle > maxAngle) maxAngle = angle;
242  }
243  return maxAngle;
244  }
245 
246  double nKLMClusterTrackMatches(const Particle* particle)
247  {
248  const KLMCluster* cluster = particle->getKLMCluster();
249  if (!cluster)
250  return std::numeric_limits<double>::quiet_NaN();
251  size_t out = cluster->getRelationsFrom<Track>().size();
252  return double(out);
253  }
254 
255  double nMatchedKLMClusters(const Particle* particle)
256  {
257  const Track* track = particle->getTrack();
258  if (!track)
259  return std::numeric_limits<double>::quiet_NaN();
260  size_t out = track->getRelationsTo<KLMCluster>().size();
261  return double(out);
262  }
263 
264  double klmClusterTrackDistance(const Particle* particle)
265  {
266  const KLMCluster* cluster = particle->getKLMCluster();
267  if (!cluster)
268  return std::numeric_limits<double>::quiet_NaN();
269  auto trackWithWeight = cluster->getRelatedFromWithWeight<Track>();
270  if (!trackWithWeight.first)
271  return std::numeric_limits<double>::quiet_NaN();
272  return 1. / trackWithWeight.second;
273  }
274 
275  VARIABLE_GROUP("KLM Cluster and KlongID");
276 
277  REGISTER_VARIABLE("klmClusterKlId", klmClusterKlId,
278  "Returns the KlId classifier output associated to the KLMCluster.");
279  REGISTER_VARIABLE("klmClusterBelleTrackFlag", klmClusterBelleTrackFlag,
280  "Returns the Belle-style Track flag.");
281  REGISTER_VARIABLE("klmClusterBelleECLFlag", klmClusterBelleECLFlag,
282  "Returns the Belle-style ECL flag.");
283  REGISTER_VARIABLE("klmClusterTiming", klmClusterTiming, R"DOC(
284 Returns the timing informationf of the associated KLMCluster.
285 
286 .. warning::
287  Currently the KLM has no time calibration. This leads to a huge discrepancy for the variable :b2:var:`klmClusterTiming` if one compares collisions and simulated data. Moreover, the distribution of the variable on collisions data has a very complicated structure, due to the different timing shifts of different KLM components.
288  It is recommended to not use it in any case until the KLM is calibrated, even for simulated data.
289 
290 )DOC");
291  REGISTER_VARIABLE("klmClusterPositionX", klmClusterPositionX, R"DOC(
292 Returns the :math:`x` position of the associated KLMCluster.
293 )DOC");
294  REGISTER_VARIABLE("klmClusterPositionY", klmClusterPositionY, R"DOC(
295 Returns the :math:`y` position of the associated KLMCluster.
296 )DOC");
297  REGISTER_VARIABLE("klmClusterPositionZ", klmClusterPositionZ, R"DOC(
298 Returns the :math:`z` position of the associated KLMCluster.
299 )DOC");
300  REGISTER_VARIABLE("klmClusterInnermostLayer", klmClusterInnermostLayer,
301  "Returns the number of the innermost KLM layer with a 2-dimensional hit of the associated KLMCluster.");
302  REGISTER_VARIABLE("klmClusterLayers", klmClusterLayers,
303  "Returns the number of KLM layers with 2-dimensional hits of the associated KLMCluster.");
304  REGISTER_VARIABLE("klmClusterEnergy", klmClusterEnergy, R"DOC(
305 Returns the energy of the associated KLMCluster.
306 
307 .. warning::
308  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`
309  (: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`).
310  It should be used with caution, and may not be physically meaningful, especially for :math:`n` candidates.
311 
312 )DOC");
313  REGISTER_VARIABLE("klmClusterMomentum", klmClusterMomentum, R"DOC(
314 Returns the momentum magnitude of the associated KLMCluster.
315 
316 .. warning::
317  This variable returns an approximation of the momentum, since it is proportional to :b2:var:`klmClusterLayers`
318  (: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`).
319  It should be used with caution, and may not be physically meaningful.
320 
321 )DOC");
322  REGISTER_VARIABLE("klmClusterIsBKLM", klmClusterIsBKLM,
323  "Returns 1 if the associated KLMCluster is in barrel KLM.");
324  REGISTER_VARIABLE("klmClusterIsEKLM", klmClusterIsEKLM,
325  "Returns 1 if the associated KLMCluster is in endcap KLM.");
326  REGISTER_VARIABLE("klmClusterIsForwardEKLM", klmClusterIsForwardEKLM,
327  "Returns 1 if the associated KLMCluster is in forward endcap KLM.");
328  REGISTER_VARIABLE("klmClusterIsBackwardEKLM", klmClusterIsBackwardEKLM,
329  "Returns 1 if the associated KLMCluster is in backward endcap KLM.");
330  REGISTER_VARIABLE("klmClusterTheta", klmClusterTheta, R"DOC(
331 Returns the polar (:math:`\theta`) angle of the associated KLMCluster.
332 )DOC");
333  REGISTER_VARIABLE("klmClusterPhi", klmClusterPhi, R"DOC(
334 Returns the azimuthal (:math:`\phi`) angle of the associated KLMCluster.
335 )DOC");
336  REGISTER_VARIABLE("maximumKLMAngleCMS", maximumKLMAngleCMS ,
337  "Returns the maximum angle in the CMS frame between the Particle and all KLMClusters in the event.");
338  REGISTER_VARIABLE("nKLMClusterTrackMatches", nKLMClusterTrackMatches, R"DOC(
339 Returns 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.
340 )DOC");
341  REGISTER_VARIABLE("nMatchedKLMClusters", nMatchedKLMClusters, R"DOC(
342 Returns the number of KLMClusters matched to the Track associated to this Particle. This variable can return only 0 or 1 and return NaN for :math:`K_{L}^0` or :math:`n` candidates originating from KLMClusters with no Tracks associated.
343 )DOC");
344  REGISTER_VARIABLE("klmClusterTrackDistance", klmClusterTrackDistance,
345  "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.");
346 
347 }
KlId
Helper functions for all klid modules to improve readability of the code.
Definition: KlId.h:28