MCTrajectory.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file MCTrajectory.cxx
3 /// \brief Container of trajectory information for a particle
4 ///
5 /// \author seligman@nevis.columbia.edu
6 ////////////////////////////////////////////////////////////////////////
7 
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 
109  return key;
110  }
111 
112  //----------------------------------------------------------------------------
113  std::string MCTrajectory::KeyToProcess(unsigned char const& key) const
114  {
115  std::string process("Unknown");
116 
117  if (key == 1) process = "hadElastic";
118  else if(key == 2) process = "pi-Inelastic";
119  else if(key == 3) process = "pi+Inelastic";
120  else if(key == 4) process = "kaon-Inelastic";
121  else if(key == 5) process = "kaon+Inelastic";
122  else if(key == 6) process = "protonInelastic";
123  else if(key == 7) process = "neutronInelastic";
124 
125  return process;
126  }
127 
128  //----------------------------------------------------------------------------
129  void MCTrajectory::Add(TLorentzVector const& p,
130  TLorentzVector const& m,
131  std::string const& process )
132  {
133  // add the the momentum and position, then get the location of the added
134  // bits to store the process
135  this->push_back(p, m);
136 
137  size_t insertLoc = ftrajectory.size() - 1;
138 
139  auto key = this->ProcessToKey(process);
140 
141  // only add a process to the list if it is defined, ie one of the values
142  // allowed in the ProcessToKey() method
143  if(key > 0)
144  fTrajectoryProcess.push_back(std::make_pair(insertLoc, key));
145 
146  return;
147  }
148 
149  //----------------------------------------------------------------------------
150  void MCTrajectory::Sparsify(double margin)
151  {
152  // This is a divide-and-conquer algorithm. If the straight line between two
153  // points is close enough to all the intermediate points, then just keep
154  // the endpoints. Otherwise, divide the range in two and try again.
155 
156  // We keep the ranges that need checking in "toCheck". If a range is good
157  // as-is, we put just the starting point in "done". The end-point will be
158  // taken care of by the next range.
159 
160  // Need at least three points to think of removing one
161  if(size() <= 2) return;
162 
163  // Deal in terms of distance-squared to save some sqrts
164  margin *= margin;
165 
166  // Deque because we add things still to check on the end, and pop things
167  // we've checked from the front.
168  std::deque<std::pair<int, int> > toCheck;
169  // Start off by trying to replace the whole trajectory with just the
170  // endpoints.
171  toCheck.push_back(std::make_pair(0, size()-1));
172 
173  // use a std::set to keep track of which indices of points we want to
174  // keep because the set does not allow duplicates and it keeps items in
175  // order
176  std::set<int> done;
177 
178  while(!toCheck.empty()){
179  const int loIdx = toCheck.front().first;
180  const int hiIdx = toCheck.front().second;
181  toCheck.pop_front();
182 
183  // Should never have been given a degenerate range
184  if(hiIdx < loIdx+2)
185  throw cet::exception("MCTrajectory") << "Degnerate range in Sparsify method";
186 
187  const TVector3 loVec = at(loIdx).first.Vect();
188  const TVector3 hiVec = at(hiIdx).first.Vect();
189 
190  const TVector3 dir = (hiVec-loVec).Unit();
191 
192  // Are all the points in between close enough?
193  bool ok = true;
194  for(int i = loIdx+1; i < hiIdx; ++i){
195  const TVector3 toHere = at(i).first.Vect()-loVec;
196  // Perpendicular distance^2 from the line joining loVec to hiVec
197  const double impact = (toHere-dir.Dot(toHere)*dir).Mag2();
198  if(impact > margin){ok = false; break;}
199  }
200 
201  if(ok){
202  // These points adequately represent this range
203  done.insert(loIdx);
204  }
205  else{
206  // Split in half
207  const int midIdx = (loIdx+hiIdx)/2;
208  // Should never have a range this small
209  if(midIdx == loIdx)
210  throw cet::exception("MCTrajectory") << "Midpoint in sparsification is same as lowpoint";
211  if(midIdx == hiIdx)
212  throw cet::exception("MCTrajectory") << "Midpoint in sparsification is same as hipoint";
213 
214  // The range can be small enough that upon splitting, the new ranges
215  // are degenerate, and should just be written out straight away. Check
216  // for those cases.
217 
218  if(midIdx == loIdx+1){
219  done.insert(loIdx);
220  }
221  else{
222  toCheck.push_back(std::make_pair(loIdx, midIdx));
223  }
224 
225  if(midIdx == hiIdx-1){
226  done.insert(midIdx);
227  }
228  else{
229  toCheck.push_back(std::make_pair(midIdx, hiIdx));
230  }
231  }
232  } // end while
233 
234  // now make sure we have not left out any of the indices where interesting
235  // processes were recorded
236  std::map< size_t, unsigned char> processMap;
237  for(auto itr : fTrajectoryProcess){
238  done.insert(itr.first);
239  processMap[itr.first] = itr.second;
240  }
241 
242  // Look up the trajectory points at the stored indices, write them to a new
243  // trajectory
244  const unsigned int I = done.size();
245  list_type newTraj;
246  newTraj.reserve(I+1);
247 
248  // make a new process map as well based on these points
249  ProcessMap newProcMap;
250 
251  for(auto idx : done){
252  newTraj.push_back(at(idx));
253  if(processMap.count(idx) > 0){
254  newProcMap.push_back(std::make_pair(newTraj.size() - 1,
255  processMap.find(idx)->second)
256  );
257  }
258  }
259 
260  // Remember to add the very last point in if it hasn't already been added
261  if(done.count(ftrajectory.size() - 1) < 1) newTraj.push_back(*rbegin());
262 
263  // Replace trajectory and fTrajectoryProcess with new versions
264  std::swap(ftrajectory, newTraj );
265  std::swap(fTrajectoryProcess, newProcMap);
266 
267  return;
268  }
269 
270 } // namespace sim
void Add(TLorentzVector const &p, TLorentzVector const &m)
ofstream output
Trajectory class.
void push_back(value_type const &v)
reverse_iterator rbegin()
Definition: MCTrajectory.h:161
const char * p
Definition: xmltok.h:285
const value_type & at(const size_type i) const
Definition: MCTrajectory.h:174
std::string KeyToProcess(unsigned char const &key) const
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::vector< std::pair< TLorentzVector, TLorentzVector > > list_type
Definition: MCTrajectory.h:63
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
unsigned char ProcessToKey(std::string const &process) const
list_type::value_type value_type
Definition: MCTrajectory.h:64
double dist
Definition: runWimpSim.h:113
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
void Sparsify(double margin=.1)
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
TVector3 Unit() const
const TLorentzVector & Position(const size_type) const
The accessor methods described above.
This class describes a particle created in the detector Monte Carlo simulation.
std::vector< std::pair< size_t, unsigned char > > ProcessMap
Definition: MCTrajectory.h:71
float Mag2() const
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
size_type size() const
Definition: MCTrajectory.h:165
T log10(T number)
Definition: d0nt_math.hpp:120
TDirectory * dir
Definition: macro.C:5
double TotalLength() const
const TLorentzVector & Momentum(const size_type) const
friend std::ostream & operator<<(std::ostream &output, const MCTrajectory &)
TVector3 Vect() const
list_type ftrajectory
The list of trajectory points.
Definition: MCTrajectory.h:77