Belle II Software  release-05-01-25
StereoHitTrackQuadTreeMatcher.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Nils Braun, Dmitrii Neverov *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 #include <tracking/trackFindingCDC/collectors/matchers/StereoHitTrackQuadTreeMatcher.h>
11 
12 #include <tracking/trackFindingCDC/hough/z0_tanLambda/HitZ0TanLambdaLegendre.h>
13 #include <tracking/trackFindingCDC/hough/quadratic/HitQuadraticLegendre.h>
14 #include <tracking/trackFindingCDC/hough/hyperbolic/HitHyperHough.h>
15 
16 #include <tracking/trackFindingCDC/eventdata/tracks/CDCTrack.h>
17 
18 #include <tracking/trackFindingCDC/eventdata/hits/CDCRLWireHit.h>
19 #include <tracking/trackFindingCDC/eventdata/hits/CDCWireHit.h>
20 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCTrajectory2D.h>
21 
22 #include <tracking/trackFindingCDC/ca/AutomatonCell.h>
23 
24 #include <framework/core/ModuleParamList.templateDetails.h>
25 
26 #include <tracking/trackFindingCDC/utilities/StringManipulation.h>
27 
28 #include <utility>
29 
30 using namespace Belle2;
31 using namespace TrackFindingCDC;
32 
33 template <class AQuadTree>
35  const std::string& prefix)
36 {
37  Super::exposeParameters(moduleParamList, prefix);
38 
39  moduleParamList->addParameter(prefixed(prefix, "level"), m_param_quadTreeLevel,
40  "The number of levels for the quad tree search.",
41  m_param_quadTreeLevel);
42 
43  moduleParamList->addParameter(prefixed(prefix, "minimumNumberOfHits"), m_param_minimumNumberOfHits,
44  "The minimum number of hits in a quad tree bin to be called as result.",
45  m_param_minimumNumberOfHits);
46 
47  moduleParamList->addParameter(prefixed(prefix, "writeDebugInformation"), m_param_writeDebugInformation,
48  "Set to true to output debug information.",
49  m_param_writeDebugInformation);
50 
51  moduleParamList->addParameter(prefixed(prefix, "checkForB2BTracks"),
52  m_param_checkForB2BTracks,
53  "Set to false to skip the check for back-2-back tracks "
54  "(good for cosmics).",
55  m_param_checkForB2BTracks);
56 
57  moduleParamList->addParameter(prefixed(prefix, "checkForInWireBoundsFactor"),
58  m_param_checkForInWireBoundsFactor,
59  "Used to scale the CDC before checking "
60  "whether hits are in the CDC z bounds.",
61  m_param_checkForInWireBoundsFactor);
62 }
63 
64 
65 template <class AQuadTree>
67 {
68  Super::initialize();
69 
70  m_quadTreeInstance.setMaxLevel(m_param_quadTreeLevel);
71  m_quadTreeInstance.initialize();
72 }
73 
74 
75 template <class AQuadTree>
77 {
78  Super::terminate();
79 
80  m_quadTreeInstance.raze();
81 }
82 
83 template <class AQuadTree>
84 void StereoHitTrackQuadTreeMatcher<AQuadTree>::match(CDCTrack& track, const std::vector<CDCRLWireHit>& rlWireHits,
85  std::vector<Super::WeightedRelationItem>& relationsForCollector)
86 {
87  // TODO: Extract this into smaller steps
88  // TODO: Split the filtering from the rest.
89  // This means this method would output WeightedRelations based on the quad tree decision and a second filter step
90  // can be applied on these weighted relations
91 
92  // Reconstruct the hits to the track
93  const CDCTrajectory2D& trajectory2D = track.getStartTrajectory3D().getTrajectory2D();
94  const bool isCurler = trajectory2D.isCurler();
95 
96  using CDCRecoHitWithRLPointer = std::pair<CDCRecoHit3D, const CDCRLWireHit*>;
97  std::vector<CDCRecoHitWithRLPointer> recoHits;
98  recoHits.reserve(rlWireHits.size() + track.size());
99 
100  /*
101  * Use the given trajectory to reconstruct the 2d hits in the vector in z direction
102  * to match the trajectory perfectly. Then add the newly created reconstructed 3D hit to the given list.
103  */
104  for (const CDCRLWireHit& rlWireHit : rlWireHits) {
105  if (rlWireHit.getWireHit().getAutomatonCell().hasTakenFlag()) continue;
106 
107  const CDCWire& wire = rlWireHit.getWire();
108  const WireLine& wireLine = wire.getWireLine();
109  double signedDriftLength = rlWireHit.getSignedRefDriftLength();
110  for (const Vector3D& recoPos3D : trajectory2D.reconstructBoth3D(wireLine, signedDriftLength)) {
111  // Skip hits out of CDC
112  if (not wire.isInCellZBounds(recoPos3D, m_param_checkForInWireBoundsFactor)) {
113  continue;
114  }
115 
116  // If the track is a curler, shift all perpS values to positive ones.
117  // Else do not use this hit if m_param_checkForB2BTracks is enabled.
118  double perpS = trajectory2D.calcArcLength2D(recoPos3D.xy());
119  if (perpS < 0) {
120  if (isCurler) {
121  perpS += trajectory2D.getArcLength2DPeriod();
122  } else if (m_param_checkForB2BTracks) {
123  continue;
124  }
125  }
126  recoHits.emplace_back(CDCRecoHit3D(rlWireHit, recoPos3D, perpS), &rlWireHit);
127  }
128  }
129 
130  // Also add already found stereo hits of the track
131  for (const CDCRecoHit3D& recoHit : track) {
132  if (not recoHit.isAxial()) {
133  recoHit.getWireHit().getAutomatonCell().setAssignedFlag();
134  const CDCRLWireHit& rlWireHit = recoHit.getRLWireHit();
135  recoHits.emplace_back(recoHit, &rlWireHit);
136  }
137  }
138 
139  // Do the tree finding
140  m_quadTreeInstance.seed(recoHits);
141 
142  if (m_param_writeDebugInformation) {
143  writeDebugInformation();
144  }
145 
146  const auto& foundStereoHitsWithNode = m_quadTreeInstance.findSingleBest(m_param_minimumNumberOfHits);
147  m_quadTreeInstance.fell();
148 
149  if (foundStereoHitsWithNode.size() != 1) {
150  return;
151  }
152 
153  // There is the possibility that we have added one cdc hits twice (as left and right one). We search for those cases here:
154  auto foundStereoHits = foundStereoHitsWithNode[0].second;
155  const auto& node = foundStereoHitsWithNode[0].first;
156 
157  if (m_param_writeDebugInformation) {
158  std::vector<CDCRecoHit3D> allHits;
159  std::vector<CDCRecoHit3D> foundHits;
160  // Turn vector of pairs into vector of first items
161  for (const CDCRecoHitWithRLPointer recoHitWithRL : recoHits) {
162  const CDCRecoHit3D& recoHit3D = recoHitWithRL.first;
163  allHits.push_back(recoHit3D);
164  }
165  for (const CDCRecoHitWithRLPointer recoHitWithRL : foundStereoHits) {
166  const CDCRecoHit3D& recoHit3D = recoHitWithRL.first;
167  foundHits.push_back(recoHit3D);
168  }
169  m_quadTreeInstance.drawDebugPlot(allHits, foundHits, node);
170  }
171 
172  // Remove all assigned hits, which where already found before (and do not need to be added again)
173  const auto& isAssignedHit = [](const CDCRecoHitWithRLPointer & recoHitWithRLPointer) {
174  const CDCRecoHit3D& recoHit3D = recoHitWithRLPointer.first;
175  const auto& automatonCell = recoHit3D.getWireHit().getAutomatonCell();
176  return automatonCell.hasAssignedFlag();
177  };
178 
179  foundStereoHits.erase(std::remove_if(foundStereoHits.begin(),
180  foundStereoHits.end(),
181  isAssignedHit),
182  foundStereoHits.end());
183 
184  // Sort the found stereo hits by same CDCHit and smaller distance to the node
185  auto sortByHitAndNodeCenterDistance = [node](const CDCRecoHitWithRLPointer & lhs,
186  const CDCRecoHitWithRLPointer & rhs) {
187 
188 
189  const CDCRecoHit3D& rhsRecoHit = rhs.first;
190  const CDCRecoHit3D& lhsRecoHit = lhs.first;
191 
192  const CDCWireHit& rhsWireHit = rhsRecoHit.getWireHit();
193  const CDCWireHit& lhsWireHit = lhsRecoHit.getWireHit();
194 
195  if (lhsWireHit < rhsWireHit) {
196  return true;
197  } else if (rhsWireHit < lhsWireHit) {
198  return false;
199  } else {
200  return AQuadTree::DecisionAlgorithm::BoxAlgorithm::compareDistances(node, lhsRecoHit, rhsRecoHit); //returns true if lhs < rhs
201  }
202  };
203 
204  const auto& sameHitComparer = [](const CDCRecoHitWithRLPointer & lhs,
205  const CDCRecoHitWithRLPointer & rhs) {
206  const CDCRecoHit3D& rhsRecoHit = rhs.first;
207  const CDCRecoHit3D& lhsRecoHit = lhs.first;
208 
209  return lhsRecoHit.getWireHit() == rhsRecoHit.getWireHit();
210  };
211 
212  std::sort(foundStereoHits.begin(),
213  foundStereoHits.end(),
214  sortByHitAndNodeCenterDistance);
215 
216  // If the same hit is added as right and left hypothesis, do only use the one with the smaller distance to the node.
217  foundStereoHits.erase(std::unique(foundStereoHits.begin(),
218  foundStereoHits.end(),
219  sameHitComparer),
220  foundStereoHits.end());
221 
222  // Add the found stereo hits to the relation vector. In the moment, all hits get the same weight (may change later).
223  for (const CDCRecoHitWithRLPointer& recoHitWithRLPointer : foundStereoHits) {
224  const CDCRLWireHit* rlWireHit = recoHitWithRLPointer.second;
225  relationsForCollector.emplace_back(&track, foundStereoHits.size(), rlWireHit);
226  }
227 }
228 
229 template <class AQuadTree>
231 {
232  std::string outputFileName = "quadTreeContent_call_" + std::to_string(m_numberOfPassedDebugCalls) + ".root";
233  m_quadTreeInstance.writeDebugInfoToFile(outputFileName);
234 
235  m_numberOfPassedDebugCalls++;
236 }
237 
Belle2::TrackFindingCDC::CDCWireHit::getAutomatonCell
AutomatonCell & getAutomatonCell() const
Mutable getter for the automaton cell.
Definition: CDCWireHit.h:294
Belle2::TrackFindingCDC::CDCRecoHit3D
Class representing a three dimensional reconstructed hit.
Definition: CDCRecoHit3D.h:62
Belle2::TrackFindingCDC::AutomatonCell::hasAssignedFlag
bool hasAssignedFlag() const
Gets the current state of the already assigned marker flag.
Definition: AutomatonCell.h:150
Belle2::TrackFindingCDC::CDCTrack
Class representing a sequence of three dimensional reconstructed hits.
Definition: CDCTrack.h:51
Belle2::TrackFindingCDC::StereoHitTrackQuadTreeMatcher::writeDebugInformation
void writeDebugInformation()
Use the writeDebugInformation function of the quad tree to write the tree into a root file with a asc...
Definition: StereoHitTrackQuadTreeMatcher.cc:230
Belle2::TrackFindingCDC::StereoHitTrackQuadTreeMatcher
A matcher algorithm for using a stereo quad tree for matching rl tagged wire hits to tracks.
Definition: StereoHitTrackQuadTreeMatcher.h:35
Belle2::TrackFindingCDC::CDCWire::isInCellZBounds
bool isInCellZBounds(const Vector3D &pos3D, const double factor=1) const
Checks whether the position is in the z bounds of the drift cell (scaled by the factor) surrounding t...
Definition: CDCWire.h:297
Belle2::TrackFindingCDC::CDCRecoHit3D::getWireHit
const CDCWireHit & getWireHit() const
Getter for the wire hit.
Definition: CDCRecoHit3D.h:248
Belle2::TrackFindingCDC::CDCTrajectory2D::isCurler
bool isCurler(double factor=1) const
Checks if the trajectory leaves the outer radius of the CDC times the given tolerance factor.
Definition: CDCTrajectory2D.cc:267
Belle2::TrackFindingCDC::CDCTrajectory2D::reconstructBoth3D
std::array< Vector3D, 2 > reconstructBoth3D(const WireLine &wireLine, double distance=0.0, double z=0) const
Gives the two three dimensional points where the drift circle touches the wire line.
Definition: CDCTrajectory2D.cc:151
Belle2::TrackFindingCDC::CDCTrajectory2D::getArcLength2DPeriod
double getArcLength2DPeriod() const
Getter for the arc length for one round trip around the trajectory.
Definition: CDCTrajectory2D.h:289
Belle2::TrackFindingCDC::CDCTrajectory2D
Particle trajectory as it is seen in xy projection represented as a circle.
Definition: CDCTrajectory2D.h:46
Belle2::ModuleParamList::addParameter
void addParameter(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module list.
Definition: ModuleParamList.templateDetails.h:38
Belle2::TrackFindingCDC::StereoHitTrackQuadTreeMatcher::exposeParameters
void exposeParameters(ModuleParamList *moduleParamList, const std::string &prefix) override
Expose the parameters to the module.
Definition: StereoHitTrackQuadTreeMatcher.cc:34
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TrackFindingCDC::CDCRLWireHit
Class representing an oriented hit wire including a hypotheses whether the causing track passes left ...
Definition: CDCRLWireHit.h:51
Belle2::TrackFindingCDC::WireLine
A three dimensional limited line represented by its closest approach to the z-axes (reference positio...
Definition: WireLine.h:41
Belle2::TrackFindingCDC::Vector3D
A three dimensional vector.
Definition: Vector3D.h:34
Belle2::TrackFindingCDC::StereoHitTrackQuadTreeMatcher::terminate
void terminate() override
Terminate the filter and the quad tree.
Definition: StereoHitTrackQuadTreeMatcher.cc:76
Belle2::TrackFindingCDC::StereoHitTrackQuadTreeMatcher::initialize
void initialize() override
Initialize the filter and the quad tree.
Definition: StereoHitTrackQuadTreeMatcher.cc:66
Belle2::TrackFindingCDC::StereoHitTrackQuadTreeMatcher::match
void match(CDCTrack &track, const std::vector< CDCRLWireHit > &rlWireHits, std::vector< Super::WeightedRelationItem > &relationsForCollector) override
Create a QuadTree and fill with each unused stereo hit (to be exact: twice for each stereo hit - righ...
Definition: StereoHitTrackQuadTreeMatcher.cc:84
Belle2::TrackFindingCDC::CDCWire::getWireLine
const WireLine & getWireLine() const
Getter for the wire line represenation of the wire.
Definition: CDCWire.h:190
Belle2::TrackFindingCDC::CDCWire
Class representing a sense wire in the central drift chamber.
Definition: CDCWire.h:60
Belle2::TrackFindingCDC::CDCWireHit
Class representing a hit wire in the central drift chamber.
Definition: CDCWireHit.h:65
Belle2::ModuleParamList
The Module parameter list class.
Definition: ModuleParamList.h:46
Belle2::TrackFindingCDC::CDCTrajectory2D::calcArcLength2D
double calcArcLength2D(const Vector2D &point) const
Calculate the travel distance from the start position of the trajectory.
Definition: CDCTrajectory2D.h:270