Belle II Software  release-05-01-25
MCParticleTrajectory.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Martin Ritter *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <simulation/dataobjects/MCParticleTrajectory.h>
12 #include <TVector3.h>
13 #include <stack>
14 #include <tuple>
15 
16 using namespace Belle2;
17 
18 namespace {
26  void setVector(TVector3& v, const MCTrajectoryPoint& a, const MCTrajectoryPoint& b, bool unit = false)
27  {
28  v.SetXYZ(a.x - b.x, a.y - b.y, a.z - b.z);
29  if (unit) {
30  v *= 1. / v.Mag();
31  }
32  }
33 }
34 
35 void MCParticleTrajectory::simplify(float distance_tolerance)
36 {
37  // cannot simplify anything, return
38  if (distance_tolerance <= 0 || m_points.size() < 3) return;
39 
40  // stack with all segments to be investigated
41  std::stack<std::pair<iterator, iterator>> stack;
42  // push full trajectory on the stack
43  stack.push(make_pair(m_points.begin(), m_points.end() - 1));
44  // next free point: we always want the starting point so start at index 1
45  iterator nextFreePoint = m_points.begin() + 1;
46  // iterators used for the segment inspection
47  iterator firstPoint, splitPoint, finalPoint;
48  // segment direction and vector between segment start and mid point
49  TVector3 n, pa;
50  // investigate all segments until all fulfill the distance requirement
51  while (!stack.empty()) {
52  //Get first and last point
53  std::tie(firstPoint, finalPoint) = stack.top();
54  //Remove segment from stack
55  stack.pop();
56  //Direction of the segment
57  setVector(n, *firstPoint, *finalPoint, true);
58  //Calculate maximum distance of all intermediate points to the segment
59  double maxDistance(0);
60  for (auto nextPoint = firstPoint + 1; nextPoint != finalPoint; ++nextPoint) {
61  //vector from segment start (p) to point (a)
62  setVector(pa, *firstPoint, *nextPoint);
63  //3D distance between point a and line p + x*n
64  const double dist = (pa - (pa * n) * n).Mag();
65  //check if this is the maximum distance so far
66  if (dist > maxDistance) {
67  splitPoint = nextPoint;
68  maxDistance = dist;
69  }
70  }
71  //Are all points close enough? if not split the segment at the largest distance
72  if (maxDistance > distance_tolerance) {
73  //If we split in this order, all points will be in correct order since we
74  //use a stack: last thing put in is first thing to get out. So initially
75  //we inspect everything, if we divide we put in the second part and then
76  //the first and in the next round we take out the first and repeat. So
77  //effectively we look at all segments in an ordered way.
78  stack.push(make_pair(splitPoint, finalPoint));
79  stack.push(make_pair(firstPoint, splitPoint));
80  continue;
81  }
82  //Ok, all points close enough, add the final point to list of points. Due
83  //to the order in which we look at the points in a ordered way we can
84  //replace the points in place
85  *(nextFreePoint++) = *finalPoint;
86  }
87 
88  //Now delete all remaining elements
89  m_points.erase(nextFreePoint, end());
90 }
Belle2::MCParticleTrajectory::iterator
std::vector< MCTrajectoryPoint >::iterator iterator
iterator definition to allow iteration
Definition: MCParticleTrajectory.h:38
Belle2::MCTrajectoryPoint
Small struct to encode a position/momentum without additional overhead.
Definition: MCTrajectoryPoint.h:21
Belle2::MCParticleTrajectory::simplify
void simplify(float distanceTolerance)
Simplify the trajectory using the Ramer-Douglas-Peuker algorithm.
Definition: MCParticleTrajectory.cc:35
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::MCParticleTrajectory::end
iterator end()
return iterator beyond the last point
Definition: MCParticleTrajectory.h:44
Belle2::MCParticleTrajectory::m_points
std::vector< MCTrajectoryPoint > m_points
Collection of points along the trajectory.
Definition: MCParticleTrajectory.h:78