Belle II Software  release-05-01-25
KlId.h
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Jo-Frederik Krohn *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #pragma once
12 
13 
14 #include <framework/datastore/StoreArray.h>
15 #include <mdst/dataobjects/MCParticle.h>
16 #include <mdst/dataobjects/KLMCluster.h>
17 #include <mdst/dataobjects/ECLCluster.h>
18 #include <mdst/dataobjects/TrackFitResult.h>
19 
20 
21 #include <framework/logging/Logger.h>
22 #include <tracking/dataobjects/RecoTrack.h>
23 #include <genfit/Exception.h>
24 #include <utility>
25 #include <math.h>
26 
28 namespace KlId {
29 
31  int BelleTrackFlag(const Belle2::KLMCluster& cluster, const float angle = 0.26)
32  {
33  const TVector3& pos = cluster.getClusterPosition();
34 
36  for (const Belle2::TrackFitResult& track : tracks) {
37  const TVector3& trackPos = track.getPosition();
38 
39  if (trackPos.Angle(pos) < angle) {
40  B2DEBUG(20, "BelleFlagTracklAngle::" << trackPos.Angle(pos));
41  return 1;
42  }
43  }
44  return 0;
45  }
46 
47 
49  int BelleECLFlag(const Belle2::KLMCluster& cluster, const float angle = 0.26)
50  {
51  const TVector3& pos = cluster.getClusterPosition();
53 
54  for (const Belle2::ECLCluster& eclcluster : eclclusters) {
55 
56  const TVector3& clusterPos = eclcluster.getClusterPosition();
57 
58  if (clusterPos.Angle(pos) < angle) {
59  B2DEBUG(20, "BelleFlagECLAngle::" << clusterPos.Angle(pos));
60  return 1;
61  }
62  }
63  return 0;
64  }
65 
68  {
69  if (part == nullptr) {
70  return 1;
71  } else {
72  return 0;
73  }
74  }
75 
78  {
79  unsigned int hirachy_counter = 0;
80  if (mcParticleIsBeamBKG(part)) {
81  return 0;
82  }
83  bool stop = false;
84  while (!stop) {
85  ++hirachy_counter;
86  if (part -> getPDG() == 130) {
87  return hirachy_counter;
88  }
89  if ((part -> getMother() == nullptr)) {
90  stop = true;
91  } else {
92  part = part -> getMother();
93  }
94  }
95  return 0;
96  }
97 
99  bool isKLMClusterSignal(const Belle2::KLMCluster& cluster, float mcWeigthCut = 0.66)
100  {
101  const auto mcParticleWeightPair = cluster.getRelatedToWithWeight<Belle2::MCParticle>();
102  Belle2::MCParticle* part = mcParticleWeightPair.first;
103  if (!part) {return false; }
104  float mcWeight = mcParticleWeightPair.second;
105  if (mcParticleIsKlong(part) && (mcWeight > mcWeigthCut)) {
106  return true;
107  } else {
108  return false;
109  }
110  }
111 
113  bool isECLClusterSignal(const Belle2::ECLCluster& cluster, float mcWeigthCut = 0.66)
114  {
115  const auto mcParticleWeightPair = cluster.getRelatedToWithWeight<Belle2::MCParticle>();
116  Belle2::MCParticle* part = mcParticleWeightPair.first;
117  if (!part) {return false; }
118  float mcWeight = mcParticleWeightPair.second;
119  if (mcParticleIsKlong(part) && (mcWeight > mcWeigthCut)) {
120  return true;
121  } else {
122  return false;
123  }
124  }
125 
127  int isMCParticlePDG(Belle2::MCParticle* part, int pdg)
128  {
129  bool stop = false;
130  while (!stop) {
131  if (std::abs(part -> getPDG()) == pdg) {
132  return true;
133  }
134  if ((part -> getMother() == nullptr)) {
135  stop = true;
136  } else {
137  part = part -> getMother();
138  }
139  }
140  return false;
141  }
142 
143 
151  {
152 
153  if (mcParticleIsBeamBKG(part)) {
154  return -999;
155  }
156 
157  bool stop = false;
158  while (!stop) {
159  if (isMCParticlePDG(part, 130)) {
160  return 130;
161  }
162  if (isMCParticlePDG(part, 310)) {
163  return 310;
164  }
165  if (isMCParticlePDG(part, 321)) {
166  return 321;
167  }
168  if (isMCParticlePDG(part, 2112)) {
169  return 2112;
170  }
171  if (isMCParticlePDG(part, 2212)) {
172  return 2212;
173  }
174  if (isMCParticlePDG(part, 211)) {
175  return 211;
176  }
177  if (isMCParticlePDG(part, 111)) {
178  return 111;
179  }
180  if (isMCParticlePDG(part, 13)) {
181  return 13;
182  }
183  if (isMCParticlePDG(part, 11)) {
184  return 11;
185  }
186  if (isMCParticlePDG(part, 22)) {
187  return 22;
188  }
189  if ((part -> getMother() == nullptr)) {
190  stop = true;
191  } else {
192  part = part -> getMother();
193  }
194  }
195  return part ->getPDG();
196  }
197 
202  std::pair<Belle2::ECLCluster*, double> findClosestECLCluster(const TVector3& klmClusterPosition,
204  {
205 
206  Belle2::ECLCluster* closestECL = nullptr;
207  double closestECLAngleDist = 1e10;
209 
210  if (eclClusters.getEntries() > 0) {
211  int index = 0;
212  int indexOfClosestCluster = -1;
213  for (Belle2::ECLCluster& eclcluster : eclClusters) {
214 
215  if (eclcluster.hasHypothesis(eclhypothesis)) {
216 
217  const TVector3& eclclusterPos = eclcluster.getClusterPosition();
218  double angularDist = eclclusterPos.Angle(klmClusterPosition);
219  if (angularDist < closestECLAngleDist) {
220  closestECLAngleDist = angularDist;
221  indexOfClosestCluster = index;
222  }
223  }
224  ++index;
225  }
226  if (indexOfClosestCluster > -1)
227  closestECL = eclClusters[indexOfClosestCluster];
228  }
229  return std::make_pair(closestECL, closestECLAngleDist);
230  }
231 
232 
234  std::tuple<const Belle2::KLMCluster*, double, double> findClosestKLMCluster(const TVector3& klmClusterPosition)
235  {
236 
238  const Belle2::KLMCluster* closestKLM = nullptr;
239  double closestKLMDist = 1e10;
240  double avInterClusterDist = 0;
241  double nKLMCluster = klmClusters.getEntries();
242 
243  if (nKLMCluster > 1) {
244 
245  unsigned int index = 0;
246  unsigned int indexOfClosestCluster = 0;
247  for (const Belle2::KLMCluster& nextCluster : klmClusters) {
248 
249  const TVector3& nextClusterPos = nextCluster.getClusterPosition();
250  const TVector3& clustDistanceVec = nextClusterPos - klmClusterPosition;
251 
252  double nextClusterDist = clustDistanceVec.Mag2();
253  avInterClusterDist = avInterClusterDist + nextClusterDist;
254 
255  if ((nextClusterDist < closestKLMDist) and not(nextClusterDist == 0)) {
256  closestKLMDist = nextClusterDist ;
257  indexOfClosestCluster = index;
258  }
259  ++index;
260  }// for next_cluster
261 
262  closestKLM = klmClusters[indexOfClosestCluster];
263  avInterClusterDist = avInterClusterDist / (1. * nKLMCluster);
264  }
265 
266  return std::make_tuple(closestKLM, closestKLMDist, avInterClusterDist);
267  }
268 
269 
271  std::tuple<Belle2::RecoTrack*, double, std::unique_ptr<const TVector3> > findClosestTrack(const TVector3& clusterPosition,
272  float cutAngle)
273  {
275  double oldDistance = 10000000;//dont wanna use infty cos that kills tmva...
276  Belle2::RecoTrack* closestTrack = nullptr;
277  TVector3 poca = TVector3(0, 0, 0);
278 
279 
280  for (Belle2::RecoTrack& track : genfitTracks) {
281 
282  try {
284  genfit::MeasuredStateOnPlane state_for_cut;
285  state_for_cut = track.getMeasuredStateOnPlaneFromFirstHit();
286 
287  // only use tracks that are close cos the extrapolation takes ages
288  if (clusterPosition.Angle(state_for_cut.getPos()) < cutAngle) {
289  state = track.getMeasuredStateOnPlaneFromLastHit();
290  state.extrapolateToPoint(clusterPosition);
291  const TVector3& trackPos = state.getPos();
292 
293  const TVector3& distanceVecCluster = clusterPosition - trackPos;
294  double newDistance = distanceVecCluster.Mag2();
295 
296  // overwrite old distance
297  if (newDistance < oldDistance) {
298  oldDistance = newDistance;
299  closestTrack = &track;
300  poca = trackPos;
301  }
302  }
303  } catch (genfit::Exception& e) {
304  }// try
305  }// for gftrack
306 
307  if (not closestTrack) {
308  return std::make_tuple(closestTrack, oldDistance, std::unique_ptr<const TVector3>(nullptr));
309  } else {
310  return std::make_tuple(closestTrack, oldDistance, std::unique_ptr<const TVector3>(new TVector3(poca)));
311  }
312  }
313 }//end namespace
genfit::Exception
Exception class for error handling in GENFIT (provides storage for diagnostic information)
Definition: Exception.h:48
KlId::BelleTrackFlag
int BelleTrackFlag(const Belle2::KLMCluster &cluster, const float angle=0.26)
get Belle stle track flag
Definition: KlId.h:39
genfit::MeasuredStateOnPlane
#StateOnPlane with additional covariance matrix.
Definition: MeasuredStateOnPlane.h:39
Belle2::ECLCluster
ECL cluster data.
Definition: ECLCluster.h:39
KlId
Helper functions for all klid modules to improve readability of the code.
Definition: KlId.h:28
Belle2::ECLCluster::EHypothesisBit
EHypothesisBit
The hypothesis bits for this ECLCluster (Connected region (CR) is split using this hypothesis.
Definition: ECLCluster.h:43
KlId::mcParticleIsKlong
int mcParticleIsKlong(Belle2::MCParticle *part)
return the mc hirachy of the klong 0:not a klong 1:final particle, 2: klong is mother etc
Definition: KlId.h:85
KlId::BelleECLFlag
int BelleECLFlag(const Belle2::KLMCluster &cluster, const float angle=0.26)
get Belle stle ECL flag
Definition: KlId.h:57
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
KlId::findClosestECLCluster
std::pair< Belle2::ECLCluster *, double > findClosestECLCluster(const TVector3 &klmClusterPosition, const Belle2::ECLCluster::EHypothesisBit eclhypothesis=Belle2::ECLCluster::EHypothesisBit::c_neutralHadron)
Find the closest ECLCluster with a neutral hadron hypothesis, and return it with its distance.
Definition: KlId.h:210
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:78
KlId::isMCParticlePDG
int isMCParticlePDG(Belle2::MCParticle *part, int pdg)
return if mc particle has a certain pdg in the decay chain
Definition: KlId.h:135
Belle2::ECLCluster::EHypothesisBit::c_neutralHadron
@ c_neutralHadron
CR is reconstructed as a neutral hadron (N2)
Belle2::KLMCluster
KLM cluster data.
Definition: KLMCluster.h:38
KlId::findClosestKLMCluster
std::tuple< const Belle2::KLMCluster *, double, double > findClosestKLMCluster(const TVector3 &klmClusterPosition)
find nearest KLMCluster, tis distance and the av intercluster distance
Definition: KlId.h:242
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::StoreArray< Belle2::TrackFitResult >
KlId::isECLClusterSignal
bool isECLClusterSignal(const Belle2::ECLCluster &cluster, float mcWeigthCut=0.66)
checks if a cluster is signal under the mcWeightcondition (mcWeight = energy deposition)
Definition: KlId.h:121
KlId::mcParticleIsBeamBKG
int mcParticleIsBeamBKG(Belle2::MCParticle *part)
return if MCparticle is beambkg
Definition: KlId.h:75
Belle2::DistanceTools::poca
TVector3 poca(TVector3 const &trackPos, TVector3 const &trackP, TVector3 const &vtxPos)
Returns the Point Of Closest Approach of a track to a vertex.
Definition: DistanceTools.cc:21
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
KlId::findClosestTrack
std::tuple< Belle2::RecoTrack *, double, std::unique_ptr< const TVector3 > > findClosestTrack(const TVector3 &clusterPosition, float cutAngle)
find nearest genfit track and return it and its distance
Definition: KlId.h:279
KlId::getPrimaryPDG
int getPrimaryPDG(Belle2::MCParticle *part)
return if mc particles primary pdg.
Definition: KlId.h:158
KlId::isKLMClusterSignal
bool isKLMClusterSignal(const Belle2::KLMCluster &cluster, float mcWeigthCut=0.66)
checks if a cluster is signal under the mcWeightcondition (mcWeight = energy deposition)
Definition: KlId.h:107