simb::MCTrajectory Class Reference

#include <MCTrajectory.h>

Public Types
typedef std::vector< std::pair< TLorentzVector, TLorentzVector > >	list_type
typedef list_type::value_type	value_type
typedef list_type::iterator	iterator
typedef list_type::const_iterator	const_iterator
typedef list_type::reverse_iterator	reverse_iterator
typedef list_type::const_reverse_iterator	const_reverse_iterator
typedef list_type::size_type	size_type
typedef list_type::difference_type	difference_type
typedef std::vector< std::pair< size_t, unsigned char > >	ProcessMap

Public Member Functions
	MCTrajectory ()
	MCTrajectory (const TLorentzVector &vertex, const TLorentzVector &momentum)
const TLorentzVector &	Position (const size_type) const
	The accessor methods described above. More...
const TLorentzVector &	Momentum (const size_type) const
double	X (const size_type i) const
double	Y (const size_type i) const
double	Z (const size_type i) const
double	T (const size_type i) const
double	Px (const size_type i) const
double	Py (const size_type i) const
double	Pz (const size_type i) const
double	E (const size_type i) const
double	TotalLength () const
iterator	begin ()
const_iterator	begin () const
iterator	end ()
const_iterator	end () const
reverse_iterator	rbegin ()
const_reverse_iterator	rbegin () const
reverse_iterator	rend ()
const_reverse_iterator	rend () const
size_type	size () const
bool	empty () const
void	swap (simb::MCTrajectory &other)
void	clear ()
const value_type &	operator[] (const size_type i) const
const value_type &	at (const size_type i) const
void	push_back (value_type const &v)
void	push_back (TLorentzVector const &p, TLorentzVector const &m)
void	Add (TLorentzVector const &p, TLorentzVector const &m)
void	Add (TLorentzVector const &p, TLorentzVector const &m, std::string const &process, bool keepTransportation=false)
unsigned char	ProcessToKey (std::string const &process) const
std::string	KeyToProcess (unsigned char const &key) const
ProcessMap const &	TrajectoryProcesses () const
void	Sparsify (double margin=.1, bool keep_second_to_last=false)

Private Attributes
list_type	ftrajectory
	The list of trajectory points. More...
ProcessMap	fTrajectoryProcess

Friends
std::ostream &	operator<< (std::ostream &output, const MCTrajectory &)

Detailed Description

Definition at line 58 of file MCTrajectory.h.

Member Typedef Documentation

typedef list_type::const_iterator simb::MCTrajectory::const_iterator

Definition at line 66 of file MCTrajectory.h.

typedef list_type::const_reverse_iterator simb::MCTrajectory::const_reverse_iterator

Definition at line 68 of file MCTrajectory.h.

typedef list_type::difference_type simb::MCTrajectory::difference_type

Definition at line 70 of file MCTrajectory.h.

typedef list_type::iterator simb::MCTrajectory::iterator

Definition at line 65 of file MCTrajectory.h.

typedef std::vector< std::pair<TLorentzVector, TLorentzVector> > simb::MCTrajectory::list_type

