ParticleNavigator.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Particle list in DetSim contains Monte Carlo particle information.
3 ///
4 /// \author seligman@nevis.columbia.edu
5 /// \date
6 ////////////////////////////////////////////////////////////////////////
7 ///
8 /// Although there's nothing in the following class that assumes
9 /// units, the standard for NOvASoft is that distances are in cm, and
10 /// energies are in GeV.
11 ////////////////////////////////////////////////////////////////////////
13 
14 #include <cmath>
15 #include <iostream>
16 #include <iterator>
17 #include <memory>
18 #include <set>
19 
20 #include "TLorentzVector.h"
21 
23 #include "cetlib_except/exception.h"
24 
26 
27 namespace sim
28 {
29 
30  // Static variable for eve ID calculator. We're using an unique_ptr,
31  // so when this object is eventually deleted (at the end of the job)
32  // it will delete the underlying pointer.
33  static std::unique_ptr<EveIdCalculator> eveIdCalculator;
34 
35  //-------------------------------------------------------------
36  // Constructor.
38  {
39  }
40 
41  //-------------------------------------------------------------
42  ParticleNavigator::ParticleNavigator(const std::vector<sim::Particle>& parts)
43  {
44  for(size_t p = 0; p < parts.size(); ++p) Add(new sim::Particle(parts[p]));
45  }
46 
47  //-------------------------------------------------------------
48  // Destructor
50  {
51  this->clear();
52  }
53 
54  //-------------------------------------------------------------
55  // Copy constructor. Note that since this class inherits from
56  // TObject, we have to copy its information explicitly.
58  {
59  // Clear any contents that we already possess.
60  this->clear();
61 
62  // Copy each entry in the other ParticleNavigator.
63  for ( const_iterator entry = rhs.fParticleList.begin();
64  entry != rhs.fParticleList.end(); ++entry ){
65  const sim::Particle* original = (*entry).second;
66  sim::Particle* copy = new sim::Particle( *original );
67  this->insert( copy );
68  }
69  }
70 
71  //-------------------------------------------------------------
72  // Assignment constructor.
74  {
75  // Usual test for self-assignment.
76  if ( this == &rhs ) return *this;
77 
78  // Clear any contents that we already possess.
79  this->clear();
80 
81  // Copy each entry in the other ParticleNavigator.
82  for ( const_iterator entry = rhs.fParticleList.begin();
83  entry != rhs.fParticleList.end(); ++entry ){
84  const sim::Particle* original = (*entry).second;
85  sim::Particle* copy = new sim::Particle( *original );
86  this->insert( copy );
87  }
88 
89  return *this;
90  }
91 
92  //-------------------------------------------------------------
93  // Apply an energy cut to the particles.
94  void ParticleNavigator::Cut( const double& cut )
95  {
96  // The safest way to do this is to create a list of track IDs that
97  // fail the cut, then delete those IDs.
98 
99  // Define a list of IDs.
100  typedef std::set< key_type > keyList_type;
101  keyList_type keyList;
102 
103  // Add each ID that fails the cut to the list.
104  for ( const_iterator i = fParticleList.begin(); i != fParticleList.end(); ++i ){
105  const sim::Particle* particle = (*i).second;
106  Double_t totalInitialEnergy = particle->E();
107  if ( totalInitialEnergy < cut ) {
108  keyList.insert( (*i).first );
109  }
110  }
111 
112  // Go through the list, deleting the particles that are on the list.
113  for ( keyList_type::const_iterator i = keyList.begin(); i != keyList.end(); ++i ){
114  this->erase( *i );
115  }
116  }
117 
118  //-------------------------------------------------
120  {
121  const_iterator i = fParticleList.begin();
122  std::advance(i,index);
123  return (*i).first;
124  }
125 
126  //-------------------------------------------------
128  {
129  const_iterator i = fParticleList.begin();
130  std::advance(i,index);
131  return (*i).second;
132  }
133 
134  //-------------------------------------------------
136  {
137  iterator i = fParticleList.begin();
138  std::advance(i,index);
139  return (*i).second;
140  }
141 
142  //-------------------------------------------------
144  {
145  return fPrimaries.find( trackID ) != fPrimaries.end();
146  }
147 
148  //-------------------------------------------------
150  {
151  return fPrimaries.size();
152  }
153 
154  //-------------------------------------------------
156  {
157  // Advance "index" entries from the beginning of the primary list.
158  primaries_const_iterator primary = fPrimaries.begin();
159  std::advance( primary, index );
160 
161  // Get the track ID from that entry in the list.
162  int trackID = *primary;
163 
164  // Find the entry in the particle list with that track ID.
165  const_iterator entry = fParticleList.find(trackID);
166 
167  // Return the Particle object in that entry.
168  return (*entry).second;
169  }
170 
171  //-------------------------------------------------
173  {
174  // Advance "index" entries from the beginning of the primary list.
175  primaries_const_iterator primary = fPrimaries.begin();
176  std::advance( primary, index );
177 
178  // Get the track ID from that entry in the list.
179  int trackID = *primary;
180 
181  // Find the entry in the particle list with that track ID.
182  iterator entry = fParticleList.find(trackID);
183 
184  // Return the Particle object in that entry.
185  return (*entry).second;
186  }
187 
188  //-------------------------------------------------
190  {
191  // Start with a fresh ParticleNavigator, the destination of the
192  // particles with adjusted track numbers.
194 
195  // For each particle in our list:
196  for ( const_iterator i = fParticleList.begin(); i != fParticleList.end(); ++i ){
197  const sim::Particle* particle = (*i).second;
198 
199  // Create a new particle with an adjusted track ID.
200  sim::Particle* adjusted = new sim::Particle( particle->TrackId() + offset,
201  particle->PdgCode(),
202  particle->Process(),
203  particle->Mother(),
204  particle->Mass() );
205 
206  adjusted->SetPolarization( particle->Polarization() );
207 
208  // Copy all the daughters, adjusting the track ID.
209  for ( int d = 0; d < particle->NumberDaughters(); ++d ){
210  int daughterID = particle->Daughter(d);
211  adjusted->AddDaughter( daughterID + offset );
212  }
213 
214  // Copy the trajectory points.
215  for ( size_t t = 0; t < particle->NumberTrajectoryPoints(); ++t ){
216  adjusted->AddTrajectoryPoint( particle->Position(t), particle->Momentum(t) );
217  }
218 
219  // Add the adjusted particle to the destination particle list.
220  // This will also adjust the destination's list of primary
221  // particles, if needed.
222  result.insert( adjusted );
223  }
224 
225  return result;
226  }
227 
228  //-------------------------------------------------------------
229  // Just in case: define the result of "scalar * ParticleNavigator" to be
230  // the same as "ParticleNavigator * scalar".
232  {
233  return list + value;
234  }
235 
236 
237  // This is the main "insertion" method for the ParticleNavigator
238  // pseudo-array pseudo-map. It does the following:
239  // - Add the Particle to the list; if the track ID is already in the
240  // list, throw an exception.
241  // - If it's a primary particle, add it to the list of primaries.
243  {
244  int trackID = particle->TrackId();
245  iterator insertion = fParticleList.lower_bound( trackID );
246  if ( insertion == fParticleList.end() ){
247  // The best "hint" we can give is that the particle will go at
248  // the end of the list.
249  fParticleList.insert( insertion, value_type( trackID, particle ) );
250  }
251  else if ( (*insertion).first == trackID ){
252  throw cet::exception("ParticleNavigator") << "sim::ParticleNavigator::insert - ERROR - "
253  << "track ID=" << trackID
254  << " is already in the list";
255  }
256  else{
257  // It turns out that the best hint we can give is one more
258  // than the result of lower_bound.
259  fParticleList.insert( ++insertion, value_type( trackID, particle ) );
260  }
261 
262  // If this is a primary particle, add it to the list. Look to see
263  // if the process name contains the string rimary - leave the
264  // "p" off cause it might be capitalized, but then again maybe not
265  if ( particle->Process().find("rimary") != std::string::npos )
266  fPrimaries.insert( trackID );
267  }
268 
269  //-------------------------------------------------------------
271  {
272  for ( iterator i = fParticleList.begin(); i != fParticleList.end(); ++i ){
273  delete (*i).second;
274  }
275 
276  fParticleList.clear();
277  fPrimaries.clear();
278 
279  eveIdCalculator.reset();
280 
281  return;
282  }
283 
284  //-------------------------------------------------------------
285  // An erase that includes the deletion of the associated Particle*.
287  {
288  iterator entry = fParticleList.find( key );
289  delete (*entry).second;
290  return fParticleList.erase( key );
291  }
292 
293  //-------------------------------------------------------------
294  std::ostream& operator<< ( std::ostream& output, const ParticleNavigator& list )
295  {
296  // Determine a field width for the particle number.
297  ParticleNavigator::size_type numberOfParticles = list.size();
298  int numberOfDigits = (int) std::log10( (double) numberOfParticles ) + 1;
299 
300  // A simple header.
301  output.width( numberOfDigits );
302  output << "#" << ": < ID, particle >" << std::endl;
303 
304  // Write each particle on a separate line.
305  ParticleNavigator::size_type nParticle = 0;
306  for ( ParticleNavigator::const_iterator particle = list.begin();
307  particle != list.end(); ++particle, ++nParticle ){
308  output.width( numberOfDigits );
309  output << nParticle << ": "
310  << "<" << (*particle).first
311  << "," << *((*particle).second)
312  << ">" << std::endl;
313  }
314 
315  return output;
316  }
317 
318  //-------------------------------------------------------------
319  // operator[] in a non-const context.
321  {
322  const_iterator i = fParticleList.find(key);
323  if(i == fParticleList.end())
324  throw cet::exception("ParticleNavigator") << "track id is not in map";
325  return (*i).second;
326  }
327 
328  //-------------------------------------------------------------
329  // operator[] in a const context; not usual for maps, but I find it useful.
331  {
332  const_iterator i = fParticleList.find(key);
333  if(i == fParticleList.end())
334  throw cet::exception("ParticleNavigator") << "track id is not in map";
335  return (*i).second;
336  }
337 
338 
339  //-------------------------------------------------------------
340  // The eve ID calculation.
341  int ParticleNavigator::EveId( const int trackID ) const
342  {
343  // If the eve ID calculator has never been initialized, use the
344  // default method.
345  if ( eveIdCalculator.get() == 0 ){
347  }
348 
349  // If the eve ID calculator has changed, or we're looking at a
350  // different ParticleNavigator, initialize the calculator.
351  static EveIdCalculator* saveEveIdCalculator = 0;
352  if ( saveEveIdCalculator != eveIdCalculator.get() ) {
353  saveEveIdCalculator = eveIdCalculator.get();
354  eveIdCalculator->Init( this );
355  }
356  if ( eveIdCalculator->ParticleNavigator() != this ){
357  eveIdCalculator->Init( this );
358  }
359 
360  // After the "bookkeeping" tests, here's where we actually do the
361  // calculation.
362  return eveIdCalculator->CalculateEveId( trackID );
363  }
364 
365  //-------------------------------------------------------------
366  // Save a new eve ID calculation method.
368  {
369  eveIdCalculator.reset(calc);
370  }
371 
372 } // end namespace sim
373 ////////////////////////////////////////////////////////////////////////
double E(const int i=0) const
Definition: MCParticle.h:232
list_type::size_type size_type
list_type::key_type key_type
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
ofstream output
const TVector3 & Polarization() const
Definition: MCParticle.h:213
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
list_type fParticleList
Sorted list of particles in the event.
void AddDaughter(const int trackID)
Definition: MCParticle.h:264
int PdgCode() const
Definition: MCParticle.h:211
friend std::ostream & operator<<(std::ostream &output, const ParticleNavigator &)
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:256
static void AdoptEveIdCalculator(EveIdCalculator *)
size_type size() const
int Mother() const
Definition: MCParticle.h:212
const char * p
Definition: xmltok.h:285
double Mass() const
Definition: MCParticle.h:238
list_type::const_iterator const_iterator
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void Cut(const double &)
std::string Process() const
Definition: MCParticle.h:214
osc::OscCalcDumb calc
list_type::iterator iterator
int NumberDaughters() const
Definition: MCParticle.h:216
size_type erase(const key_type &key)
list_type::mapped_type mapped_type
int TrackId() const
Definition: MCParticle.h:209
int Daughter(const int i) const
Definition: MCParticle.cxx:112
const key_type & TrackId(const size_type) const
void SetPolarization(const TVector3 &p)
Definition: MCParticle.h:265
ParticleNavigator operator+(const int &value) const
void Init(const sim::ParticleNavigator *list)
Initialize this calculator for a particular ParticleNavigator.
const sim::Particle * Primary(const int) const
const XML_Char int const XML_Char * value
Definition: expat.h:331
Float_t d
Definition: plot.C:236
primaries_type::const_iterator primaries_const_iterator
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
Definition: FillTruth.h:15
bool IsPrimary(int trackID) const
const Cut cut
Definition: exporter_fd.C:30
static std::unique_ptr< EveIdCalculator > eveIdCalculator
T log10(T number)
Definition: d0nt_math.hpp:120
void insert(sim::Particle *value)
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
mapped_type operator[](const key_type &key) const
ParticleNavigator Add(const int &offset) const
Int_t trackID
Definition: plot.C:84
ParticleNavigator & operator=(const ParticleNavigator &rhs)
def entry(str)
Definition: HTMLTools.py:26
list_type::value_type value_type
primaries_type fPrimaries
Sorted list of the track IDs of primary particles.
mapped_type Particle(const size_type) const
int EveId(const int trackID) const