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