Some type definitions to make life easier, and to help "hide" the implementation details. (If you're not familiar with STL, you can ignore these definitions.)

Definition at line 63 of file MCTrajectory.h.

typedef std::vector< std::pair<size_t, unsigned char> > simb::MCTrajectory::ProcessMap

Definition at line 71 of file MCTrajectory.h.

typedef list_type::reverse_iterator simb::MCTrajectory::reverse_iterator

Definition at line 67 of file MCTrajectory.h.

typedef list_type::size_type simb::MCTrajectory::size_type

Definition at line 69 of file MCTrajectory.h.

typedef list_type::value_type simb::MCTrajectory::value_type

Definition at line 64 of file MCTrajectory.h.

Constructor & Destructor Documentation

simb::MCTrajectory::MCTrajectory

(

)

Standard constructor: Start with initial position and momentum of the particle.

Definition at line 24 of file MCTrajectory.cxx.

    : ftrajectory()
  {}

simb::MCTrajectory::MCTrajectory	(	const TLorentzVector &	vertex,
		const TLorentzVector &	momentum
	)

Definition at line 29 of file MCTrajectory.cxx.

  {
    ftrajectory.push_back( value_type( position, momentum ) );
  }

Member Function Documentation

void simb::MCTrajectory::Add	(	TLorentzVector const &	p,
		TLorentzVector const &	m
	)

void simb::MCTrajectory::Add	(	TLorentzVector const &	p,
		TLorentzVector const &	m,
		std::string const &	process,
		bool	keepTransportation = false
	)

Definition at line 136 of file MCTrajectory.cxx.

  {
    // add the the momentum and position, then get the location of the added
    // bits to store the process
    this->push_back(p, m);
    
    size_t insertLoc = ftrajectory.size() - 1;
    
    auto key = this->ProcessToKey(process);
    
    // only add a process to the list if it is defined, ie one of the values
    // allowed in the ProcessToKey() method
    //
    // Also, keep 10 (transportation) if the flag allows
    if(key > 0 && (key != 10 || keepTransportation))
      fTrajectoryProcess.push_back(std::make_pair(insertLoc, key));
    
    return;
  }

const simb::MCTrajectory::value_type & simb::MCTrajectory::at ( const size_type  i ) const

inline

Definition at line 175 of file MCTrajectory.h.

                                                                                             { return ftrajectory.at(i); }

simb::MCTrajectory::iterator simb::MCTrajectory::begin ( )

inline

Standard STL methods, to make this class look like an STL map. Again, if you don't know STL, you can just ignore these methods.

Definition at line 158 of file MCTrajectory.h.

{ return ftrajectory.begin();  }

simb::MCTrajectory::const_iterator simb::MCTrajectory::begin ( ) const

inline

Definition at line 159 of file MCTrajectory.h.

{ return ftrajectory.begin();  }

void simb::MCTrajectory::clear ( )

inline

Definition at line 168 of file MCTrajectory.h.

{ ftrajectory.clear();         }

double simb::MCTrajectory::E ( const size_type  i ) const

inline

Definition at line 156 of file MCTrajectory.h.

{ return Momentum(i).E();      }

bool simb::MCTrajectory::empty ( ) const

inline

Definition at line 167 of file MCTrajectory.h.

{ return ftrajectory.empty();  }

simb::MCTrajectory::iterator simb::MCTrajectory::end ( void  )

inline

Definition at line 160 of file MCTrajectory.h.

{ return ftrajectory.end();    }

simb::MCTrajectory::const_iterator simb::MCTrajectory::end ( void  ) const

inline

Definition at line 161 of file MCTrajectory.h.

{ return ftrajectory.end();    }

std::string simb::MCTrajectory::KeyToProcess

(

unsigned char const &

key

)

const

Definition at line 117 of file MCTrajectory.cxx.

  {
    std::string process("Unknown");
    
    if     (key == 1) process = "hadElastic";
    else if(key == 2) process = "pi-Inelastic";
    else if(key == 3) process = "pi+Inelastic";
    else if(key == 4) process = "kaon-Inelastic";
    else if(key == 5) process = "kaon+Inelastic";
    else if(key == 6) process = "protonInelastic";
    else if(key == 7) process = "neutronInelastic";
    else if(key == 8) process = "CoulombScat";
    else if(key == 9) process = "nCapture";
    else if(key == 10) process = "Transportation";
    
    return process;
  }

const TLorentzVector & simb::MCTrajectory::Momentum

(

const size_type

index

)

const

Definition at line 44 of file MCTrajectory.cxx.

  {
    const_iterator i = ftrajectory.begin();
    std::advance(i,index);
    return (*i).second;
  }

const simb::MCTrajectory::value_type & simb::MCTrajectory::operator[] ( const size_type  i ) const

inline

Definition at line 172 of file MCTrajectory.h.

                                                                                                { return ftrajectory[i]; }

const TLorentzVector & simb::MCTrajectory::Position

(

const size_type

index

)

const

The accessor methods described above.

Definition at line 36 of file MCTrajectory.cxx.

  {
    const_iterator i = ftrajectory.begin();
    std::advance(i,index);
    return (*i).first;
  }

unsigned char simb::MCTrajectory::ProcessToKey

(

std::string const &

process

)

const

Definition at line 97 of file MCTrajectory.cxx.

  {
    unsigned char key = 0;
    
    if     (process.compare("hadElastic")       == 0) key = 1;
    else if(process.compare("pi-Inelastic")     == 0) key = 2;
    else if(process.compare("pi+Inelastic")     == 0) key = 3;
    else if(process.compare("kaon-Inelastic")   == 0) key = 4;
    else if(process.compare("kaon+Inelastic")   == 0) key = 5;
    else if(process.compare("protonInelastic")  == 0) key = 6;
    else if(process.compare("neutronInelastic") == 0) key = 7;
    else if(process.compare("CoulombScat")      == 0) key = 8;
    else if(process.compare("nCapture")         == 0) key = 9;
    else if(process.compare("Transportation") == 0)
      key = 10;
    
    return key;
  }

void simb::MCTrajectory::push_back

(

value_type const &

v

)

The only "set" methods for this class; once you've added a trajectory point, you can't take it back.

void simb::MCTrajectory::push_back	(	TLorentzVector const &	p,
		TLorentzVector const &	m
	)

double simb::MCTrajectory::Px ( const size_type  i ) const

inline

Definition at line 153 of file MCTrajectory.h.

{ return Momentum(i).Px();     }

double simb::MCTrajectory::Py ( const size_type  i ) const

inline

Definition at line 154 of file MCTrajectory.h.

{ return Momentum(i).Py();     }

double simb::MCTrajectory::Pz ( const size_type  i ) const

inline

Definition at line 155 of file MCTrajectory.h.

{ return Momentum(i).Pz();     }

simb::MCTrajectory::reverse_iterator simb::MCTrajectory::rbegin ( )

inline

Definition at line 162 of file MCTrajectory.h.

{ return ftrajectory.rbegin(); }

simb::MCTrajectory::const_reverse_iterator simb::MCTrajectory::rbegin ( ) const

inline

Definition at line 163 of file MCTrajectory.h.

{ return ftrajectory.rbegin(); }

simb::MCTrajectory::reverse_iterator simb::MCTrajectory::rend ( )

inline

Definition at line 164 of file MCTrajectory.h.

{ return ftrajectory.rend();   }

simb::MCTrajectory::const_reverse_iterator simb::MCTrajectory::rend ( ) const

inline

Definition at line 165 of file MCTrajectory.h.

{ return ftrajectory.rend();   }

simb::MCTrajectory::size_type simb::MCTrajectory::size ( void  ) const

inline

Definition at line 166 of file MCTrajectory.h.

{ return ftrajectory.size();   }

void simb::MCTrajectory::Sparsify	(	double	margin = .1,
		bool	keep_second_to_last = false
	)

Remove points from trajectory. Straight line interpolation between the remaining points will pass no further than margin from removed points.

Definition at line 160 of file MCTrajectory.cxx.

  {
    // This is a divide-and-conquer algorithm. If the straight line between two
    // points is close enough to all the intermediate points, then just keep
    // the endpoints. Otherwise, divide the range in two and try again.

    // We keep the ranges that need checking in "toCheck". If a range is good
    // as-is, we put just the starting point in "done". The end-point will be
    // taken care of by the next range.

    // Need at least three points to think of removing one
    // D.R. -- let's up this to four points before we start removing
    //      -- this is helpful when retrieving the energy of the particle prior
    //      -- to a final interaction : (Start, p1, ..., p_(n-1), End)
    if(size() <= 3 && keep_second_to_last) return;
    else if(size() <= 2) return;

    // Deal in terms of distance-squared to save some sqrts
    margin *= margin;

    // Deque because we add things still to check on the end, and pop things
    // we've checked from the front.
    std::deque<std::pair<int, int> > toCheck;
    // Start off by trying to replace the whole trajectory with just the
    // endpoints.
    toCheck.push_back(std::make_pair(0, size()-1));

    // use a std::set to keep track of which indices of points we want to
    // keep because the set does not allow duplicates and it keeps items in
    // order
    std::set<int> done;
    if (keep_second_to_last)
      done.insert(size()-2); // -- D.R. store the penultimate point

    while(!toCheck.empty()){
      const int loIdx = toCheck.front().first;
      const int hiIdx = toCheck.front().second;
      toCheck.pop_front();

      // Should never have been given a degenerate range
      if(hiIdx < loIdx+2)
        throw cet::exception("MCTrajectory") << "Degnerate range in Sparsify method";

      const TVector3 loVec = at(loIdx).first.Vect();
      const TVector3 hiVec = at(hiIdx).first.Vect();

      const TVector3 dir = (hiVec-loVec).Unit();

      // Are all the points in between close enough?
      bool ok = true;
      for(int i = loIdx+1; i < hiIdx; ++i){
        const TVector3 toHere = at(i).first.Vect()-loVec;
          // Perpendicular distance^2 from the line joining loVec to hiVec
        const double impact = (toHere-dir.Dot(toHere)*dir).Mag2();
        if(impact > margin){ok = false; break;}
      }
      
      if(ok){
        // These points adequately represent this range
        done.insert(loIdx);
      }
      else{
        // Split in half
        const int midIdx = (loIdx+hiIdx)/2;
        // Should never have a range this small
        if(midIdx == loIdx)
          throw cet::exception("MCTrajectory") << "Midpoint in sparsification is same as lowpoint";
        if(midIdx == hiIdx)
          throw cet::exception("MCTrajectory") << "Midpoint in sparsification is same as hipoint";

        // The range can be small enough that upon splitting, the new ranges
        // are degenerate, and should just be written out straight away. Check
        // for those cases.

        if(midIdx == loIdx+1){
          done.insert(loIdx);
        }
        else{
          toCheck.push_back(std::make_pair(loIdx, midIdx));
        }

        if(midIdx == hiIdx-1){
          done.insert(midIdx);
        }
        else{
          toCheck.push_back(std::make_pair(midIdx, hiIdx));
        }
      }
    } // end while

    // now make sure we have not left out any of the indices where interesting
    // processes were recorded
    std::map< size_t, unsigned char> processMap;
    for(auto itr : fTrajectoryProcess){
      done.insert(itr.first);
      processMap[itr.first] = itr.second;
    }

    // Look up the trajectory points at the stored indices, write them to a new
    // trajectory
    const unsigned int I = done.size();
    list_type newTraj;
    newTraj.reserve(I+1);

    // make a new process map as well based on these points
    ProcessMap newProcMap;
    
    for(auto idx : done){
      newTraj.push_back(at(idx));
      if(processMap.count(idx) > 0){
        newProcMap.push_back(std::make_pair(newTraj.size() - 1,
                                            processMap.find(idx)->second)
                             );
      }
    }

    // Remember to add the very last point in if it hasn't already been added
    if(done.count(ftrajectory.size() - 1) < 1) newTraj.push_back(*rbegin());

    // Replace trajectory and fTrajectoryProcess with new versions
    std::swap(ftrajectory,        newTraj   );
    std::swap(fTrajectoryProcess, newProcMap);
    
    return;
  }

void simb::MCTrajectory::swap ( simb::MCTrajectory &  other )

inline

Definition at line 169 of file MCTrajectory.h.

                                                                                { ftrajectory.swap( other.ftrajectory ); }

double simb::MCTrajectory::T ( const size_type  i ) const

inline

Definition at line 152 of file MCTrajectory.h.

{ return Position(i).T();      }

double simb::MCTrajectory::TotalLength

(

)

const

Definition at line 52 of file MCTrajectory.cxx.

  {
    const int N = size();
    if(N < 2) return 0;

    // We take the sum of the straight lines between the trajectory points
    double dist = 0;
    for(int n = 0; n < N-1; ++n){
      dist += (Position(n+1)-Position(n)).Vect().Mag();
    }

    return dist;
  }

simb::MCTrajectory::ProcessMap const & simb::MCTrajectory::TrajectoryProcesses ( ) const

inline

Definition at line 188 of file MCTrajectory.h.

{ return fTrajectoryProcess; }

double simb::MCTrajectory::X ( const size_type  i ) const

inline

Definition at line 149 of file MCTrajectory.h.

{ return Position(i).X();      }

double simb::MCTrajectory::Y ( const size_type  i ) const

inline

Definition at line 150 of file MCTrajectory.h.

{ return Position(i).Y();      }

double simb::MCTrajectory::Z ( const size_type  i ) const

inline

Definition at line 151 of file MCTrajectory.h.

{ return Position(i).Z();      }

Friends And Related Function Documentation

std::ostream& operator<< ( std::ostream &  output, const MCTrajectory &  list  )

friend

Definition at line 67 of file MCTrajectory.cxx.

  {
    // Determine a field width for the voxel number.
    MCTrajectory::size_type numberOfTrajectories = list.size();
    int numberOfDigits = (int) std::log10( (double) numberOfTrajectories ) + 1;

    // A simple header.
    output.width( numberOfDigits );
    output << "#" << ": < position (x,y,z,t), momentum (Px,Py,Pz,E) >" << std::endl; 

    // Write each trajectory point on a separate line.
    MCTrajectory::size_type nTrajectory = 0;
    for ( MCTrajectory::const_iterator trajectory = list.begin(); trajectory != list.end(); ++trajectory, ++nTrajectory ){
        output.width( numberOfDigits );
        output << nTrajectory << ": "
               << "< (" << (*trajectory).first.X() 
               << "," << (*trajectory).first.Y() 
               << "," << (*trajectory).first.Z() 
               << "," << (*trajectory).first.T() 
               << ") , (" << (*trajectory).second.Px() 
               << "," << (*trajectory).second.Py() 
               << "," << (*trajectory).second.Pz() 
               << "," << (*trajectory).second.E() 
               << ") >" << std::endl;
      }

    return output;
  }

Member Data Documentation

list_type simb::MCTrajectory::ftrajectory

private

The list of trajectory points.

Definition at line 77 of file MCTrajectory.h.

ProcessMap simb::MCTrajectory::fTrajectoryProcess

private

map of the scattering process to index in ftrajectory for a given point

Definition at line 78 of file MCTrajectory.h.

---

The documentation for this class was generated from the following files:

• nusimdata/nusimdata/SimulationBase/MCTrajectory.h
• nusimdata/nusimdata/SimulationBase/MCTrajectory.cxx