LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
MCTrajectory.cxx
Go to the documentation of this file.
1 
8 #include "cetlib_except/exception.h"
9 
11 
12 #include <TLorentzVector.h>
13 
14 #include <cmath>
15 #include <deque>
16 #include <iterator>
17 #include <vector>
18 #include <set>
19 #include <map>
20 
21 namespace simb {
22 
23  // Nothing special need be done for the default constructor or destructor.
25  : ftrajectory()
26  {}
27 
28  //----------------------------------------------------------------------------
29  MCTrajectory::MCTrajectory( const TLorentzVector& position,
30  const TLorentzVector& momentum )
31  {
32  ftrajectory.push_back( value_type( position, momentum ) );
33  }
34 
35  //----------------------------------------------------------------------------
36  const TLorentzVector& MCTrajectory::Position( const size_type index ) const
37  {
38  const_iterator i = ftrajectory.begin();
39  std::advance(i,index);
40  return (*i).first;
41  }
42 
43  //----------------------------------------------------------------------------
44  const TLorentzVector& MCTrajectory::Momentum( const size_type index ) const
45  {
46  const_iterator i = ftrajectory.begin();
47  std::advance(i,index);
48  return (*i).second;
49  }
50 
51  //----------------------------------------------------------------------------
53  {
54  const int N = size();
55  if(N < 2) return 0;
56 
57  // We take the sum of the straight lines between the trajectory points
58  double dist = 0;
59  for(int n = 0; n < N-1; ++n){
60  dist += (Position(n+1)-Position(n)).Vect().Mag();
61  }
62 
63  return dist;
64  }
65 
66  //----------------------------------------------------------------------------
67  std::ostream& operator<< ( std::ostream& output, const MCTrajectory& list )
68  {
69  // Determine a field width for the voxel number.
70  MCTrajectory::size_type numberOfTrajectories = list.size();
71  int numberOfDigits = (int) std::log10( (double) numberOfTrajectories ) + 1;
72 
73  // A simple header.
74  output.width( numberOfDigits );
75  output << "#" << ": < position (x,y,z,t), momentum (Px,Py,Pz,E) >" << std::endl;
76 
77  // Write each trajectory point on a separate line.
78  MCTrajectory::size_type nTrajectory = 0;
79  for ( MCTrajectory::const_iterator trajectory = list.begin(); trajectory != list.end(); ++trajectory, ++nTrajectory ){
80  output.width( numberOfDigits );
81  output << nTrajectory << ": "
82  << "< (" << (*trajectory).first.X()
83  << "," << (*trajectory).first.Y()
84  << "," << (*trajectory).first.Z()
85  << "," << (*trajectory).first.T()
86  << ") , (" << (*trajectory).second.Px()
87  << "," << (*trajectory).second.Py()
88  << "," << (*trajectory).second.Pz()
89  << "," << (*trajectory).second.E()
90  << ") >" << std::endl;
91  }
92 
93  return output;
94  }
95 
96  //----------------------------------------------------------------------------
97  unsigned char MCTrajectory::ProcessToKey(std::string const& process) const
98  {
99  unsigned char key = 0;
100 
101  if (process.compare("hadElastic") == 0) key = 1;
102  else if(process.compare("pi-Inelastic") == 0) key = 2;
103  else if(process.compare("pi+Inelastic") == 0) key = 3;
104  else if(process.compare("kaon-Inelastic") == 0) key = 4;
105  else if(process.compare("kaon+Inelastic") == 0) key = 5;
106  else if(process.compare("protonInelastic") == 0) key = 6;
107  else if(process.compare("neutronInelastic") == 0) key = 7;
108  else if(process.compare("CoulombScat") == 0) key = 8;
109  else if(process.compare("nCapture") == 0) key = 9;
110  else if(process.compare("Transportation") == 0)
111  key = 10;
112 
113  return key;
114  }
115 
116  //----------------------------------------------------------------------------
117  std::string MCTrajectory::KeyToProcess(unsigned char const& key) const
118  {
119  std::string process("Unknown");
120 
121  if (key == 1) process = "hadElastic";
122  else if(key == 2) process = "pi-Inelastic";
123  else if(key == 3) process = "pi+Inelastic";
124  else if(key == 4) process = "kaon-Inelastic";
125  else if(key == 5) process = "kaon+Inelastic";
126  else if(key == 6) process = "protonInelastic";
127  else if(key == 7) process = "neutronInelastic";
128  else if(key == 8) process = "CoulombScat";
129  else if(key == 9) process = "nCapture";
130  else if(key == 10) process = "Transportation";
131 
132  return process;
133  }
134 
135  //----------------------------------------------------------------------------
136  void MCTrajectory::Add(TLorentzVector const& p,
137  TLorentzVector const& m,
138  std::string const& process,
139  bool keepTransportation)
140  {
141  // add the the momentum and position, then get the location of the added
142  // bits to store the process
143  this->push_back(p, m);
144 
145  size_t insertLoc = ftrajectory.size() - 1;
146 
147  auto key = this->ProcessToKey(process);
148 
149  // only add a process to the list if it is defined, ie one of the values
150  // allowed in the ProcessToKey() method
151  //
152  // Also, keep 10 (transportation) if the flag allows
153  if(key > 0 && (key != 10 || keepTransportation))
154  fTrajectoryProcess.push_back(std::make_pair(insertLoc, key));
155 
156  return;
157  }
158 
159  //----------------------------------------------------------------------------
160  void MCTrajectory::Sparsify(double margin, bool keep_second_to_last)
161  {
162  // This is a divide-and-conquer algorithm. If the straight line between two
163  // points is close enough to all the intermediate points, then just keep
164  // the endpoints. Otherwise, divide the range in two and try again.
165 
166  // We keep the ranges that need checking in "toCheck". If a range is good
167  // as-is, we put just the starting point in "done". The end-point will be
168  // taken care of by the next range.
169 
170  // Need at least three points to think of removing one
171  // D.R. -- let's up this to four points before we start removing
172  // -- this is helpful when retrieving the energy of the particle prior
173  // -- to a final interaction : (Start, p1, ..., p_(n-1), End)
174  if(size() <= 3 && keep_second_to_last) return;
175  else if(size() <= 2) return;
176 
177  // Deal in terms of distance-squared to save some sqrts
178  margin *= margin;
179 
180  // Deque because we add things still to check on the end, and pop things
181  // we've checked from the front.
182  std::deque<std::pair<int, int> > toCheck;
183  // Start off by trying to replace the whole trajectory with just the
184  // endpoints.
185  toCheck.push_back(std::make_pair(0, size()-1));
186 
187  // use a std::set to keep track of which indices of points we want to
188  // keep because the set does not allow duplicates and it keeps items in
189  // order
190  std::set<int> done;
191  if (keep_second_to_last)
192  done.insert(size()-2); // -- D.R. store the penultimate point
193 
194  while(!toCheck.empty()){
195  const int loIdx = toCheck.front().first;
196  const int hiIdx = toCheck.front().second;
197  toCheck.pop_front();
198 
199  // Should never have been given a degenerate range
200  if(hiIdx < loIdx+2)
201  throw cet::exception("MCTrajectory") << "Degnerate range in Sparsify method";
202 
203  const TVector3 loVec = at(loIdx).first.Vect();
204  const TVector3 hiVec = at(hiIdx).first.Vect();
205 
206  const TVector3 dir = (hiVec-loVec).Unit();
207 
208  // Are all the points in between close enough?
209  bool ok = true;
210  for(int i = loIdx+1; i < hiIdx; ++i){
211  const TVector3 toHere = at(i).first.Vect()-loVec;
212  // Perpendicular distance^2 from the line joining loVec to hiVec
213  const double impact = (toHere-dir.Dot(toHere)*dir).Mag2();
214  if(impact > margin){ok = false; break;}
215  }
216 
217  if(ok){
218  // These points adequately represent this range
219  done.insert(loIdx);
220  }
221  else{
222  // Split in half
223  const int midIdx = (loIdx+hiIdx)/2;
224  // Should never have a range this small
225  if(midIdx == loIdx)
226  throw cet::exception("MCTrajectory") << "Midpoint in sparsification is same as lowpoint";
227  if(midIdx == hiIdx)
228  throw cet::exception("MCTrajectory") << "Midpoint in sparsification is same as hipoint";
229 
230  // The range can be small enough that upon splitting, the new ranges
231  // are degenerate, and should just be written out straight away. Check
232  // for those cases.
233 
234  if(midIdx == loIdx+1){
235  done.insert(loIdx);
236  }
237  else{
238  toCheck.push_back(std::make_pair(loIdx, midIdx));
239  }
240 
241  if(midIdx == hiIdx-1){
242  done.insert(midIdx);
243  }
244  else{
245  toCheck.push_back(std::make_pair(midIdx, hiIdx));
246  }
247  }
248  } // end while
249 
250  // now make sure we have not left out any of the indices where interesting
251  // processes were recorded
252  std::map< size_t, unsigned char> processMap;
253  for(auto itr : fTrajectoryProcess){
254  done.insert(itr.first);
255  processMap[itr.first] = itr.second;
256  }
257 
258  // Look up the trajectory points at the stored indices, write them to a new
259  // trajectory
260  const unsigned int I = done.size();
261  list_type newTraj;
262  newTraj.reserve(I+1);
263 
264  // make a new process map as well based on these points
265  ProcessMap newProcMap;
266 
267  for(auto idx : done){
268  newTraj.push_back(at(idx));
269  if(processMap.count(idx) > 0){
270  newProcMap.push_back(std::make_pair(newTraj.size() - 1,
271  processMap.find(idx)->second)
272  );
273  }
274  }
275 
276  // Remember to add the very last point in if it hasn't already been added
277  if(done.count(ftrajectory.size() - 1) < 1) newTraj.push_back(*rbegin());
278 
279  // Replace trajectory and fTrajectoryProcess with new versions
280  std::swap(ftrajectory, newTraj );
281  std::swap(fTrajectoryProcess, newProcMap);
282 
283  return;
284  }
285 
286 } // namespace sim
void Add(TLorentzVector const &p, TLorentzVector const &m)
Trajectory class.
void push_back(value_type const &v)
reverse_iterator rbegin()
Definition: MCTrajectory.h:162
const value_type & at(const size_type i) const
Definition: MCTrajectory.h:175
std::string KeyToProcess(unsigned char const &key) const
std::vector< std::pair< TLorentzVector, TLorentzVector > > list_type
Definition: MCTrajectory.h:63
unsigned char ProcessToKey(std::string const &process) const
list_type::value_type value_type
Definition: MCTrajectory.h:64
list_type::size_type size_type
Definition: MCTrajectory.h:69
list_type::const_iterator const_iterator
Definition: MCTrajectory.h:66
ProcessMap fTrajectoryProcess
Definition: MCTrajectory.h:78
const TLorentzVector & Position(const size_type) const
The accessor methods described above.
ART objects.
std::vector< std::pair< size_t, unsigned char > > ProcessMap
Definition: MCTrajectory.h:71
size_type size() const
Definition: MCTrajectory.h:166
void Sparsify(double margin=.1, bool keep_second_to_last=false)
TDirectory * dir
Definition: macro.C:5
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
double TotalLength() const
Char_t n[5]
const TLorentzVector & Momentum(const size_type) const
friend std::ostream & operator<<(std::ostream &output, const MCTrajectory &)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
list_type ftrajectory
The list of trajectory points.
Definition: MCTrajectory.h:77