GHepParticle.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Sep 15, 2009 - CA
14  Remove redundant IsFake(), IsNucleus() and IsParticle() methods. Removed
15  the corresponding priv data members. Use pdg::IsPseudoParticle(),
16  pdg::IsIon() and pdg::IsParticle instead.
17  @ Sep 15, 2009 - CA
18  Added 'rescattering code' to allow intranuclear hadron transport MCs to
19  store a hadronic reaction code which can not always be easily recreated
20  from the particle list.
21  @ May 05, 2010 - CR
22  Adding special ctor for ROOT I/O purposes so as to avoid memory leak due to
23  memory allocated in the default ctor when objects of this class are read by
24  the ROOT Streamer.
25 
26 */
27 //____________________________________________________________________________
28 
29 #include <cstdlib>
30 #include <cassert>
31 #include <iomanip>
32 
33 #include <TMath.h>
34 #include <TRootIOCtor.h>
35 
44 
45 using namespace genie;
46 using namespace genie::constants;
47 
48 using std::setw;
49 using std::setprecision;
50 using std::setfill;
51 using std::ios;
52 using std::cout;
53 
54 const double kPCutOff = 1e-15;
55 const double kOffShellDm = 0.002; // 2 MeV
56 
58 
59 //____________________________________________________________________________
61  ostream & operator<< (ostream& stream, const GHepParticle & particle)
62  {
63  particle.Print(stream);
64  return stream;
65  }
66 }
67 //___________________________________________________________________________
69 TObject()
70 {
71  this->Init();
72 }
73 //___________________________________________________________________________
74 // TParticle-like constructor
76  int mother1, int mother2, int daughter1, int daughter2,
77  const TLorentzVector & p, const TLorentzVector & v) :
78 TObject(),
79 fStatus(status),
80 fFirstMother(mother1),
81 fLastMother(mother2),
82 fFirstDaughter(daughter1),
83 fLastDaughter(daughter2)
84 {
85  this->SetPdgCode(pdg);
86 
87  fP4 = new TLorentzVector(p);
88  fX4 = new TLorentzVector(v);
89 
90  fRescatterCode = -1;
91  fPolzTheta = -999;
92  fPolzPhi = -999;
93  fIsBound = false;
94  fRemovalEnergy = 0.;
95 }
96 //___________________________________________________________________________
97 // TParticle-like constructor
99  int mother1, int mother2, int daughter1, int daughter2,
100  double px, double py, double pz, double En,
101  double x, double y, double z, double t) :
102 TObject(),
103 fStatus(status),
104 fFirstMother(mother1),
105 fLastMother(mother2),
106 fFirstDaughter(daughter1),
107 fLastDaughter(daughter2)
108 {
109  this->SetPdgCode(pdg);
110 
111  fP4 = new TLorentzVector(px,py,pz,En);
112  fX4 = new TLorentzVector(x,y,z,t);
113 
114  fRescatterCode = -1;
115  fPolzTheta = -999;
116  fPolzPhi = -999;
117  fIsBound = false;
118  fRemovalEnergy = 0.;
119 }
120 //___________________________________________________________________________
121 // Copy constructor
123 TObject()
124 {
125  this->Init();
126  this->Copy(particle);
127 }
128 //___________________________________________________________________________
130 TObject(),
131 fPdgCode(0),
133 fRescatterCode(-1),
134 fFirstMother(-1),
135 fLastMother(-1),
136 fFirstDaughter(-1),
137 fLastDaughter(-1),
138 fP4(0),
139 fX4(0),
140 fPolzTheta(-999.),
141 fPolzPhi(-999.),
142 fRemovalEnergy(0),
143 fIsBound(false)
144 {
145 
146 }
147 //___________________________________________________________________________
149 {
150  this->CleanUp();
151 }
152 //___________________________________________________________________________
153 string GHepParticle::Name(void) const
154 {
155  this->AssertIsKnownParticle();
156 
157  TParticlePDG * p = PDGLibrary::Instance()->Find(fPdgCode);
158  return p->GetName();
159 }
160 //___________________________________________________________________________
161 double GHepParticle::Mass(void) const
162 {
163  this->AssertIsKnownParticle();
164 
165  TParticlePDG * p = PDGLibrary::Instance()->Find(fPdgCode);
166  return p->Mass();
167 }
168 //___________________________________________________________________________
169 double GHepParticle::Charge(void) const
170 {
171  this->AssertIsKnownParticle();
172 
173  TParticlePDG * p = PDGLibrary::Instance()->Find(fPdgCode);
174  return p->Charge();
175 }
176 //___________________________________________________________________________
177 double GHepParticle::KinE(bool mass_from_pdg) const
178 {
179  if(!fP4) {
180  LOG("GHepParticle", pWARN) << "4-momentum not yet set!";
181  return 0;
182  }
183 
184  double En = fP4->Energy();
185  double M = ( (mass_from_pdg) ? this->Mass() : fP4->M() );
186  double K = En - M;
187 
188  K = TMath::Max(K,0.);
189  return K;
190 }
191 //___________________________________________________________________________
192 int GHepParticle::Z(void) const
193 {
194 // Decoding Z from the PDG code
195 
196  if(!pdg::IsIon(fPdgCode))
197  return -1;
198  else
200 }
201 //___________________________________________________________________________
202 int GHepParticle::A(void) const
203 {
204 // Decoding A from the PDG code
205 
206  if(!pdg::IsIon(fPdgCode))
207  return -1;
208  else
210 }
211 //___________________________________________________________________________
212 TLorentzVector * GHepParticle::GetP4(void) const
213 {
214 // see GHepParticle::P4() for a method that does not create a new object and
215 // transfers its ownership
216 
217  if(fP4) {
218  TLorentzVector * p4 = new TLorentzVector(*fP4);
219 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
220  LOG("GHepParticle", pDEBUG)
221  << "Return vp = " << utils::print::P4AsShortString(p4);
222 #endif
223  return p4;
224  } else {
225  LOG("GHepParticle", pWARN) << "NULL 4-momentum TLorentzVector";
226  return 0;
227  }
228 }
229 //___________________________________________________________________________
230 TLorentzVector * GHepParticle::GetX4(void) const
231 {
232 // see GHepParticle::X4() for a method that does not create a new object and
233 // transfers its ownership
234 
235  if(fX4) {
236  TLorentzVector * x4 = new TLorentzVector(*fX4);
237 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
238  LOG("GHepParticle", pDEBUG)
239  << "Return x4 = " << utils::print::X4AsString(x4);
240 #endif
241  return x4;
242  } else {
243  LOG("GHepParticle", pWARN) << "NULL 4-position TLorentzVector";
244  return 0;
245  }
246 }
247 //___________________________________________________________________________
249 {
250  fPdgCode = code;
251  this->AssertIsKnownParticle();
252 }
253 //___________________________________________________________________________
254 void GHepParticle::SetMomentum(const TLorentzVector & p4)
255 {
256  if(fP4)
257  fP4->SetPxPyPzE( p4.Px(), p4.Py(), p4.Pz(), p4.Energy() );
258  else
259  fP4 = new TLorentzVector(p4);
260 }
261 //___________________________________________________________________________
262 void GHepParticle::SetMomentum(double px, double py, double pz, double En)
263 {
264  if(fP4)
265  fP4->SetPxPyPzE(px, py, pz, En);
266  else
267  fP4 = new TLorentzVector(px, py, pz, En);
268 }
269 //___________________________________________________________________________
270 void GHepParticle::SetPosition(const TLorentzVector & v4)
271 {
272  this->SetPosition(v4.X(), v4.Y(), v4.Z(), v4.T());
273 }
274 //___________________________________________________________________________
275 void GHepParticle::SetPosition(double x, double y, double z, double t)
276 {
277 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
278  LOG("GHepParticle", pDEBUG)
279  << "Setting position to (x = " << x << ", y = "
280  << y << ", z = " << z << ", t = " << t << ")";
281 #endif
282 
283  if(fX4) fX4->SetXYZT(x,y,z,t);
284  else fX4 = new TLorentzVector(x,y,z,t);
285 }
286 //___________________________________________________________________________
287 void GHepParticle::SetEnergy(double En)
288 {
289  this->SetMomentum(this->Px(), this->Py(), this->Pz(), En);
290 }
291 //___________________________________________________________________________
292 void GHepParticle::SetPx(double px)
293 {
294  this->SetMomentum(px, this->Py(), this->Pz(), this->E());
295 }
296 //___________________________________________________________________________
297 void GHepParticle::SetPy(double py)
298 {
299  this->SetMomentum(this->Px(), py, this->Pz(), this->E());
300 }
301 //___________________________________________________________________________
302 void GHepParticle::SetPz(double pz)
303 {
304  this->SetMomentum(this->Px(), this->Py(), pz, this->E());
305 }
306 //___________________________________________________________________________
308 {
309  this->AssertIsKnownParticle();
310 
311  TParticlePDG * p = PDGLibrary::Instance()->Find(fPdgCode);
312 
313  double Mpdg = p->Mass();
314  double M4p = (fP4) ? fP4->M() : 0.;
315 
316 // return utils::math::AreEqual(Mpdg, M4p);
317 
318  return (TMath::Abs(M4p-Mpdg) < kOffShellDm);
319 }
320 //___________________________________________________________________________
322 {
323  return (! this->IsOnMassShell());
324 }
325 //___________________________________________________________________________
326 bool GHepParticle::PolzIsSet(void) const
327 {
328 // checks whether the polarization angles have been set
329 
330  return (fPolzTheta > -999 && fPolzPhi > -999);
331 }
332 //___________________________________________________________________________
333 void GHepParticle::GetPolarization(TVector3 & polz)
334 {
335 // gets the polarization vector
336 
337  if(! this->PolzIsSet() ) {
338  polz.SetXYZ(0.,0.,0.);
339  return;
340  }
341  polz.SetX( TMath::Sin(fPolzTheta) * TMath::Cos(fPolzPhi) );
342  polz.SetY( TMath::Sin(fPolzTheta) * TMath::Sin(fPolzPhi) );
343  polz.SetZ( TMath::Cos(fPolzTheta) );
344 }
345 //___________________________________________________________________________
346 void GHepParticle::SetPolarization(double theta, double phi)
347 {
348 // sets the polarization angles
349 
350  if(theta>=0 && theta<=kPi && phi>=0 && phi<2*kPi)
351  {
352  fPolzTheta = theta;
353  fPolzPhi = phi;
354 
355  } else {
356  LOG("GHepParticle", pERROR)
357  << "Invalid polarization angles (polar = " << theta
358  << ", azimuthal = " << phi << ")";
359  }
360 }
361 //___________________________________________________________________________
362 void GHepParticle::SetPolarization(const TVector3 & polz)
363 {
364 // sets the polarization angles
365 
366  double p = polz.Mag();
367  if(! (p>0) ) {
368  LOG("GHepParticle", pERROR)
369  << "Input polarization vector has non-positive norm! Ignoring it";
370  return;
371  }
372 
373  double theta = TMath::ACos(polz.z()/p);
374  double phi = kPi + TMath::ATan2(-polz.y(), -polz.x());
375 
376  this->SetPolarization(theta,phi);
377 }
378 //___________________________________________________________________________
379 void GHepParticle::SetBound(bool bound)
380 {
381  // only set it for p or n
382  bool is_nucleon = pdg::IsNeutronOrProton(fPdgCode);
383  if(!is_nucleon && bound) {
384  LOG("GHepParticle", pERROR)
385  << "Refusing to set the bound flag for particles other than nucleons";
386  LOG("GHepParticle", pERROR)
387  << "(Requested for pdg = " << fPdgCode << ")";
388  return;
389  }
390  // if the particles isn't bound then make sure that its removal energy = 0
391  if(!bound) {
392  fRemovalEnergy = 0;
393  }
394  // set the flag
395  fIsBound = bound;
396 }
397 //___________________________________________________________________________
399 {
400  fRemovalEnergy = TMath::Max(Erm, 0.); // non-negative
401 
402  // if a value was set, make sure that the IsBound flag is turned on
403  if(fRemovalEnergy>0) this->SetBound(true);
404 }
405 //___________________________________________________________________________
407 {
408  fPdgCode = 0;
410  fRescatterCode = -1;
411  fFirstMother = -1;
412  fLastMother = -1;
413  fFirstDaughter = -1;
414  fLastDaughter = -1;
415  fPolzTheta = -999;
416  fPolzPhi = -999;
417  fIsBound = false;
418  fRemovalEnergy = 0.;
419  fP4 = new TLorentzVector(0,0,0,0);
420  fX4 = new TLorentzVector(0,0,0,0);
421 }
422 //___________________________________________________________________________
424 {
425 // deallocate memory
426 
427  if(fP4) delete fP4;
428  if(fX4) delete fX4;
429  fP4 = 0;
430  fX4 = 0;
431 }
432 //___________________________________________________________________________
434 {
435 // deallocate memory + initialize
436 
437  this->CleanUp();
438  this->Init();
439 }
440 //___________________________________________________________________________
441 void GHepParticle::Clear(Option_t * /*option*/)
442 {
443 // implement the Clear(Option_t *) method so that the GHepParticle when is a
444 // member of a GHepRecord, gets deleted properly when calling TClonesArray's
445 // Clear("C")
446 
447  this->CleanUp();
448 }
449 //___________________________________________________________________________
450 void GHepParticle::Print(ostream & stream) const
451 {
452  stream << "\n |";
453  stream << setfill(' ') << setw(14) << this->Name() << " | ";
454  stream << setfill(' ') << setw(14) << this->Pdg() << " | ";
455  stream << setfill(' ') << setw(6) << this->Status() << " | ";
456  stream << setfill(' ') << setw(3) << this->FirstMother() << " | ";
457  stream << setfill(' ') << setw(3) << this->LastMother() << " | ";
458  stream << setfill(' ') << setw(3) << this->FirstDaughter() << " | ";
459  stream << setfill(' ') << setw(3) << this->LastDaughter() << " | ";
460  stream << std::fixed << setprecision(3);
461  stream << setfill(' ') << setw(6) << this->Px() << " | ";
462  stream << setfill(' ') << setw(6) << this->Py() << " | ";
463  stream << setfill(' ') << setw(6) << this->Pz() << " | ";
464  stream << setfill(' ') << setw(6) << this->E() << " | ";
465  stream << setfill(' ') << setw(6) << this->Mass() << " | ";
466 
467  int rescat_code = this->RescatterCode();
468  if( rescat_code != -1 )
469  {
470  stream << setfill(' ') << setw(5) << rescat_code << " | ";
471  }
472 }
473 //___________________________________________________________________________
474 void GHepParticle::Print(Option_t * /*opt*/) const
475 {
476 // implement the TObject's Print(Option_t *) method
477 
478  this->Print(cout);
479 }
480 //___________________________________________________________________________
482 {
483 // Do the comparisons in steps & put the ones that cost most
484 // in the inner-most {}
485 
486  bool same_particle = (this->fPdgCode == p->fPdgCode);
487  bool same_status = (this->fStatus == p->fStatus );
488 
489  if( !same_particle || !same_status ) return false;
490  else {
491  if ( ! this->CompareFamily(p) ) return false;
492  else {
493  if( ! this->CompareMomentum(p) ) return false;
494  else return true;
495  }
496  }
497 }
498 //___________________________________________________________________________
500 {
501  return (this->fPdgCode == p->fPdgCode);
502 }
503 //___________________________________________________________________________
505 {
506  return (this->fStatus == p->fStatus);
507 }
508 //___________________________________________________________________________
510 {
511  bool same_family = (
512  this->fFirstMother == p->fFirstMother &&
513  this->fLastMother == p->fLastMother &&
514  this->fFirstDaughter == p->fFirstDaughter &&
515  this->fLastDaughter == p->fLastDaughter
516  );
517  return same_family;
518 }
519 //___________________________________________________________________________
521 {
522  double dE = TMath::Abs( this->E() - p->E() );
523  double dPx = TMath::Abs( this->Px() - p->Px() );
524  double dPy = TMath::Abs( this->Py() - p->Py() );
525  double dPz = TMath::Abs( this->Pz() - p->Pz() );
526 
527  bool same_momentum =
528  (dE < kPCutOff && dPx < kPCutOff && dPy < kPCutOff && dPz < kPCutOff);
529 
530  return same_momentum;
531 }
532 //___________________________________________________________________________
533 void GHepParticle::Copy(const GHepParticle & particle)
534 {
535  this->SetStatus (particle.Status() );
536  this->SetPdgCode (particle.Pdg() );
537  this->SetRescatterCode (particle.RescatterCode() );
538  this->SetFirstMother (particle.FirstMother() );
539  this->SetLastMother (particle.LastMother() );
540  this->SetFirstDaughter (particle.FirstDaughter() );
541  this->SetLastDaughter (particle.LastDaughter() );
542 
543  this->SetMomentum (*particle.P4());
544  this->SetPosition (*particle.X4());
545 
546  this->fPolzTheta = particle.fPolzTheta;
547  this->fPolzPhi = particle.fPolzPhi;
548 
549  this->fIsBound = particle.fIsBound;
550  this->fRemovalEnergy = particle.fRemovalEnergy;
551 }
552 //___________________________________________________________________________
554 {
555  TParticlePDG * p = PDGLibrary::Instance()->Find(fPdgCode);
556  if(!p) {
557  LOG("GHepParticle", pFATAL)
558  << "\n** You are attempting to insert particle with PDG code = "
559  << fPdgCode << " into the event record."
560  << "\n** This particle can not be found in "
561  << "$GENIE/data/evgen/catalogues/pdg/genie_pdg_table.txt";
562  gAbortingInErr = true;
563  exit(1);
564  }
565 }
566 //___________________________________________________________________________
568 {
569  return (this->Compare(&p));
570 }
571 //___________________________________________________________________________
573 {
574  this->Copy(p);
575  return (*this);
576 }
577 //___________________________________________________________________________
578 
const double kPi
const double kPCutOff
int Z(void) const
void SetFirstMother(int m)
Definition: GHepParticle.h:133
Basic constants.
int RescatterCode(void) const
Definition: GHepParticle.h:66
TLorentzVector * GetX4(void) const
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:52
int status
Definition: fabricate.py:1613
bool IsOnMassShell(void) const
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
double E(void) const
Get energy.
Definition: GHepParticle.h:92
void SetPz(double pz)
double fRemovalEnergy
removal energy for bound nucleons (GeV)
Definition: GHepParticle.h:186
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
int FirstDaughter(void) const
Definition: GHepParticle.h:69
int fLastDaughter
last daughter idx
Definition: GHepParticle.h:181
enum genie::EGHepStatus GHepStatus_t
int fPdgCode
particle PDG code
Definition: GHepParticle.h:175
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:61
void AssertIsKnownParticle(void) const
void SetBound(bool bound)
void SetMomentum(const TLorentzVector &p4)
GHepParticle & operator=(const GHepParticle &p)
double Mass(void) const
Mass that corresponds to the PDG code.
void SetPolarization(double theta, double phi)
int fFirstMother
first mother idx
Definition: GHepParticle.h:178
double fPolzPhi
azimuthal polarization angle (rad)
Definition: GHepParticle.h:185
bool fIsBound
&#39;is it a bound particle?&#39; flag
Definition: GHepParticle.h:187
GHepStatus_t fStatus
particle status
Definition: GHepParticle.h:176
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
double dE
GHepStatus_t Status(void) const
Definition: GHepParticle.h:65
void SetPx(double px)
int fFirstDaughter
first daughter idx
Definition: GHepParticle.h:180
Double_t K
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
bool CompareMomentum(const GHepParticle *p) const
bool Compare(const GHepParticle *p) const
int LastMother(void) const
Definition: GHepParticle.h:68
int Pdg(void) const
Definition: GHepParticle.h:64
int FirstMother(void) const
Definition: GHepParticle.h:67
string Name(void) const
Name that corresponds to the PDG code.
void Clear(Option_t *option)
TLorentzVector * fP4
momentum 4-vector (GeV)
Definition: GHepParticle.h:182
void SetPosition(const TLorentzVector &v4)
int LastDaughter(void) const
Definition: GHepParticle.h:70
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
bool operator==(const GHepParticle &p) const
void SetLastDaughter(int d)
Definition: GHepParticle.h:136
int fRescatterCode
rescattering code
Definition: GHepParticle.h:177
void SetRescatterCode(int code)
Definition: GHepParticle.h:130
void Copy(const GHepParticle &particle)
TLorentzVector * GetP4(void) const
bool ComparePdgCodes(const GHepParticle *p) const
void Print(ostream &stream) const
double Charge(void) const
Chrg that corresponds to the PDG code.
z
Definition: test.py:28
const double kOffShellDm
bool PolzIsSet(void) const
#define pWARN
Definition: Messenger.h:61
bool IsOffMassShell(void) const
OStream cout
Definition: OStream.cxx:6
void GetPolarization(TVector3 &polz)
Definition: inftrees.h:24
bool CompareStatusCodes(const GHepParticle *p) const
void SetRemovalEnergy(double Erm)
TLorentzVector * fX4
position 4-vector (in the target nucleus coordinate system / x,y,z in fm / t=0)
Definition: GHepParticle.h:183
bool CompareFamily(const GHepParticle *p) const
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
void SetLastMother(int m)
Definition: GHepParticle.h:134
bool IsIon(int pdgc)
Definition: PDGUtils.cxx:40
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
ostream & operator<<(ostream &stream, const AlgConfigPool &config_pool)
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:127
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:314
int fLastMother
last mother idx
Definition: GHepParticle.h:179
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:53
double fPolzTheta
polar polarization angle (rad)
Definition: GHepParticle.h:184
int A(void) const
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:64
void SetEnergy(double E)
Float_t e
Definition: plot.C:35
void SetFirstDaughter(int d)
Definition: GHepParticle.h:135
bool gAbortingInErr
Definition: Messenger.cxx:56
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
ClassImp(GHepParticle) namespace genie
void SetPy(double py)
#define pDEBUG
Definition: Messenger.h:64
void SetPdgCode(int c)
double Py(void) const
Get Py.
Definition: GHepParticle.h:90