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