Belle II Software  release-05-02-19
ECLTrackBremFinderModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * Contributors: Thomas Hauth, Patrick Ecker *
5  * *
6  * This software is provided "as is" without any warranty. *
7  **************************************************************************/
8 //This module
9 #include <ecl/modules/eclTrackBremFinder/ECLTrackBremFinderModule.h>
10 
11 //Framework
12 #include <framework/dataobjects/EventMetaData.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/gearbox/Const.h>
15 #include <framework/datastore/StoreArray.h>
16 
17 //MDST
18 #include <mdst/dataobjects/ECLCluster.h>
19 #include <mdst/dataobjects/Track.h>
20 #include <mdst/dataobjects/PIDLikelihood.h>
21 #include <mdst/dataobjects/TrackFitResult.h>
22 
23 //Tracking
24 #include <tracking/dataobjects/RecoTrack.h>
25 #include <tracking/dataobjects/BremHit.h>
26 
27 //ECL
28 #include <ecl/modules/eclTrackBremFinder/BestMatchContainer.h>
29 #include <ecl/modules/eclTrackBremFinder/BremFindingMatchCompute.h>
30 
31 
32 using namespace Belle2;
33 
34 REG_MODULE(ECLTrackBremFinder)
35 
37  Module()
38 {
39  setDescription("Use Track direction to pick up possible ECL Brem Cluster");
40  setPropertyFlags(c_ParallelProcessingCertified);
41 
42  addParam("eclClustersStoreArrayName", m_param_eclClustersStoreArrayName, "StoreArray name of the ECLClusters for brem matching",
43  m_param_eclClustersStoreArrayName);
44 
45  addParam("tracksStoreArrayName", m_param_tracksStoreArrayName, "StoreArray name of the Tracks for brem matching",
46  m_param_tracksStoreArrayName);
47 
48  addParam("recoTracksStoreArrayName", m_param_recoTracksStoreArrayName, "StoreArray name of the reco tracks used for brem search",
49  m_param_recoTracksStoreArrayName);
50 
51  addParam("clusterAcceptanceFactor", m_clusterAcceptanceFactor,
52  "Factor which is multiplied onto the cluster position error to check for matches",
53  m_clusterAcceptanceFactor);
54 
55  addParam("virtualHitRadii", m_virtualHitRadii, "Radii where virtual hits for the extrapolation will be generated",
56  m_virtualHitRadii);
57 
58  addParam("relativeClusterEnergy", m_relativeClusterEnergy, "Fraction of the tracks energy the ECL cluster has "
59  "to possess to be considered for bremsstrahlung finding",
60  m_relativeClusterEnergy);
61 
62  addParam("requestedNumberOfCDCHits", m_requestedNumberOfCDCHits, "Minimal/Maximal number of CDC hits, the track has to possess "
63  "to be considered for bremsstrahlung finding",
64  m_requestedNumberOfCDCHits);
65 
66  addParam("electronProbabilityCut", m_electronProbabilityCut, "Cut on the electron probability (from pid) of track",
67  m_electronProbabilityCut);
68 
69  addParam("clusterDistanceCut", m_clusterDistanceCut,
70  "Cut on the distance between the cluster position angle and the extrapolation angle",
71  m_clusterDistanceCut);
72 }
73 
75 {
77  m_eclClusters.registerRelationTo(m_eclClusters);
78 
80 
81  const std::string relationName = "Bremsstrahlung";
82  m_eclClusters.registerRelationTo(m_tracks, DataStore::c_Event, DataStore::c_WriteOut, relationName);
83 
84  m_recoTracks.isRequired();
85 
86  m_bremHits.registerInDataStore();
87  m_bremHits.registerRelationTo(m_eclClusters);
88  m_bremHits.registerRelationTo(m_recoTracks);
89 }
90 
92 {
93 
94 
95  // either use the Clusters matched to tracks (non-neutral) or use the smarter decision
96  // done by the neutral / non-neutral classification code
97  // todo: there needs to be a global (all tracks vs. all clusters) conflict resolution,
98  // sort tracks by Pt, as high energy tracks are more likely to radiate bremsstrahlung photons
99  // with sufficient energy to be detected and reconstructed
100  for (auto& track : m_tracks) {
101 
102  // since the module runs after the reconstruction the pid likelihood can be checked to sort out tracks,
103  // which are not expected to be from electrons
104  const PIDLikelihood* pid = track.getRelated<PIDLikelihood>();
105  if (pid) {
107  Const::ChargedStable mostLikelyPDG = Const::electron;
108  double highestProb = 0;
109  for (Const::ChargedStable pdg : possiblePDGs) {
110  double probability = pid->getProbability(pdg);
111  if (probability > highestProb) {
112  highestProb = probability;
113  mostLikelyPDG = pdg;
114  }
115  }
116  if (mostLikelyPDG != Const::electron || highestProb <= m_electronProbabilityCut) {
117  B2DEBUG(20, "Track is expected not to be from electron");
118  continue;
119  }
120  }
121 
122  const TrackFitResult* trackFitResult = track.getTrackFitResult(Const::ChargedStable(211));
123  double trackMomentum;
124  if (trackFitResult) {
125  trackMomentum = trackFitResult->getMomentum().Mag();
126  // if the momentum of the track is higher than 5 GeV, do not use this track
127  if (trackMomentum > 5.0) {
128  B2DEBUG(20, "Track momentum higher than 5GeV! Track is not used for bremsstrahlung finding");
129  continue;
130  }
131  } else {
132  continue;
133  }
134 
135 
136  B2DEBUG(20, "Checking track for related ECLCluster");
137 
138  // searching for an assigned primary cluster
139  ECLCluster* primaryClusterOfTrack = nullptr;
140  auto relatedClustersToTrack =
141  track.getRelationsWith<ECLCluster>
142  (m_param_eclClustersStoreArrayName); //check the cluster hypothesis ID here (take c_nPhotons hypothesis)!!
143  for (auto& relatedCluster : relatedClustersToTrack) {
144  if (relatedCluster.hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
145  primaryClusterOfTrack = &relatedCluster;
146  }
147  }
148 
149  // get the RecoTrack to have easy access to individual hits and
150  // their fit state
151  auto recoTrack = track.getRelatedTo<RecoTrack>(m_param_recoTracksStoreArrayName);
152 
153  if (!recoTrack) {
154  // no reco track
155  B2DEBUG(20, "No RecoTrack for this Track");
156  continue;
157  }
158 
159  // check if the RecoTrack has the requested number of CDC hits
160  if (recoTrack->getNumberOfCDCHits() < std::get<0>(m_requestedNumberOfCDCHits) ||
161  recoTrack->getNumberOfCDCHits() > std::get<1>(m_requestedNumberOfCDCHits)) {
162  B2DEBUG(20, "RecoTrack has not requested number of CDC hits");
163  continue;
164  }
165 
166  // set the params for the virtual hits
167  std::vector<std::pair<float, RecoHitInformation*>> extrapolationParams = {};
168  std::vector<genfit::MeasuredStateOnPlane> extrapolatedStates = {};
169  try {
170  for (auto virtualHitRadius : m_virtualHitRadii) {
172  for (auto hit : recoTrack->getRecoHitInformations(true)) {
173  if (hit->useInFit() && recoTrack->hasTrackFitStatus()) {
174  try {
175  auto measState = recoTrack->getMeasuredStateOnPlaneFromRecoHit(hit);
176  float hitRadius = measState.getPos().Perp();
177  float distance = abs(hitRadius - virtualHitRadius);
178  // for higher values the extrapolation will be too bad
179  if (distance < 3) {
180  nearestHitContainer.add(hit, distance);
181  }
182  } catch (NoTrackFitResult&) {
183  B2DEBUG(29, "No track fit result available for this hit! Event: " << m_evtPtr->getEvent());
184  }
185  }
186  if (nearestHitContainer.hasMatch()) {
187  auto nearestHit = nearestHitContainer.getBestMatch();
188  extrapolationParams.push_back({virtualHitRadius, nearestHit});
189  }
190  }
191  }
192 
193  for (auto param : extrapolationParams) {
194  auto fitted_state = recoTrack->getMeasuredStateOnPlaneFromRecoHit(param.second);
195  try {
196  fitted_state.extrapolateToCylinder(param.first);
197  extrapolatedStates.push_back(fitted_state);
198  } catch (genfit::Exception& exception1) {
199  B2DEBUG(20, "Extrapolation failed!");
200  }
201  }
202  } catch (const genfit::Exception& e) {
203  B2WARNING("Exception" << e.what());
204  }
205 
206  // possible improvement: use fast lookup using kd-tree, this is nasty
207  // iterate over full cluster list to find possible compatible clusters
208  for (ECLCluster& cluster : m_eclClusters) {
209  //check if the cluster belongs to a photon or electron
210  if (!cluster.hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
211  B2DEBUG(20, "Cluster has wrong hypothesis!");
212  continue;
213  }
214 
215  //check if cluster is already related to a track -> if true, can't be bremsstrahlung cluster
216  auto relatedTrack = cluster.getRelatedFrom<Track>();
217  if (relatedTrack) {
218  B2DEBUG(20, "Cluster already related to track, bailing out");
219  continue;
220  }
221 
222  //check if the cluster is already related to a BremHit
223  //procedure: first come, first served
224  auto relatedBremHit = cluster.getRelated<BremHit>();
225  if (relatedBremHit) {
226  B2DEBUG(20, "Cluster already assumed to be bremsstrahlung cluster!");
227  continue;
228  }
229 
230  // check if the cluster energy is higher than the track momentum itself
231  // also check that cluster has more than 2% of the track momentum
232  double relativeEnergy = cluster.getEnergy(ECLCluster::EHypothesisBit::c_nPhotons) / trackMomentum;
233  if (relativeEnergy > 1 || relativeEnergy < m_relativeClusterEnergy) {
234  B2DEBUG(20, "Relative energy of cluster higher than 1 or below threshold!");
235  continue;
236  }
237 
238  typedef std::tuple<ECLCluster*, genfit::MeasuredStateOnPlane, double, double > ClusterMSoPPair;
240 
241  // iterate over all track points and see whether this cluster matches
242  // the track points direction
243  // only VXD hits shall be used
244  for (auto hit : recoTrack->getRecoHitInformations(true)) {
245  if (hit->getTrackingDetector() == RecoHitInformation::c_PXD || hit->getTrackingDetector() == RecoHitInformation::c_SVD) {
246  try {
247  if (!recoTrack->hasTrackFitStatus()) {
248  continue;
249  }
250  auto measState = recoTrack->getMeasuredStateOnPlaneFromRecoHit(hit);
251  auto bremFinder = BremFindingMatchCompute(m_clusterAcceptanceFactor, cluster, measState);
252  if (bremFinder.isMatch()) {
253  ClusterMSoPPair match_pair = std::make_tuple(&cluster, measState, bremFinder.getDistanceHitCluster(),
254  bremFinder.getEffAcceptanceFactor());
255  matchContainer.add(match_pair, bremFinder.getDistanceHitCluster());
256  }
257  } catch (NoTrackFitResult&) {
258  B2DEBUG(29, "No track fit result available for this hit! Event: " << m_evtPtr->getEvent());
259  } catch (genfit::Exception& e) {
260  B2WARNING("Exception" << e.what());
261  }
262  }
263  }
264 
265  // check for matches of the extrapolation of the virtual hits with the cluster position
266  for (auto fitted_state : extrapolatedStates) {
267  auto bremFinder = BremFindingMatchCompute(m_clusterAcceptanceFactor, cluster, fitted_state);
268  if (bremFinder.isMatch()) {
269  ClusterMSoPPair match_pair = std::make_tuple(&cluster, fitted_state, bremFinder.getDistanceHitCluster(),
270  bremFinder.getEffAcceptanceFactor());
271  matchContainer.add(match_pair, bremFinder.getDistanceHitCluster());
272  }
273  }
274 
275 
276  // have we found the best possible track point for this cluster
277  if (matchContainer.hasMatch()) {
278  auto matchClustermSoP = matchContainer.getBestMatch();
279  const auto fitted_state = std::get<1>(matchClustermSoP);
280 
281  const auto fitted_pos = fitted_state.getPos();
282  const auto fitted_mom = fitted_state.getMom();
283  const auto fitted_dir = fitted_state.getDir();
284 
285  const auto hit_theta = fitted_mom.Theta();
286  const auto hit_phi = fitted_mom.Phi();
287 
288  B2DEBUG(20, "Best Cluster" << std::endl
289  << " Cluster Phi=" << std::get<0>(matchClustermSoP)->getPhi() << " Theta=" << std::get<0>(matchClustermSoP)->getTheta()
290  << " TrackHit Phi=" << hit_phi << " Theta=" << hit_theta);
291 
292  // create a BremHit if a match is found
293  // relate this BremHit to the bremsstrahlung cluster and the recoTrack
294  // if the track has a primary cluster, add a relation between bremsstrahlung cluster and primary cluster
295 
296  double effAcceptanceFactor = std::get<3>(matchClustermSoP);
297  ECLCluster* bremCluster = std::get<0>(matchClustermSoP);
298  double clusterDistance = std::get<2>(matchClustermSoP);
299 
300  if (fitted_pos.Perp() <= 16 && clusterDistance <= m_clusterDistanceCut) {
301  m_bremHits.appendNew(recoTrack, bremCluster,
302  fitted_pos, bremCluster->getEnergy(ECLCluster::EHypothesisBit::c_nPhotons),
303  clusterDistance, effAcceptanceFactor);
304 
305  // add a relation between the bremsstrahlung cluster and the track to transfer the information to the analysis
306  // set the acceptance factor as weight
307  bremCluster->addRelationTo(&track, effAcceptanceFactor, "Bremsstrahlung");
308 
309  if (primaryClusterOfTrack) {
310  primaryClusterOfTrack->addRelationTo(bremCluster, effAcceptanceFactor);
311  }
312  }
313  }
314  }
315  }
316 }
genfit::Exception
Exception class for error handling in GENFIT (provides storage for diagnostic information)
Definition: Exception.h:48
Belle2::ECLTrackBremFinderModule::m_evtPtr
StoreObjPtr< EventMetaData > m_evtPtr
StoreObjPtr EventMetaData.
Definition: ECLTrackBremFinderModule.h:75
Belle2::TrackFitResult::getMomentum
TVector3 getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:116
Belle2::ECLCluster
ECL cluster data.
Definition: ECLCluster.h:39
Belle2::Const::electron
static const ChargedStable electron
electron particle
Definition: Const.h:533
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::ECLTrackBremFinderModule::event
virtual void event() override
Called once for each event.
Definition: ECLTrackBremFinderModule.cc:91
Belle2::ECLTrackBremFinderModule::m_param_tracksStoreArrayName
std::string m_param_tracksStoreArrayName
StoreArray name of the Tracks for brem matching.
Definition: ECLTrackBremFinderModule.h:116
Belle2::BestMatchContainer::getBestMatch
TItem const & getBestMatch() const
Definition: BestMatchContainer.h:82
Belle2::BestMatchContainer::hasMatch
bool hasMatch() const
Definition: BestMatchContainer.h:73
Belle2::ECLCluster::EHypothesisBit::c_nPhotons
@ c_nPhotons
CR is split into n photons (N1)
Belle2::RelationsInterface::addRelationTo
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
Definition: RelationsObject.h:144
Belle2::PIDLikelihood
Class to collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Definition: PIDLikelihood.h:37
Belle2::ECLTrackBremFinderModule::m_virtualHitRadii
std::vector< float > m_virtualHitRadii
Radii where virtual hits for the extrapolation will be generated The default values are taken from br...
Definition: ECLTrackBremFinderModule.h:86
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
Belle2::DataStore::c_WriteOut
@ c_WriteOut
Object/array should be saved by output modules.
Definition: DataStore.h:72
Belle2::Const::kaon
static const ChargedStable kaon
charged kaon particle
Definition: Const.h:536
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::ECLTrackBremFinderModule::m_recoTracks
StoreArray< RecoTrack > m_recoTracks
StoreArray RecoTrack.
Definition: ECLTrackBremFinderModule.h:71
Belle2::ECLCluster::getEnergy
double getEnergy(const EHypothesisBit &hypothesis) const
Return Energy (GeV).
Definition: ECLCluster.cc:21
Belle2::ECLTrackBremFinderModule::m_eclClusters
StoreArray< ECLCluster > m_eclClusters
StoreArray ECLCluster.
Definition: ECLTrackBremFinderModule.h:67
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:78
Belle2::Const::pion
static const ChargedStable pion
charged pion particle
Definition: Const.h:535
Belle2::ECLTrackBremFinderModule::m_relativeClusterEnergy
double m_relativeClusterEnergy
Fraction of the tracks energy the ECL cluster has to possess to be considered for bremsstrahlung find...
Definition: ECLTrackBremFinderModule.h:91
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECLTrackBremFinderModule::m_bremHits
StoreArray< BremHit > m_bremHits
StoreArray BremHits.
Definition: ECLTrackBremFinderModule.h:73
Belle2::Const::deuteron
static const ChargedStable deuteron
deuteron particle
Definition: Const.h:538
Belle2::ECLTrackBremFinderModule::m_electronProbabilityCut
float m_electronProbabilityCut
Cut on the electron probability (from pid) of track.
Definition: ECLTrackBremFinderModule.h:101
Belle2::ECLTrackBremFinderModule
Module to assign ECL Clusters resulting from Bremsstrahlung to the primary electron track.
Definition: ECLTrackBremFinderModule.h:44
Belle2::ECLTrackBremFinderModule::m_clusterDistanceCut
float m_clusterDistanceCut
Cut on the distance between the cluster position angle and the extrapolation angle.
Definition: ECLTrackBremFinderModule.h:106
Belle2::ECLTrackBremFinderModule::m_param_eclClustersStoreArrayName
std::string m_param_eclClustersStoreArrayName
StoreArray name of the ECLClusters for brem matching.
Definition: ECLTrackBremFinderModule.h:111
Belle2::ECLTrackBremFinderModule::m_tracks
StoreArray< Track > m_tracks
StoreArray Track.
Definition: ECLTrackBremFinderModule.h:69
Belle2::BremFindingMatchCompute
Module to compute if an extrapolation to the ECL matches the position of an secondary ECLCLuster to f...
Definition: BremFindingMatchCompute.h:28
Belle2::ECLTrackBremFinderModule::m_param_recoTracksStoreArrayName
std::string m_param_recoTracksStoreArrayName
StoreArray name of the RecoTracks for brem matching.
Definition: ECLTrackBremFinderModule.h:121
Belle2::BremHit
A bremsstrahlung hit that correlates an ECLCluster with a RecoTrack.
Definition: BremHit.h:30
Belle2::Const::proton
static const ChargedStable proton
proton particle
Definition: Const.h:537
Belle2::Const::ChargedStable
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:465
Belle2::ECLTrackBremFinderModule::m_clusterAcceptanceFactor
double m_clusterAcceptanceFactor
Factor which is multiplied onto the cluster position error to check for matches.
Definition: ECLTrackBremFinderModule.h:79
Belle2::ECLTrackBremFinderModule::initialize
virtual void initialize() override
Use this to initialize resources or memory your module needs.
Definition: ECLTrackBremFinderModule.cc:74
Belle2::Const::muon
static const ChargedStable muon
muon particle
Definition: Const.h:534
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
Belle2::DataStore::c_Event
@ c_Event
Different object in each event, all objects/arrays are invalidated after event() function has been ca...
Definition: DataStore.h:61
Belle2::BestMatchContainer::add
bool add(TItem const &item, TEstimator est, EstimatorComparison estComparison=[](TEstimator currentBest, TEstimator newEst) {return newEst< currentBest;})
Add a new item with an estimator value.
Definition: BestMatchContainer.h:50
Belle2::BestMatchContainer
Multiple entries can be added, but only the one will be kept, which has the best quality estimator.
Definition: BestMatchContainer.h:35
Belle2::ECLTrackBremFinderModule::m_requestedNumberOfCDCHits
std::tuple< unsigned int, unsigned int > m_requestedNumberOfCDCHits
Minimal/Maximal number of CDC hits, the track has to possess to be considered for bremsstrahlung find...
Definition: ECLTrackBremFinderModule.h:96