Belle II Software  light-2212-foldex
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 /* Belle 2 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 using namespace std;
28 
29 namespace Belle2::Variable {
30 
31  double klmClusterKlId(const Particle* particle)
32  {
33  const KLMCluster* cluster = particle->getKLMCluster();
34  if (!cluster) {
35  return std::numeric_limits<double>::quiet_NaN();
36  }
37  const KlId* klid = cluster->getRelatedTo<KlId>();
38  if (!klid) {
39  return std::numeric_limits<double>::quiet_NaN();
40  }
41  return klid->getKlId();
42  }
43 
44  int klmClusterBelleTrackFlag(const Particle* particle)
45  {
46  const float angle = 0.24;
47  const KLMCluster* cluster = particle->getKLMCluster();
48  if (!cluster) {
49  return 0;
50  }
51  const B2Vector3D& pos = cluster->getClusterPosition();
52  StoreArray<TrackFitResult> tracks;
53  for (const TrackFitResult& track : tracks) {
54  const B2Vector3D& trackPos = track.getPosition();
55  if (trackPos.Angle(pos) < angle) {
56  return 1;
57  }
58  }
59  return 0;
60  }
61 
62  int klmClusterBelleECLFlag(const Particle* particle)
63  {
64  const float angle = 0.24;
65  const KLMCluster* klmCluster = particle->getKLMCluster();
66  if (!klmCluster) {
67  return 0;
68  }
69  const B2Vector3D& klmClusterPos = klmCluster->getClusterPosition();
70  StoreArray<ECLCluster> eclClusters;
71  for (const ECLCluster& eclCluster : eclClusters) {
72  const B2Vector3D& eclClusterPos = eclCluster.getClusterPosition();
73  if (eclClusterPos.Angle(klmClusterPos) < angle) {
74  return 1;
75  }
76  }
77  return 0;
78  }
79 
80  double klmClusterTiming(const Particle* particle)
81  {
82  const KLMCluster* cluster = particle->getKLMCluster();
83  if (!cluster) {
84  return std::numeric_limits<double>::quiet_NaN();
85  }
86  return cluster->getTime();
87  }
88 
89 
90  double klmClusterPositionX(const Particle* particle)
91  {
92  const KLMCluster* cluster = particle->getKLMCluster();
93  if (!cluster) {
94  return std::numeric_limits<double>::quiet_NaN();
95  }
96  return cluster->getClusterPosition().X();
97  }
98 
99 
100  double klmClusterPositionY(const Particle* particle)
101  {
102  const KLMCluster* cluster = particle->getKLMCluster();
103  if (!cluster) {
104  return std::numeric_limits<double>::quiet_NaN();
105  }
106  return cluster->getClusterPosition().Y();
107  }
108 
109 
110  double klmClusterPositionZ(const Particle* particle)
111  {
112  const KLMCluster* cluster = particle->getKLMCluster();
113  if (!cluster) {
114  return std::numeric_limits<double>::quiet_NaN();
115  }
116  return cluster->getClusterPosition().Z();
117  }
118 
119 
120  double klmClusterInnermostLayer(const Particle* particle)
121  {
122  const KLMCluster* cluster = particle->getKLMCluster();
123  if (!cluster) {
124  return std::numeric_limits<double>::quiet_NaN();
125  }
126  return cluster->getInnermostLayer();
127  }
128 
129 
130  double klmClusterLayers(const Particle* particle)
131  {
132  const KLMCluster* cluster = particle->getKLMCluster();
133  if (!cluster) {
134  return std::numeric_limits<double>::quiet_NaN();
135  }
136  return cluster->getLayers();
137  }
138 
139  double klmClusterEnergy(const Particle* particle)
140  {
141  const KLMCluster* cluster = particle->getKLMCluster();
142  if (!cluster) {
143  return std::numeric_limits<double>::quiet_NaN();
144  }
145  return cluster->getEnergy();
146  }
147 
148  double klmClusterMomentum(const Particle* particle)
149  {
150  const KLMCluster* cluster = particle->getKLMCluster();
151  if (!cluster) {
152  return std::numeric_limits<double>::quiet_NaN();
153  }
154  return cluster->getMomentumMag();
155  }
156 
157  double klmClusterIsBKLM(const Particle* particle)
158  {
159  const KLMCluster* cluster = particle->getKLMCluster();
160  if (!cluster) {
161  return std::numeric_limits<double>::quiet_NaN();
162  }
163  float clusterZ = cluster->getClusterPosition().Z();
164  if ((clusterZ > -180) && (clusterZ < 275)) {
165  return 1;
166  }
167  return 0;
168  }
169 
170  double klmClusterIsEKLM(const Particle* particle)
171  {
172  const KLMCluster* cluster = particle->getKLMCluster();
173  if (!cluster) {
174  return std::numeric_limits<double>::quiet_NaN();
175  }
176  float clusterZ = cluster->getClusterPosition().Z();
177  if ((clusterZ < -180) || (clusterZ > 275)) {
178  return 1;
179  }
180  return 0;
181  }
182 
183  double klmClusterIsForwardEKLM(const Particle* particle)
184  {
185  const KLMCluster* cluster = particle->getKLMCluster();
186  if (!cluster) {
187  return std::numeric_limits<double>::quiet_NaN();
188  }
189  float clusterZ = cluster->getClusterPosition().Z();
190  if (clusterZ > 275) {
191  return 1;
192  }
193  return 0;
194  }
195 
196  double klmClusterIsBackwardEKLM(const Particle* particle)
197  {
198  const KLMCluster* cluster = particle->getKLMCluster();
199  if (!cluster) {
200  return std::numeric_limits<double>::quiet_NaN();
201  }
202  float clusterZ = cluster->getClusterPosition().Z();
203  if (clusterZ < -180) {
204  return 1;
205  }
206  return 0;
207  }
208 
209  double klmClusterTheta(const Particle* particle)
210  {
211  const KLMCluster* cluster = particle->getKLMCluster();
212  if (!cluster) {
213  return std::numeric_limits<double>::quiet_NaN();
214  }
215  return cluster->getClusterPosition().Theta();
216  }
217 
218  double klmClusterPhi(const Particle* particle)
219  {
220  const KLMCluster* cluster = particle->getKLMCluster();
221  if (!cluster) {
222  return std::numeric_limits<double>::quiet_NaN();
223  }
224  return cluster->getClusterPosition().Phi();
225  }
226 
227  double maximumKLMAngleCMS(const Particle* particle)
228  {
229  // check there actually are KLM clusters in the event
230  StoreArray<KLMCluster> clusters;
231  if (clusters.getEntries() == 0) return std::numeric_limits<double>::quiet_NaN();
232 
233  // get the input particle's vector momentum in the CMS frame
234  PCmsLabTransform T;
235  const B2Vector3D pCms = (T.rotateLabToCms() * particle->get4Vector()).Vect();
236 
237  // find the KLM cluster with the largest angle
238  double maxAngle = 0.0;
239  for (int iKLM = 0; iKLM < clusters.getEntries(); iKLM++) {
240  const B2Vector3D clusterMomentumCms = (T.rotateLabToCms() * clusters[iKLM]->getMomentum()).Vect();
241  double angle = pCms.Angle(clusterMomentumCms);
242  if (angle > maxAngle) maxAngle = angle;
243  }
244  return maxAngle;
245  }
246 
247  double nKLMClusterTrackMatches(const Particle* particle)
248  {
249  const KLMCluster* cluster = particle->getKLMCluster();
250  if (!cluster)
251  return std::numeric_limits<double>::quiet_NaN();
252  size_t out = cluster->getRelationsFrom<Track>().size();
253  return double(out);
254  }
255 
256  double nMatchedKLMClusters(const Particle* particle)
257  {
258  Belle2::Particle::EParticleSourceObject particleSource = particle->getParticleSource();
259  if (particleSource == Particle::EParticleSourceObject::c_Track) {
260  return particle->getTrack()->getRelationsTo<KLMCluster>().size();
261  } else if (particleSource == Particle::EParticleSourceObject::c_ECLCluster) {
262  return particle->getECLCluster()->getRelationsTo<KLMCluster>().size();
263  } else {
264  return std::numeric_limits<double>::quiet_NaN();
265  }
266  }
267 
268  double nKLMClusterECLClusterMatches(const Particle* particle)
269  {
270  const KLMCluster* cluster = particle->getKLMCluster();
271  if (!cluster)
272  return std::numeric_limits<double>::quiet_NaN();
273  size_t out = cluster->getRelationsFrom<ECLCluster>().size();
274  return double(out);
275  }
276 
277  double klmClusterTrackDistance(const Particle* particle)
278  {
279  const KLMCluster* cluster = particle->getKLMCluster();
280  if (!cluster)
281  return std::numeric_limits<double>::quiet_NaN();
282  auto trackWithWeight = cluster->getRelatedFromWithWeight<Track>();
283  if (!trackWithWeight.first)
284  return std::numeric_limits<double>::quiet_NaN();
285  return 1. / trackWithWeight.second;
286  }
287 
288  VARIABLE_GROUP("KLM Cluster and KlongID");
289 
290  REGISTER_VARIABLE("klmClusterKlId", klmClusterKlId,
291  "Returns the KlId classifier output associated to the KLMCluster.");
292  REGISTER_VARIABLE("klmClusterBelleTrackFlag", klmClusterBelleTrackFlag,
293  "Returns the Belle-style Track flag.");
294  REGISTER_VARIABLE("klmClusterBelleECLFlag", klmClusterBelleECLFlag,
295  "Returns the Belle-style ECL flag.");
296  REGISTER_VARIABLE("klmClusterTiming", klmClusterTiming, R"DOC(
297 Returns the timing information of the associated KLMCluster.
298 
299 .. warning::
300  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.
301  It is recommended to not use it in any case until the KLM is calibrated, even for simulated data.
302 
303 )DOC","ns");
304  REGISTER_VARIABLE("klmClusterPositionX", klmClusterPositionX, R"DOC(
305 Returns the :math:`x` position of the associated KLMCluster.
306 )DOC","cm");
307  REGISTER_VARIABLE("klmClusterPositionY", klmClusterPositionY, R"DOC(
308 Returns the :math:`y` position of the associated KLMCluster.
309 )DOC","cm");
310  REGISTER_VARIABLE("klmClusterPositionZ", klmClusterPositionZ, R"DOC(
311 Returns the :math:`z` position of the associated KLMCluster.
312 )DOC","cm");
313  REGISTER_VARIABLE("klmClusterInnermostLayer", klmClusterInnermostLayer,
314  "Returns the number of the innermost KLM layer with a 2-dimensional hit of the associated KLMCluster.");
315  REGISTER_VARIABLE("klmClusterLayers", klmClusterLayers,
316  "Returns the number of KLM layers with 2-dimensional hits of the associated KLMCluster.");
317  REGISTER_VARIABLE("klmClusterEnergy", klmClusterEnergy, R"DOC(
318 Returns the energy of the associated KLMCluster.
319 
320 .. warning::
321  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`
322  (: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`).
323  It should be used with caution, and may not be physically meaningful, especially for :math:`n` candidates.
324 
325 )DOC","GeV");
326  REGISTER_VARIABLE("klmClusterMomentum", klmClusterMomentum, R"DOC(
327 Returns the momentum magnitude of the associated KLMCluster.
328 
329 .. warning::
330  This variable returns an approximation of the momentum, since it is proportional to :b2:var:`klmClusterLayers`
331  (: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`).
332  It should be used with caution, and may not be physically meaningful.
333 
334 )DOC","GeV/c");
335  REGISTER_VARIABLE("klmClusterIsBKLM", klmClusterIsBKLM,
336  "Returns 1 if the associated KLMCluster is in barrel KLM.");
337  REGISTER_VARIABLE("klmClusterIsEKLM", klmClusterIsEKLM,
338  "Returns 1 if the associated KLMCluster is in endcap KLM.");
339  REGISTER_VARIABLE("klmClusterIsForwardEKLM", klmClusterIsForwardEKLM,
340  "Returns 1 if the associated KLMCluster is in forward endcap KLM.");
341  REGISTER_VARIABLE("klmClusterIsBackwardEKLM", klmClusterIsBackwardEKLM,
342  "Returns 1 if the associated KLMCluster is in backward endcap KLM.");
343  REGISTER_VARIABLE("klmClusterTheta", klmClusterTheta, R"DOC(
344 Returns the polar (:math:`\theta`) angle of the associated KLMCluster.
345 )DOC","rad");
346  REGISTER_VARIABLE("klmClusterPhi", klmClusterPhi, R"DOC(
347 Returns the azimuthal (:math:`\phi`) angle of the associated KLMCluster.
348 )DOC","rad");
349  REGISTER_VARIABLE("maximumKLMAngleCMS", maximumKLMAngleCMS ,
350  "Returns the maximum angle in the CMS frame between the Particle and all KLMClusters in the event.","rad");
351  REGISTER_VARIABLE("nKLMClusterTrackMatches", nKLMClusterTrackMatches, R"DOC(
352 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.
353 )DOC");
354  REGISTER_VARIABLE("nMatchedKLMClusters", nMatchedKLMClusters, R"DOC(
355  Returns the number of KLMClusters matched to the particle. It only works for
356  Particles created either from Tracks or from ECLCluster, while it returns NaN
357  for :math:`K_{L}^0` or :math:`n` candidates originating from KLMClusters.
358  )DOC");
359  REGISTER_VARIABLE("nKLMClusterECLClusterMatches", nKLMClusterECLClusterMatches, R"DOC(
360  Returns the number of ECLClusters matched to the KLMCluster associated to this Particle.
361  )DOC");
362  REGISTER_VARIABLE("klmClusterTrackDistance", klmClusterTrackDistance,
363  "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.",
364  "cm");
365 
366 }
DataType Angle(const B2Vector3< DataType > &q) const
The angle w.r.t.
Definition: B2Vector3.h:302
EParticleSourceObject
particle source enumerators
Definition: Particle.h:81
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516