MCParticle.h
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file MCParticle.h
3 /// \brief Particle class
4 /// \version $Id: MCParticle.h,v 1.16 2012-11-20 17:39:38 brebel Exp $
5 /// \author brebel@fnal.gov
6 ////////////////////////////////////////////////////////////////////////
7 
8 /// This class describes a particle created in the detector Monte
9 /// Carlo simulation.
10 
11 #ifndef SIMB_MCPARTICLE_H
12 #define SIMB_MCPARTICLE_H
13 
15 
16 #include <set>
17 #include <string>
18 #include <iostream>
19 #include "TVector3.h"
20 #include "TLorentzVector.h"
21 
22 namespace simb {
23 
24  class MCParticle {
25  public:
26 
27  // An indicator for an uninitialized variable (see MCParticle.cxx).
28  static const int s_uninitialized; //! Don't write this as ROOT output
29 
30  MCParticle();
31 
32  protected:
33  typedef std::set<int> daughters_type;
34 
35  int fstatus; ///< Status code from generator, geant, etc
36  int ftrackId; ///< TrackId
37  int fpdgCode; ///< PDG code
38  int fmother; ///< Mother
39  std::string fprocess; ///< Detector-simulation physics process that created the particle
40  std::string fendprocess; ///< end process for the particle
41  simb::MCTrajectory ftrajectory; ///< particle trajectory (position,momentum)
42  double fmass; ///< Mass; from PDG unless overridden Should be in GeV
43  TVector3 fpolarization; ///< Polarization
44  daughters_type fdaughters; ///< Sorted list of daughters of this particle.
45  double fWeight; ///< Assigned weight to this particle for MC tests
46  TLorentzVector fGvtx; ///< Vertex needed by generater (genie) to rebuild
47  ///< genie::EventRecord for event reweighting
48  int frescatter; ///< rescatter code
49 
50  public:
51 
52  // Standard constructor. If the mass is not supplied in the
53  // argument, then the PDG mass is used.
54  // status code = 1 means the particle is to be tracked, default it to be tracked
55  // mother = -1 means that this particle has no mother
56  MCParticle(const int trackId,
57  const int pdg,
58  const std::string process,
59  const int mother = -1,
60  const double mass = s_uninitialized,
61  const int status = 1);
62 
63 
64  // our own copy and move assignment constructors (default)
65  MCParticle(MCParticle const &) = default; // Copy constructor.
66  MCParticle& operator=( const MCParticle&) = default;
67  MCParticle(MCParticle&&) = default;
68  MCParticle& operator= (MCParticle&&) = default;
69 
70 
71  // constructor for copy from MCParticle, but with offset trackID
72  MCParticle(MCParticle const&, int);
73 
74  // Accessors.
75  //
76  // The track ID number assigned by the Monte Carlo. This will be
77  // unique for each Particle in an event. - 0 for primary particles
78  int TrackId() const;
79 
80  // Get at the status code returned by GENIE, Geant4, etc
81  int StatusCode() const;
82 
83  // The PDG code of the particle. Note that Geant4 uses the
84  // "extended" system for encoding nuclei; e.g., 1000180400 is an
85  // Argon nucleus. See "Monte Carlo PArticle Numbering Scheme" in
86  // any Review of Particle Physics.
87  int PdgCode() const;
88 
89  // The track ID of the mother particle. Note that it's possible
90  // for a particle to have a mother that's not recorded in the
91  // Particle list; e.g., an excited nucleus with low kinetic energy
92  // emits a photon with high kinetic energy.
93  int Mother() const;
94 
95  const TVector3& Polarization() const;
96  void SetPolarization( const TVector3& p );
97 
98  // The detector-simulation physics process that created the
99  // particle. If this is a primary particle, it will have the
100  // value "primary"
101  std::string Process() const;
102 
103  std::string EndProcess() const;
105 
106  // Accessors for daughter information. Note that it's possible
107  // (even likely) for a daughter track not to be included in a
108  // Particle list, if that daughter particle falls below the energy cut.
109  void AddDaughter( const int trackID );
110  int NumberDaughters() const;
111  int Daughter(const int i) const; //> Returns the track ID for the "i-th" daughter.
112 
113  // Accessors for trajectory information.
114  unsigned int NumberTrajectoryPoints() const;
115 
116  // To avoid confusion with the X() and Y() methods of MCTruth
117  // (which return Feynmann x and y), use "Vx,Vy,Vz" for the
118  // vertex.
119  const TLorentzVector& Position( const int i = 0 ) const;
120  double Vx(const int i = 0) const;
121  double Vy(const int i = 0) const;
122  double Vz(const int i = 0) const;
123  double T(const int i = 0) const;
124 
125  const TLorentzVector& EndPosition() const;
126  double EndX() const;
127  double EndY() const;
128  double EndZ() const;
129  double EndT() const;
130 
131  const TLorentzVector& Momentum( const int i = 0 ) const;
132  double Px(const int i = 0) const;
133  double Py(const int i = 0) const;
134  double Pz(const int i = 0) const;
135  double E(const int i = 0) const;
136  double P(const int i = 0) const;
137  double Pt(const int i = 0) const;
138  double Mass() const;
139 
140  const TLorentzVector& EndMomentum() const;
141  double EndPx() const;
142  double EndPy() const;
143  double EndPz() const;
144  double EndE() const;
145 
146  // Getters and setters for the generator vertex
147  // These are for setting the generator vertex. In the case of genie
148  // the generator assumes a cooridnate system with origin at the nucleus.
149  // These variables save the particle vertexs in this cooridnate system.
150  // After genie generates the event, a cooridnate transformation is done
151  // to place the event in the detector cooridnate system. These variables
152  // store the vertex before that cooridnate transformation happens.
153  void SetGvtx(double *v);
154  void SetGvtx(float *v);
155  void SetGvtx(TLorentzVector v);
156  void SetGvtx(double x,
157  double y,
158  double z,
159  double t);
160  TLorentzVector GetGvtx() const;
161  double Gvx() const;
162  double Gvy() const;
163  double Gvz() const;
164  double Gvt() const;
165 
166  //Getters and setters for first and last daughter data members
167  int FirstDaughter() const;
168  int LastDaughter() const;
169 
170  //Getters and setters for rescatter status
171  void SetRescatter(int code);
172  int Rescatter() const;
173 
174  // Access to the trajectory in both a const and non-const context.
175  const simb::MCTrajectory& Trajectory() const;
176 
177  // Make it easier to add a (position,momentum) point to the
178  // trajectory. You must add this information for every point you wish to keep
179  void AddTrajectoryPoint(TLorentzVector const& position,
180  TLorentzVector const& momentum );
181  void AddTrajectoryPoint(TLorentzVector const& position,
182  TLorentzVector const& momentum,
183  std::string const& process);
184 
185  // methods for giving/accessing a weight to this particle for use
186  // in studies of rare processes, etc
187  double Weight() const;
188  void SetWeight(double wt);
189 
190  void SparsifyTrajectory();
191 
192  // Define a comparison operator for particles. This allows us to
193  // keep them in sets or maps. It makes sense to order a list of
194  // particles by track ID... but take care! After we get past the
195  // primary particles in an event, it is NOT safe to assume that a
196  // particle with a lower track ID is "closer" to the event
197  // vertex.
198  bool operator<( const simb::MCParticle& other ) const;
199 
200  friend std::ostream& operator<< ( std::ostream& output, const simb::MCParticle& );
201  };
202 
203 } // namespace simb
204 
205 #include <functional> // so we can redefine less<> below
206 #include <math.h>
207 
208 // methods to access data members and other information
209 inline int simb::MCParticle::TrackId() const { return ftrackId; }
210 inline int simb::MCParticle::StatusCode() const { return fstatus; }
211 inline int simb::MCParticle::PdgCode() const { return fpdgCode; }
212 inline int simb::MCParticle::Mother() const { return fmother; }
213 inline const TVector3& simb::MCParticle::Polarization() const { return fpolarization; }
216 inline int simb::MCParticle::NumberDaughters() const { return fdaughters.size(); }
217 inline unsigned int simb::MCParticle::NumberTrajectoryPoints() const { return ftrajectory.size(); }
218 inline const TLorentzVector& simb::MCParticle::Position( const int i ) const { return ftrajectory.Position(i); }
219 inline const TLorentzVector& simb::MCParticle::Momentum( const int i ) const { return ftrajectory.Momentum(i); }
220 inline double simb::MCParticle::Vx(const int i) const { return Position(i).X(); }
221 inline double simb::MCParticle::Vy(const int i) const { return Position(i).Y(); }
222 inline double simb::MCParticle::Vz(const int i) const { return Position(i).Z(); }
223 inline double simb::MCParticle::T(const int i) const { return Position(i).T(); }
224 inline const TLorentzVector& simb::MCParticle::EndPosition() const { return Position(ftrajectory.size()-1); }
225 inline double simb::MCParticle::EndX() const { return Position(ftrajectory.size()-1).X(); }
226 inline double simb::MCParticle::EndY() const { return Position(ftrajectory.size()-1).Y(); }
227 inline double simb::MCParticle::EndZ() const { return Position(ftrajectory.size()-1).Z(); }
228 inline double simb::MCParticle::EndT() const { return Position(ftrajectory.size()-1).T(); }
229 inline double simb::MCParticle::Px(const int i) const { return Momentum(i).Px(); }
230 inline double simb::MCParticle::Py(const int i) const { return Momentum(i).Py(); }
231 inline double simb::MCParticle::Pz(const int i) const { return Momentum(i).Pz(); }
232 inline double simb::MCParticle::E(const int i) const { return Momentum(i).E(); }
233 inline double simb::MCParticle::P(const int i) const { return std::sqrt(std::pow(Momentum(i).E(),2.)
234  - std::pow(fmass,2.)); }
235 inline double simb::MCParticle::Pt(const int i) const { return std::sqrt( std::pow(Momentum(i).Px(),2.)
236  + std::pow(Momentum(i).Py(),2.)); }
237 
238 inline double simb::MCParticle::Mass() const { return fmass; }
239 inline const TLorentzVector& simb::MCParticle::EndMomentum() const { return Momentum(ftrajectory.size()-1); }
240 inline double simb::MCParticle::EndPx() const { return Momentum(ftrajectory.size()-1).X(); }
241 inline double simb::MCParticle::EndPy() const { return Momentum(ftrajectory.size()-1).Y(); }
242 inline double simb::MCParticle::EndPz() const { return Momentum(ftrajectory.size()-1).Z(); }
243 inline double simb::MCParticle::EndE() const { return Momentum(ftrajectory.size()-1).T(); }
244 inline TLorentzVector simb::MCParticle::GetGvtx() const { return fGvtx; }
245 inline double simb::MCParticle::Gvx() const { return fGvtx.X(); }
246 inline double simb::MCParticle::Gvy() const { return fGvtx.Y(); }
247 inline double simb::MCParticle::Gvz() const { return fGvtx.Z(); }
248 inline double simb::MCParticle::Gvt() const { return fGvtx.T(); }
249 inline int simb::MCParticle::FirstDaughter() const { return *(fdaughters.begin()); }
250 inline int simb::MCParticle::LastDaughter() const { return *(fdaughters.rbegin()); }
251 inline int simb::MCParticle::Rescatter() const { return frescatter; }
253 inline double simb::MCParticle::Weight() const { return fWeight; }
254 
255 // methods to set information
256 inline void simb::MCParticle::AddTrajectoryPoint(TLorentzVector const& position,
257  TLorentzVector const& momentum )
258  { ftrajectory.Add( position, momentum ); }
259 inline void simb::MCParticle::AddTrajectoryPoint(TLorentzVector const& position,
260  TLorentzVector const& momentum,
261  std::string const& process)
262  { ftrajectory.Add( position, momentum, process); }
264 inline void simb::MCParticle::AddDaughter(int const trackID) { fdaughters.insert(trackID); }
265 inline void simb::MCParticle::SetPolarization(TVector3 const& p) { fpolarization = p; }
266 inline void simb::MCParticle::SetRescatter(int code) { frescatter = code; }
267 inline void simb::MCParticle::SetWeight(double wt) { fWeight = wt; }
268 
269 // definition of the < operator
270 inline bool simb::MCParticle::operator<( const simb::MCParticle& other ) const { return ftrackId < other.ftrackId; }
271 
272 // A potentially handy definition: At this stage, I'm not sure
273 // whether I'm going to be keeping a list based on Particle or on
274 // Particle*. We've already defined operator<(Particle,Particle),
275 // that is, how to compare two Particle objects; by default that also
276 // defines less<Particle>, which is what the STL containers use for
277 // comparisons.
278 
279 // The following defines less<Particle*>, that is, how to compare two
280 // Particle*: by looking at the objects, not at the pointer
281 // addresses. The result is that, e.g., a set<Particle*> will be
282 // sorted in the order I expect.
283 
284 namespace std {
285  template <>
286  class less<simb::MCParticle*>
287  {
288  public:
289  bool operator()( const simb::MCParticle* lhs, const simb::MCParticle* rhs )
290  {
291  return (*lhs) < (*rhs);
292  }
293  };
294 } // std
295 
296 #endif // SIMB_MCPARTICLE_H
double E(const int i=0) const
Definition: MCParticle.h:232
void Add(TLorentzVector const &p, TLorentzVector const &m)
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
ofstream output
const TVector3 & Polarization() const
Definition: MCParticle.h:213
void SparsifyTrajectory()
Definition: MCParticle.h:263
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
void AddDaughter(const int trackID)
Definition: MCParticle.h:264
int PdgCode() const
Definition: MCParticle.h:211
double Py(const int i=0) const
Definition: MCParticle.h:230
double Gvz() const
Definition: MCParticle.h:247
int status
Definition: fabricate.py:1613
Trajectory class.
friend std::ostream & operator<<(std::ostream &output, const simb::MCParticle &)
Definition: MCParticle.cxx:151
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:256
double EndE() const
Definition: MCParticle.h:243
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:224
double EndZ() const
Definition: MCParticle.h:227
double Weight() const
Definition: MCParticle.h:253
static const int s_uninitialized
Definition: MCParticle.h:28
int FirstDaughter() const
Definition: MCParticle.h:249
TLorentzVector fGvtx
Definition: MCParticle.h:46
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:252
int Mother() const
Definition: MCParticle.h:212
int fmother
Mother.
Definition: MCParticle.h:38
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
double Mass() const
Definition: MCParticle.h:238
constexpr T pow(T x)
Definition: pow.h:75
double fmass
Mass; from PDG unless overridden Should be in GeV.
Definition: MCParticle.h:42
double Px(const int i=0) const
Definition: MCParticle.h:229
std::string fendprocess
end process for the particle
Definition: MCParticle.h:40
int StatusCode() const
Definition: MCParticle.h:210
double Gvx() const
Definition: MCParticle.h:245
double Gvy() const
Definition: MCParticle.h:246
std::string Process() const
Definition: MCParticle.h:214
TLorentzVector GetGvtx() const
Definition: MCParticle.h:244
double EndY() const
Definition: MCParticle.h:226
int NumberDaughters() const
Definition: MCParticle.h:216
int ftrackId
TrackId.
Definition: MCParticle.h:36
int TrackId() const
Definition: MCParticle.h:209
int Daughter(const int i) const
Definition: MCParticle.cxx:112
Float_t Y
Definition: plot.C:38
Float_t Z
Definition: plot.C:38
void Sparsify(double margin=.1)
bool operator()(const simb::MCParticle *lhs, const simb::MCParticle *rhs)
Definition: MCParticle.h:289
const XML_Char * s
Definition: expat.h:262
double Pt(const int i=0) const
Definition: MCParticle.h:235
void SetPolarization(const TVector3 &p)
Definition: MCParticle.h:265
int frescatter
rescatter code
Definition: MCParticle.h:48
std::string EndProcess() const
Definition: MCParticle.h:215
int fpdgCode
PDG code.
Definition: MCParticle.h:37
MCParticle()
Don&#39;t write this as ROOT output.
Definition: MCParticle.cxx:32
double EndPz() const
Definition: MCParticle.h:242
double P(const int i=0) const
Definition: MCParticle.h:233
daughters_type fdaughters
Sorted list of daughters of this particle.
Definition: MCParticle.h:44
double T(const int i=0) const
Definition: MCParticle.h:223
bool operator<(const simb::MCParticle &other) const
Definition: MCParticle.h:270
std::string fprocess
Detector-simulation physics process that created the particle.
Definition: MCParticle.h:39
double EndT() const
Definition: MCParticle.h:228
const TLorentzVector & Position(const size_type) const
The accessor methods described above.
std::set< int > daughters_type
Definition: MCParticle.h:33
z
Definition: test.py:28
This class describes a particle created in the detector Monte Carlo simulation.
void SetWeight(double wt)
Definition: MCParticle.h:267
void SetEndProcess(std::string s)
Definition: MCParticle.cxx:105
simb::MCTrajectory ftrajectory
particle trajectory (position,momentum)
Definition: MCParticle.h:41
Definition: inftrees.h:24
double Vx(const int i=0) const
Definition: MCParticle.h:220
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetGvtx(double *v)
Definition: MCParticle.cxx:120
size_type size() const
Definition: MCTrajectory.h:165
double EndPy() const
Definition: MCParticle.h:241
double fWeight
Assigned weight to this particle for MC tests.
Definition: MCParticle.h:45
int LastDaughter() const
Definition: MCParticle.h:250
double Gvt() const
Definition: MCParticle.h:248
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
double Pz(const int i=0) const
Definition: MCParticle.h:231
int Rescatter() const
Definition: MCParticle.h:251
double Vz(const int i=0) const
Definition: MCParticle.h:222
int fstatus
Status code from generator, geant, etc.
Definition: MCParticle.h:35
double EndPx() const
Definition: MCParticle.h:240
TVector3 fpolarization
Polarization.
Definition: MCParticle.h:43
Int_t trackID
Definition: plot.C:84
double EndX() const
Definition: MCParticle.h:225
void SetRescatter(int code)
Definition: MCParticle.h:266
Float_t X
Definition: plot.C:38
const TLorentzVector & Momentum(const size_type) const
const TLorentzVector & EndMomentum() const
Definition: MCParticle.h:239
MCParticle & operator=(const MCParticle &)=default
double Vy(const int i=0) const
Definition: MCParticle.h:221