10#include <analysis/variables/KLMClusterVariables.h>
13#include <analysis/VariableManager/Manager.h>
16#include <analysis/dataobjects/Particle.h>
17#include <analysis/utility/PCmsLabTransform.h>
18#include <analysis/utility/DetectorSurface.h>
21#include <framework/datastore/StoreArray.h>
22#include <framework/dataobjects/Helix.h>
23#include <framework/gearbox/Const.h>
24#include <framework/geometry/VectorUtil.h>
25#include <mdst/dataobjects/ECLCluster.h>
26#include <mdst/dataobjects/KlId.h>
27#include <mdst/dataobjects/KLMCluster.h>
28#include <mdst/dataobjects/Track.h>
29#include <mdst/dataobjects/TrackFitResult.h>
32#include <Math/Vector3D.h>
33#include <Math/Vector4D.h>
34#include <Math/VectorUtil.h>
41namespace Belle2::Variable {
43 double klmClusterKlId(
const Particle* particle)
45 const KLMCluster* cluster = particle->getKLMCluster();
49 const KlId* klid = cluster->getRelatedTo<KlId>();
53 return klid->getKlId();
56 int klmClusterBelleTrackFlag(
const Particle* particle)
58 const float angle = 0.24;
59 const KLMCluster* cluster = particle->getKLMCluster();
63 const ROOT::Math::XYZVector& pos = cluster->getClusterPosition();
64 StoreArray<TrackFitResult> tracks;
65 for (
const TrackFitResult& track : tracks) {
66 const ROOT::Math::XYZVector& trackPos = track.getPosition();
67 if (ROOT::Math::VectorUtil::Angle(trackPos, pos) < angle) {
74 int klmClusterBelleECLFlag(
const Particle* particle)
76 const float angle = 0.24;
77 const KLMCluster* klmCluster = particle->getKLMCluster();
81 const ROOT::Math::XYZVector& klmClusterPos = klmCluster->getClusterPosition();
82 StoreArray<ECLCluster> eclClusters;
83 for (
const ECLCluster& eclCluster : eclClusters) {
84 const ROOT::Math::XYZVector& eclClusterPos = eclCluster.getClusterPosition();
85 if (ROOT::Math::VectorUtil::Angle(eclClusterPos, klmClusterPos) < angle) {
92 double klmClusterTiming(
const Particle* particle)
94 const KLMCluster* cluster = particle->getKLMCluster();
98 return cluster->getTime();
102 double klmClusterPositionX(
const Particle* particle)
104 const KLMCluster* cluster = particle->getKLMCluster();
108 return cluster->getClusterPosition().X();
112 double klmClusterPositionY(
const Particle* particle)
114 const KLMCluster* cluster = particle->getKLMCluster();
118 return cluster->getClusterPosition().Y();
122 double klmClusterPositionZ(
const Particle* particle)
124 const KLMCluster* cluster = particle->getKLMCluster();
128 return cluster->getClusterPosition().Z();
132 double klmClusterInnermostLayer(
const Particle* particle)
134 const KLMCluster* cluster = particle->getKLMCluster();
138 return cluster->getInnermostLayer();
142 double klmClusterLayers(
const Particle* particle)
144 const KLMCluster* cluster = particle->getKLMCluster();
148 return cluster->getLayers();
151 double klmClusterEnergy(
const Particle* particle)
153 const KLMCluster* cluster = particle->getKLMCluster();
157 return cluster->getEnergy();
160 double klmClusterMomentum(
const Particle* particle)
162 const KLMCluster* cluster = particle->getKLMCluster();
166 return cluster->getMomentumMag();
169 double klmClusterIsBKLM(
const Particle* particle)
171 const KLMCluster* cluster = particle->getKLMCluster();
175 float clusterZ = cluster->getClusterPosition().Z();
176 if ((clusterZ > -180) && (clusterZ < 275)) {
182 double klmClusterIsEKLM(
const Particle* particle)
184 const KLMCluster* cluster = particle->getKLMCluster();
188 float clusterZ = cluster->getClusterPosition().Z();
189 if ((clusterZ < -180) || (clusterZ > 275)) {
195 double klmClusterIsForwardEKLM(
const Particle* particle)
197 const KLMCluster* cluster = particle->getKLMCluster();
201 float clusterZ = cluster->getClusterPosition().Z();
202 if (clusterZ > 275) {
208 double klmClusterIsBackwardEKLM(
const Particle* particle)
210 const KLMCluster* cluster = particle->getKLMCluster();
214 float clusterZ = cluster->getClusterPosition().Z();
215 if (clusterZ < -180) {
221 double klmClusterTheta(
const Particle* particle)
223 const KLMCluster* cluster = particle->getKLMCluster();
227 return cluster->getClusterPosition().Theta();
230 double klmClusterPhi(
const Particle* particle)
232 const KLMCluster* cluster = particle->getKLMCluster();
236 return cluster->getClusterPosition().Phi();
239 double maximumKLMAngleCMS(
const Particle* particle)
242 StoreArray<KLMCluster> clusters;
247 const ROOT::Math::PxPyPzEVector pCms = T.rotateLabToCms() * particle->get4Vector();
250 double maxAngle = 0.0;
251 for (
int iKLM = 0; iKLM < clusters.getEntries(); iKLM++) {
252 const ROOT::Math::PxPyPzEVector clusterMomentumCms = T.rotateLabToCms() * clusters[iKLM]->getMomentum();
253 double angle = ROOT::Math::VectorUtil::Angle(pCms, clusterMomentumCms);
254 if (angle > maxAngle) maxAngle = angle;
259 double nKLMClusterTrackMatches(
const Particle* particle)
261 const KLMCluster* cluster = particle->getKLMCluster();
264 size_t out = cluster->getRelationsFrom<Track>().size();
268 double nMatchedKLMClusters(
const Particle* particle)
271 if (particleSource == Particle::EParticleSourceObject::c_Track) {
272 return particle->getTrack()->getRelationsTo<KLMCluster>().size();
273 }
else if (particleSource == Particle::EParticleSourceObject::c_ECLCluster) {
274 return particle->getECLCluster()->getRelationsTo<KLMCluster>().size();
280 double nKLMClusterECLClusterMatches(
const Particle* particle)
282 const KLMCluster* cluster = particle->getKLMCluster();
285 size_t out = cluster->getRelationsFrom<ECLCluster>().size();
289 double klmClusterTrackDistance(
const Particle* particle)
291 const KLMCluster* cluster = particle->getKLMCluster();
294 return cluster->getClusterTrackSeparation();
297 double ClusterTrackDistance_usingHelixExtrapolate(
const KLMCluster* cluster,
const Belle2::Track* specificTrack =
nullptr)
301 const double r_BKLM = klmBounds.m_rho;
302 const double z_EFWD = klmBounds.m_zfwd;
303 const double z_EBWD = klmBounds.m_zbwd;
307 auto getPositionOnHelix_trackfit = [&](
const TrackFitResult * trackFit) -> ROOT::Math::XYZVector {
308 if (!trackFit)
return vecNaN;
311 const double z0 = trackFit->getZ0();
312 const double tanlambda = trackFit->getTanLambda();
313 const Helix h = trackFit->getHelix();
316 const double arcLength = h.getArcLength2DAtCylindricalR(r_BKLM);
317 const double lHelixRadius = arcLength > 0 ? arcLength : std::numeric_limits<double>::max();
320 const double lFWD = (z_EFWD - z0) / tanlambda > 0 ? (z_EFWD - z0) / tanlambda : std::numeric_limits<double>::max();
323 const double lBWD = (z_EBWD - z0) / tanlambda > 0 ? (z_EBWD - z0) / tanlambda : std::numeric_limits<double>::max();
326 const double l = std::min({lHelixRadius, lFWD, lBWD});
328 ROOT::Math::XYZVector ext_helix = h.getPositionAtArcLength2D(l);
329 double helixExtR_surface = r_BKLM / sin(ext_helix.Theta());
333 helixExtR_surface = z_EFWD / cos(ext_helix.Theta());
334 }
else if (l == lBWD)
336 helixExtR_surface = z_EBWD / cos(ext_helix.Theta());
340 if (!std::isfinite(helixExtR_surface))
return vecNaN;
342 ROOT::Math::XYZVector helixExt_surface_position(0, 0, 0);
343 VectorUtil::setMagThetaPhi(helixExt_surface_position, helixExtR_surface, ext_helix.Theta(), ext_helix.Phi());
345 return helixExt_surface_position;
348 bool isEFWD = (cluster->getClusterPosition().Z() > klmBounds.m_zfwd);
349 bool isEBWD = (cluster->getClusterPosition().Z() < klmBounds.m_zbwd);
351 double klmClusterR_surface = r_BKLM / sin(cluster->getClusterPosition().Theta());
352 if (isEFWD) { klmClusterR_surface = z_EFWD / cos(cluster->getClusterPosition().Theta()); }
353 else if (isEBWD) { klmClusterR_surface = z_EBWD / cos(cluster->getClusterPosition().Theta()); }
355 ROOT::Math::XYZVector klmCluster_surface_position(0, 0, 0);
356 VectorUtil::setMagThetaPhi(klmCluster_surface_position,
357 klmClusterR_surface, cluster->getClusterPosition().Theta(), cluster->getClusterPosition().Phi());
359 double min_dist = 1e10;
361 if (specificTrack !=
nullptr) {
362 const TrackFitResult* trackfit = specificTrack->getTrackFitResultWithBestPValue();
363 ROOT::Math::XYZVector helixPos = getPositionOnHelix_trackfit(trackfit);
365 if (!std::isnan(helixPos.X()) && !std::isnan(helixPos.Y()) && !std::isnan(helixPos.Z())) {
366 min_dist = (klmCluster_surface_position - helixPos).
R();
369 StoreArray<Track> tracks;
370 for (
const Track& track : tracks) {
371 const TrackFitResult* trackfit = track.getTrackFitResultWithBestPValue();
372 ROOT::Math::XYZVector helixPos = getPositionOnHelix_trackfit(trackfit);
374 if (std::isnan(helixPos.X()) || std::isnan(helixPos.Y()) || std::isnan(helixPos.Z()))
continue;
376 double dist = (klmCluster_surface_position - helixPos).
R();
377 if (dist < min_dist) min_dist = dist;
384 double klmClusterTrackDistance_helix_extrapolation(
const Particle* particle)
386 const KLMCluster* cluster = particle->getKLMCluster();
389 double min_dist = ClusterTrackDistance_usingHelixExtrapolate(cluster);
394 double klmClusterTrackSeparationAngle(
const Particle* particle)
396 const KLMCluster* cluster = particle->getKLMCluster();
399 return cluster->getClusterTrackSeparationAngle();
402 double klmClusterTrackRotationAngle(
const Particle* particle)
404 const KLMCluster* cluster = particle->getKLMCluster();
407 return cluster->getClusterTrackRotationAngle();
410 double klmClusterShapeStdDev1(
const Particle* particle)
412 const KLMCluster* cluster = particle->getKLMCluster();
415 return cluster->getShapeStdDev1();
418 double klmClusterShapeStdDev2(
const Particle* particle)
420 const KLMCluster* cluster = particle->getKLMCluster();
423 return cluster->getShapeStdDev2();
426 double klmClusterShapeStdDev3(
const Particle* particle)
428 const KLMCluster* cluster = particle->getKLMCluster();
431 return cluster->getShapeStdDev3();
434 VARIABLE_GROUP(
"KLM Cluster and KlongID");
436 REGISTER_VARIABLE(
"klmClusterKlId", klmClusterKlId,
437 "Returns the KlId classifier output associated to the KLMCluster.");
438 REGISTER_VARIABLE(
"klmClusterBelleTrackFlag", klmClusterBelleTrackFlag,
439 "Returns the Belle-style Track flag.");
440 REGISTER_VARIABLE(
"klmClusterBelleECLFlag", klmClusterBelleECLFlag,
441 "Returns the Belle-style ECL flag.");
442 REGISTER_VARIABLE(
"klmClusterTiming", klmClusterTiming, R
"DOC(
443Returns the timing information of the associated KLMCluster.
446 REGISTER_VARIABLE("klmClusterPositionX", klmClusterPositionX, R"DOC(
447Returns the :math:`x` position of the associated KLMCluster.
450 REGISTER_VARIABLE("klmClusterPositionY", klmClusterPositionY, R"DOC(
451Returns the :math:`y` position of the associated KLMCluster.
454 REGISTER_VARIABLE("klmClusterPositionZ", klmClusterPositionZ, R"DOC(
455Returns the :math:`z` position of the associated KLMCluster.
458 REGISTER_VARIABLE("klmClusterInnermostLayer", klmClusterInnermostLayer,
459 "Returns the number of the innermost KLM layer with a 2-dimensional hit of the associated KLMCluster.");
460 REGISTER_VARIABLE("klmClusterLayers", klmClusterLayers,
461 "Returns the number of KLM layers with 2-dimensional hits of the associated KLMCluster.");
462 REGISTER_VARIABLE("klmClusterEnergy", klmClusterEnergy, R"DOC(
463Returns the energy of the associated KLMCluster. This variable returns an approximation of the energy: it uses `klmClusterMomentum` as momentum and the hypothesis that the KLMCluster is originated by a :math:`K_{L}^0`
465 (: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 `klmClusterMomentum`).
468 MAKE_DEPRECATED("klmClusterEnergy", true, "light-2511-gacrux", R"DOC(
469 As this variable is deemed not physically meaningful it has been deprecated to avoid further use.
471 REGISTER_VARIABLE("klmClusterMomentum", klmClusterMomentum, R
"DOC(
472Returns the momentum magnitude of the associated KLMCluster. This variable returns an approximation of the momentum, since it is proportional to `klmClusterLayers`
474 (:math:`p_{\text{KLM}} = 0.215 \cdot N_{\text{layers}}`, where :math:`p_{\text{KLM}}` is this variable and :math:`N_{\text{layers}}` is `klmClusterLayers`).
477 MAKE_DEPRECATED("klmClusterMomentum", true, "light-2511-gacrux", R"DOC(
478 As this variable is deemed not physically meaningful it has been deprecated to avoid further use.
480 REGISTER_VARIABLE("klmClusterIsBKLM", klmClusterIsBKLM,
481 "Returns 1 if the associated KLMCluster is in barrel KLM.");
482 REGISTER_VARIABLE(
"klmClusterIsEKLM", klmClusterIsEKLM,
483 "Returns 1 if the associated KLMCluster is in endcap KLM.");
484 REGISTER_VARIABLE(
"klmClusterIsForwardEKLM", klmClusterIsForwardEKLM,
485 "Returns 1 if the associated KLMCluster is in forward endcap KLM.");
486 REGISTER_VARIABLE(
"klmClusterIsBackwardEKLM", klmClusterIsBackwardEKLM,
487 "Returns 1 if the associated KLMCluster is in backward endcap KLM.");
488 REGISTER_VARIABLE(
"klmClusterTheta", klmClusterTheta, R
"DOC(
489Returns the polar (:math:`\theta`) angle of the associated KLMCluster.
492 REGISTER_VARIABLE("klmClusterPhi", klmClusterPhi, R"DOC(
493Returns the azimuthal (:math:`\phi`) angle of the associated KLMCluster.
496 REGISTER_VARIABLE("maximumKLMAngleCMS", maximumKLMAngleCMS ,
497 "Returns the maximum angle in the CMS frame between the Particle and all KLMClusters in the event.\n\n","rad");
498 REGISTER_VARIABLE("nKLMClusterTrackMatches", nKLMClusterTrackMatches, R"DOC(
499Returns 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.
501 REGISTER_VARIABLE("nMatchedKLMClusters", nMatchedKLMClusters, R
"DOC(
502 Returns the number of KLMClusters matched to the particle. It only works for
503 Particles created either from Tracks or from ECLCluster, while it returns NaN
504 for :math:`K_{L}^0` or :math:`n` candidates originating from KLMClusters.
506 REGISTER_VARIABLE("nKLMClusterECLClusterMatches", nKLMClusterECLClusterMatches, R
"DOC(
507 Returns the number of ECLClusters matched to the KLMCluster associated to this Particle.
509 REGISTER_VARIABLE("klmClusterTrackDistance_helix_extrapolation", klmClusterTrackDistance_helix_extrapolation, R
"DOC(
510This also returns the distance between KLM cluster and its closes track, but calculated using helix extrapolation. This variable returns ``NAN`` if there is no Track-to-KLMCluster relationship.
513 This is only an approximation of the cluster-track distance that is intended to be used for mDSTs processed with release-09 or earlier. For mDSTs created using release-10 or later please use `klmClusterTrackDistance`.
516 REGISTER_VARIABLE("klmClusterTrackDistance", klmClusterTrackDistance, R"DOC(
517Returns the distance between the KLMCluster associated to this Particle and the closest track to this cluster. This variable returns NaN if there is no Track-to-KLMCluster relationship.
520 This variable only works on mdsts created with release-10 and onwards. For mdsts created with a previous release this variable would
521 return NaN. Please check the release version of the input mdst using
522 `b2file-metadata-show <https://software.belle2.org/development/sphinx/framework/doc/tools/02-b2file.html#b2file-metadata-show-show-the-metadata-of-a-basf2-output-file>`_ with the ``--all`` option.
525 REGISTER_VARIABLE("klmClusterTrackRotationAngle", klmClusterTrackRotationAngle, R"DOC(
526Returns the angle between momenta of the track initially at IP and then at POCA to the KLMCluster associated to this Particle for the closest Track. This variable returns NaN if there is no Track-to-KLMCluster relationship.
529 This variable only works on mdsts created with release-10 and onwards. For mdsts created with a previous release this variable would
530 return NaN. Please check the release version of the input mdst using
531 `b2file-metadata-show <https://software.belle2.org/development/sphinx/framework/doc/tools/02-b2file.html#b2file-metadata-show-show-the-metadata-of-a-basf2-output-file>`_ with the ``--all`` option.
534 REGISTER_VARIABLE("klmClusterTrackSeparationAngle", klmClusterTrackSeparationAngle, R"DOC(
535Returns 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.
538 This variable only works on mdsts created with release-10 and onwards. For mdsts created with a previous release this variable would
539 return NaN. Please check the release version of the input mdst using
540 `b2file-metadata-show <https://software.belle2.org/development/sphinx/framework/doc/tools/02-b2file.html#b2file-metadata-show-show-the-metadata-of-a-basf2-output-file>`_ with the ``--all`` option.
544 REGISTER_VARIABLE("klmClusterShapeStdDev1", klmClusterShapeStdDev1, R"DOC(
545Returns 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 2D Hit.
548 This variable only works on mdsts created with release-10 and onwards. For mdsts created with a previous release this variable would
549 return NaN. Please check the release version of the input mdst using
550 `b2file-metadata-show <https://software.belle2.org/development/sphinx/framework/doc/tools/02-b2file.html#b2file-metadata-show-show-the-metadata-of-a-basf2-output-file>`_ with the ``--all`` option.
553 REGISTER_VARIABLE("klmClusterShapeStdDev2", klmClusterShapeStdDev2, R"DOC(
554Returns the std deviation of the 2nd axis from a PCA of the KLMCluster associated to this Particle. This variable returns 0 if this KLMCluster has all the 2D hits lying on a single straight line.
557 This variable only works on mdsts created with release-10 and onwards. For mdsts created with a previous release this variable would
558 return NaN. Please check the release version of the input mdst using
559 `b2file-metadata-show <https://software.belle2.org/development/sphinx/framework/doc/tools/02-b2file.html#b2file-metadata-show-show-the-metadata-of-a-basf2-output-file>`_ with the ``--all`` option.
562 REGISTER_VARIABLE("klmClusterShapeStdDev3", klmClusterShapeStdDev3, R"DOC(
563Returns the std deviation of the 3nd axis from a PCA of the KLMCluster associated to this Particle. This variable returns 0 if this KLMCluster has all the 2D hits lying on a 2D plane.
566 This variable only works on mdsts created with release-10 and onwards. For mdsts created with a previous release this variable would
567 return NaN. Please check the release version of the input mdst using
568 `b2file-metadata-show <https://software.belle2.org/development/sphinx/framework/doc/tools/02-b2file.html#b2file-metadata-show-show-the-metadata-of-a-basf2-output-file>`_ with the ``--all`` option.
static const double doubleNaN
quiet_NaN
EParticleSourceObject
particle source enumerators
static const std::unordered_map< std::string, DetSurfCylBoundaries > detToSurfBoundaries
Map that associates to each detector its valid cylindrical surface's boundaries.