Belle II Software  release-05-01-25
AxialLegendreLeafProcessor.icc.h
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2015 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Bastian Kronenbitter, Thomas Hauth, Viktor Trusov, *
7  * Nils Braun, Oliver Frost *
8  * *
9  * This software is provided "as is" without any warranty. *
10  **************************************************************************/
11 #pragma once
12 #include <tracking/trackFindingCDC/hough/perigee/AxialLegendreLeafProcessor.h>
13 
14 #include <tracking/trackFindingCDC/processing/AxialTrackUtil.h>
15 
16 #include <tracking/trackFindingCDC/hough/perigee/StereoHitContained.h>
17 #include <tracking/trackFindingCDC/hough/perigee/OffOrigin.h>
18 #include <tracking/trackFindingCDC/hough/algorithms/InPhi0ImpactCurvBox.h>
19 #include <tracking/trackFindingCDC/hough/baseelements/WithSharedMark.h>
20 
21 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCTrajectory2D.h>
22 #include <tracking/trackFindingCDC/fitting/CDCObservations2D.h>
23 #include <tracking/trackFindingCDC/fitting/CDCKarimakiFitter.h>
24 #include <tracking/trackFindingCDC/geometry/PerigeeCircle.h>
25 
26 #include <tracking/trackFindingCDC/legendre/precisionFunctions/PrecisionUtil.h>
27 
28 #include <tracking/trackFindingCDC/utilities/StringManipulation.h>
29 
30 #include <framework/core/ModuleParamList.templateDetails.h>
31 
32 namespace Belle2 {
37  namespace TrackFindingCDC {
38 
39  template <class ANode>
41  {
42  // Special post processing looking for more hits
43  // Start off by fitting the items of the leaf with a general trajectory
44  // that may have a distinct impact != 0
46  auto fitPos = EFitPos::c_RecoPos;
47  auto fitVariance = EFitVariance::c_DriftLength;
48  // Other combinations of fit information
49  // EFitPos::c_RLDriftCircle x EFitVariance::(c_Nominal, c_Pseudo, c_Proper)
50  // have been tried, but found to be worse, which is
51  // not intutive. Probably the perfect circle trajectory
52  // is not as good of a model on the full CDC volume.
53  CDCObservations2D observations2D(fitPos, fitVariance);
54 
55  CDCKarimakiFitter fitter;
56  // Tested alternative CDCRiemannFitter with only marginal differences;
57 
58  std::vector<WithSharedMark<CDCRLWireHit>> hits(leaf->begin(), leaf->end());
59  std::sort(hits.begin(), hits.end()); // Hits should be naturally sorted
60  observations2D.appendRange(hits);
61  CDCTrajectory2D trajectory2D = fitter.fit(observations2D);
62  {
63  const double curv = trajectory2D.getCurvature();
64  std::array<DiscreteCurv, 2> curvs = leaf->template getBounds<DiscreteCurv>();
65 
66  float lowerCurv = *(curvs[0]);
67  float upperCurv = *(curvs[1]);
68  if (ESignUtil::common(lowerCurv, upperCurv) * curv < 0) {
69  trajectory2D.reverse();
70  }
71  }
72  trajectory2D.setLocalOrigin(Vector2D(0, 0));
73 
74  // Look for more hits near the found trajectory
76  if (m_param_nRoadSearches > 0) {
77  // Acquire all available items in some parent node for the road search
78  // Somewhat bypasses the logic of the tree walk - Must handle with care!
79  ANode* roadNode = leaf;
80  while (roadNode->getParent() and roadNode->getLevel() > m_param_roadLevel) {
81  roadNode = roadNode->getParent();
82  }
83 
84  for (int iRoadSearch = 0; iRoadSearch < m_param_nRoadSearches; ++iRoadSearch) {
85  // Use a road search to find new hits from the trajectory
86  // hits = this->searchRoad(*roadNode, trajectory2D); // In case you only want the road hits
87  // if (hits.size() < 5) return;
88 
89  // Second version always holding on to the originally found hits
90  int nHitsBefore = hits.size();
91  std::vector<WithSharedMark<CDCRLWireHit>> roadHits = this->searchRoad(*roadNode, trajectory2D);
92  std::sort(roadHits.begin(), roadHits.end());
93  hits.insert(hits.end(), roadHits.begin(), roadHits.end());
94  std::inplace_merge(hits.begin(), hits.begin() + nHitsBefore, hits.end());
95  hits.erase(std::unique(hits.begin(), hits.end()), hits.end());
96 
97  // Update the current fit
98  observations2D.clear();
99  observations2D.appendRange(hits);
100  trajectory2D = fitter.fit(observations2D);
101  trajectory2D.setLocalOrigin(Vector2D(0.0, 0.0));
102  }
103  }
104 
105  // Mark found hit as used and safe them with the trajectory
107  std::vector<const CDCWireHit*> foundWireHits;
108  for (CDCRLWireHit& rlWireHit : hits) {
109  foundWireHits.push_back(&rlWireHit.getWireHit());
110  }
111 
112  AxialTrackUtil::addCandidateFromHits(foundWireHits, m_axialWireHits, m_tracks, true);
113 
114  // Sync up the marks with the used hits
115  for (WithSharedMark<CDCRLWireHit>& markableRLWireHit : leaf->getTree()->getTopNode()) {
116  const AutomatonCell& automatonCell = markableRLWireHit.getWireHit().getAutomatonCell();
117  if (automatonCell.hasTakenFlag() or automatonCell.hasMaskedFlag()) {
118  markableRLWireHit.mark();
119  } else {
120  markableRLWireHit.unmark();
121  }
122  }
123  }
124 
125  template <class ANode>
126  std::vector<WithSharedMark<CDCRLWireHit> >
127  AxialLegendreLeafProcessor<ANode>::searchRoad(const ANode& node, const CDCTrajectory2D& trajectory2D)
128  {
129  PerigeeCircle circle = trajectory2D.getGlobalCircle();
130  Vector2D support = trajectory2D.getGlobalPerigee();
131  const float curv = circle.curvature();
132  const float phi0 = circle.phi0();
133 
134  StereoHitContained<OffOrigin<InPhi0ImpactCurvBox> > hitInPhi0CurvBox(m_param_curlCurv);
135  hitInPhi0CurvBox.setLocalOrigin(support);
136  using RoadHoughBox = StereoHitContained<OffOrigin<InPhi0ImpactCurvBox> >::HoughBox;
137 
138  // Determine a precision that we expect to achieve at the fitted momentum
139  // There certainly is some optimizsation potential here.
140  // Spread in the impact parameter is made available here but is not activated yet.
141  const float levelPrecision = 9.0;
142  // Earlier version
143  // const float levelPrecision = 10.5 - 0.24 * exp(-4.13118 * PrecisionUtil::convertRhoToPt(curv) + 2.74);
144  const float phi0Precision = 3.1415 / std::pow(2.0, levelPrecision + 1.0);
145  const float impactPrecision = 0.0 * std::sqrt(CDCWireHit::c_simpleDriftLengthVariance);
146  const float curvPrecision = 0.15 / std::pow(2.0, levelPrecision);
147 
148  DiscreteCurv::Array curvBounds{{curv - curvPrecision, curv + curvPrecision}};
149  ContinuousImpact::Array impactBounds{{ -impactPrecision, impactPrecision}};
150  DiscretePhi0::Array phi0Bounds{{Vector2D::Phi(phi0 - phi0Precision), Vector2D::Phi(phi0 + phi0Precision)}};
151  RoadHoughBox precisionPhi0CurvBox(DiscretePhi0::getRange(phi0Bounds),
152  ContinuousImpact::getRange(impactBounds),
153  DiscreteCurv::getRange(curvBounds));
154 
155  std::vector<WithSharedMark<CDCRLWireHit>> hitsInPrecisionBox;
156 
157  // Explicitly making a copy here to ensure that we do not change the node content
158  for (WithSharedMark<CDCRLWireHit> markableRLWireHit : node) {
159  // Skip marked hits
160  if (markableRLWireHit.isMarked()) continue;
161  Weight weight = hitInPhi0CurvBox(markableRLWireHit, &precisionPhi0CurvBox);
162  if (not std::isnan(weight)) hitsInPrecisionBox.push_back(markableRLWireHit);
163  }
164 
165  return hitsInPrecisionBox;
166  }
167  }
169 }
170 
171 namespace Belle2 {
176  namespace TrackFindingCDC {
177  template <class ANode>
178  std::vector<typename AxialLegendreLeafProcessor<ANode>::Candidate>
180  {
181  std::vector<Candidate> result;
182  for (const CDCTrack& track : m_tracks) {
183  std::vector<CDCRLWireHit> rlWireHits;
184  for (const CDCRecoHit3D& recoHit3D : track) {
185  rlWireHits.push_back(recoHit3D.getRLWireHit());
186  }
187  CDCTrajectory2D trajectory2D = track.getStartTrajectory3D().getTrajectory2D();
188  result.emplace_back(std::move(trajectory2D), std::move(rlWireHits));
189  }
190  return result;
191  }
192 
193  template <class ANode>
195  const std::string& prefix)
196  {
197  moduleParamList->addParameter(prefixed(prefix, "maxLevel"),
198  m_param_maxLevel,
199  "Level of divisions in the hough space at which a leaf is reached",
200  m_param_maxLevel);
201 
202  moduleParamList->addParameter(prefixed(prefix, "minWeight"),
203  m_param_minWeight,
204  "Minimal exceptable weight of a leaf to be considered",
205  m_param_minWeight);
206 
207  moduleParamList->addParameter(prefixed(prefix, "maxCurv"),
208  m_param_maxCurv,
209  "Maximal curvature of a leaf to be considered",
210  m_param_maxCurv);
211 
212  moduleParamList->addParameter(prefixed(prefix, "curlCurv"),
213  m_param_curlCurv,
214  "Curvature below which hits on both arms of the trajectory are allowed",
215  m_param_curlCurv);
216 
217  moduleParamList->addParameter(prefixed(prefix, "nRoadSearches"),
218  m_param_nRoadSearches,
219  "How often the road search should be performed to find new hits",
220  m_param_nRoadSearches);
221 
222  moduleParamList->addParameter(prefixed(prefix, "roadLevel"),
223  m_param_roadLevel,
224  "Level of the read from which additional hits in the road search can be taken",
225  m_param_roadLevel);
226 
227  moduleParamList->addParameter(prefixed(prefix, "curvResolution"),
228  m_param_curvResolution,
229  "The name of the resolution function to be used. "
230  "Valid values are 'none', 'const', 'basic', 'origin', 'nonOrigin'",
231  m_param_curvResolution);
232  }
233 
234  template <class ANode>
236  {
237  // Setup the requested precision function
238  if (m_param_curvResolution == "none") {
239  m_curvResolution = [](double curv __attribute__((unused))) { return NAN; };
240  } else if (m_param_curvResolution == "const") {
241  m_curvResolution = [](double curv __attribute__((unused))) { return 0.0008; };
242  } else if (m_param_curvResolution == "basic") {
243  m_curvResolution = &PrecisionUtil::getBasicCurvPrecision;
244  } else if (m_param_curvResolution == "origin") {
245  m_curvResolution = &PrecisionUtil::getOriginCurvPrecision;
246  } else if (m_param_curvResolution == "nonOrigin") {
247  m_curvResolution = &PrecisionUtil::getNonOriginCurvPrecision;
248  } else {
249  B2WARNING("Unknown curvature resolution function " << m_param_curvResolution);
250  m_curvResolution = [](double curv __attribute__((unused))) { return NAN; };
251  }
252  }
253  }
255 }
Belle2::TrackFindingCDC::ESignUtil::common
static ESign common(ESign n1, ESign n2)
Check if two values have a common sign.
Definition: ESign.h:67
Belle2::TrackFindingCDC::CDCRecoHit3D
Class representing a three dimensional reconstructed hit.
Definition: CDCRecoHit3D.h:62
Belle2::TrackFindingCDC::ContinuousValue::getRange
static std::array< This, 2 > getRange(const Array &bounds)
Extract the range from an array providing the discrete values.
Definition: ContinuousValue.h:50
Belle2::TrackFindingCDC::CDCTrack
Class representing a sequence of three dimensional reconstructed hits.
Definition: CDCTrack.h:51
Belle2::TrackFindingCDC::PrecisionUtil::getOriginCurvPrecision
static double getOriginCurvPrecision(double curv)
Function which estimates desired curvature resolution of quadtree node in the given pt region paramet...
Definition: PrecisionUtil.h:75
Belle2::TrackFindingCDC::ContinuousValue::Array
std::array< T, 2 > Array
Mock array type to be a drop in replacement for the discrete values.
Definition: ContinuousValue.h:46
Belle2::TrackFindingCDC::CDCTrajectory2D::reverse
void reverse()
Reverses the trajectory in place.
Definition: CDCTrajectory2D.cc:97
Belle2::TrackFindingCDC::CDCKarimakiFitter
Class implementing the fitter using Karimakis method.
Definition: CDCKarimakiFitter.h:33
Belle2::TrackFindingCDC::CDCTrajectory2D::setLocalOrigin
double setLocalOrigin(const Vector2D &localOrigin)
Setter for the origin of the local coordinate system.
Definition: CDCTrajectory2D.cc:357
Belle2::TrackFindingCDC::CDCTrajectory2D::getCurvature
double getCurvature() const
Getter for the curvature as seen from the xy projection.
Definition: CDCTrajectory2D.h:432
Belle2::TrackFindingCDC::AxialLegendreLeafProcessor::searchRoad
std::vector< WithSharedMark< CDCRLWireHit > > searchRoad(const ANode &node, const CDCTrajectory2D &trajectory2D)
Look for more hits near a ftted trajectory from hits available in the give node.
Definition: AxialLegendreLeafProcessor.icc.h:136
Belle2::TrackFindingCDC::PrecisionUtil::getBasicCurvPrecision
static double getBasicCurvPrecision(double __attribute__((unused)) curv)
Basic function to estimate the curvature precision Takes a curvature value and returns a width that.
Definition: PrecisionUtil.h:59
Belle2::TrackFindingCDC::CDCTrajectory2D
Particle trajectory as it is seen in xy projection represented as a circle.
Definition: CDCTrajectory2D.h:46
Belle2::TrackFindingCDC::DiscreteValue::getRange
static std::array< This, 2 > getRange(Array &values)
Extract the range from an array providing the discrete values.
Definition: DiscreteValue.h:69
Belle2::TrackFindingCDC::PrecisionUtil::getNonOriginCurvPrecision
static double getNonOriginCurvPrecision(double curv)
Function which estimates desired curvature resolution of quadtree node in the given pt region paramet...
Definition: PrecisionUtil.h:100
Belle2::TrackFindingCDC::Vector2D::Phi
static Vector2D Phi(const double phi)
Constucts a unit vector with azimuth angle equal to phi.
Definition: Vector2D.h:73
Belle2::TrackFindingCDC::CDCObservations2D
Class serving as a storage of observed drift circles to present to the Riemann fitter.
Definition: CDCObservations2D.h:53
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::AxialTrackUtil::addCandidateFromHits
static void addCandidateFromHits(const std::vector< const CDCWireHit * > &foundAxialWireHits, const std::vector< const CDCWireHit * > &allAxialWireHits, std::vector< CDCTrack > &axialTracks, bool withPostprocessing=true)
Create CDCTrack using CDCWireHit hits and store it in the list. Then call the postprocessing on it.
Definition: AxialTrackUtil.cc:26
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TrackFindingCDC::AxialLegendreLeafProcessor::processLeaf
void processLeaf(ANode *leaf)
A valuable leaf has been found in the hough tree walk.
Definition: AxialLegendreLeafProcessor.icc.h:49
Belle2::TrackFindingCDC::AxialLegendreLeafProcessor::getCandidates
std::vector< Candidate > getCandidates() const
Getter for the candidates structure still used in some tests.
Definition: AxialLegendreLeafProcessor.icc.h:179
Belle2::TrackFindingCDC::AxialLegendreLeafProcessor::exposeParameters
void exposeParameters(ModuleParamList *moduleParamList, const std::string &prefix)
Expose the parameters as a module parameter list.
Definition: AxialLegendreLeafProcessor.icc.h:194
Belle2::TrackFindingCDC::AxialLegendreLeafProcessor::beginWalk
void beginWalk()
Function to notify the leaf processor about changes in parameters before a new walk.
Definition: AxialLegendreLeafProcessor.icc.h:235
Belle2::ModuleParamList
The Module parameter list class.
Definition: ModuleParamList.h:46
Belle2::TrackFindingCDC::CDCWireHit::c_simpleDriftLengthVariance
static constexpr const double c_simpleDriftLengthVariance
A default value for the drift length variance if no variance from the drift length translation is ava...
Definition: CDCWireHit.h:74
Belle2::TrackFindingCDC::StereoHitContained
Predicate class to check for the containment of axial and stereo hits in some hough space part.
Definition: StereoHitContained.h:45
Belle2::TrackFindingCDC::DiscreteValue::Array
std::vector< T > Array
The type of the array which contains the underlying values.
Definition: DiscreteValue.h:65