GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
InitialState.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  Changes required to implement the GENIE Boosted Dark Matter module
10  were installed by Josh Berger (Univ. of Wisconsin)
11 
12  CMEnergy() method added by Andy Furmanski (Univ. of Manchester)
13  and Joe Johnston (Univ of Pittsburgh)
14 */
15 //____________________________________________________________________________
16 
17 #include <cassert>
18 #include <sstream>
19 #include <iomanip>
20 
21 #include <TRootIOCtor.h>
22 
27 
28 using namespace genie;
29 
30 using std::endl;
31 using std::setprecision;
32 using std::setw;
33 using std::ostringstream;
34 
36 
37 //____________________________________________________________________________
38 namespace genie {
39  ostream & operator << (ostream & stream, const InitialState & init_state)
40  {
41  init_state.Print(stream);
42  return stream;
43  }
44 }
45 //___________________________________________________________________________
47 TObject()
48 {
49  this->Init();
50 }
51 //___________________________________________________________________________
52 InitialState::InitialState(int target_pdgc, int probe_pdgc) :
53 TObject()
54 {
55  this->Init(target_pdgc, probe_pdgc);
56 }
57 //___________________________________________________________________________
58 InitialState::InitialState(int Z, int A, int probe_pdgc) :
59 TObject()
60 {
61  int target_pdgc = pdg::IonPdgCode(A,Z);
62  this->Init(target_pdgc, probe_pdgc);
63 }
64 //___________________________________________________________________________
65 InitialState::InitialState(const Target & tgt, int probe_pdgc) :
66 TObject()
67 {
68  int target_pdgc = tgt.Pdg();
69  this->Init(target_pdgc, probe_pdgc);
70 }
71 //___________________________________________________________________________
73 TObject()
74 {
75  this->Init();
76  this->Copy(init_state);
77 }
78 //___________________________________________________________________________
80 TObject(),
81 fProbePdg(0),
82 fTgt(0),
83 fProbeP4(0),
84 fTgtP4(0)
85 {
86 
87 }
88 //___________________________________________________________________________
90 {
91  this->CleanUp();
92 }
93 //___________________________________________________________________________
95 {
96  fProbePdg = 0;
97  fTgt = new Target();
98  fProbeP4 = new TLorentzVector(0, 0, 0, 0);
99  fTgtP4 = new TLorentzVector(0, 0, 0, 0);
100 }
101 //___________________________________________________________________________
102 void InitialState::Init(int target_pdgc, int probe_pdgc)
103 {
104  TParticlePDG * t = PDGLibrary::Instance()->Find(target_pdgc);
105  TParticlePDG * p = PDGLibrary::Instance()->Find(probe_pdgc );
106 
107  assert(t && p);
108 
109  double m = p->Mass();
110  double M = t->Mass();
111 
112  fProbePdg = probe_pdgc;
113  fTgt = new Target(target_pdgc);
114  fProbeP4 = new TLorentzVector(0, 0, 0, m);
115  fTgtP4 = new TLorentzVector(0, 0, 0, M);
116 }
117 //___________________________________________________________________________
119 {
120  delete fTgt;
121  delete fProbeP4;
122  delete fTgtP4;
123 }
124 //___________________________________________________________________________
126 {
127  this->CleanUp();
128  this->Init();
129 }
130 //___________________________________________________________________________
131 void InitialState::Copy(const InitialState & init_state)
132 {
133  fProbePdg = init_state.fProbePdg;
134 
135  fTgt->Copy(*init_state.fTgt);
136 
137  this -> SetProbeP4 ( *init_state.fProbeP4 );
138  this -> SetTgtP4 ( *init_state.fTgtP4 );
139 }
140 //___________________________________________________________________________
141 int InitialState::TgtPdg(void) const
142 {
143  assert(fTgt);
144  return fTgt->Pdg();
145 }
146 //___________________________________________________________________________
147 TParticlePDG * InitialState::Probe(void) const
148 {
149  TParticlePDG * p = PDGLibrary::Instance()->Find(fProbePdg);
150  return p;
151 }
152 //___________________________________________________________________________
153 void InitialState::SetPdgs(int tgt_pdgc, int probe_pdgc)
154 {
155  this->CleanUp();
156  this->Init(tgt_pdgc, probe_pdgc);
157 }
158 //___________________________________________________________________________
159 void InitialState::SetTgtPdg(int tgt_pdgc)
160 {
161  int probe_pdgc = this->ProbePdg();
162 
163  this->CleanUp();
164  this->Init(tgt_pdgc, probe_pdgc);
165 }
166 //___________________________________________________________________________
167 void InitialState::SetProbePdg(int probe_pdgc)
168 {
169  int tgt_pdgc = this->TgtPdg();
170 
171  this->CleanUp();
172  this->Init(tgt_pdgc, probe_pdgc);
173 }
174 //___________________________________________________________________________
176 {
177  fProbeP4 -> SetE ( E );
178  fProbeP4 -> SetPx ( 0.);
179  fProbeP4 -> SetPy ( 0.);
180  fProbeP4 -> SetPz ( E );
181 }
182 //___________________________________________________________________________
183 void InitialState::SetProbeP4(const TLorentzVector & P4)
184 {
185  fProbeP4 -> SetE ( P4.E() );
186  fProbeP4 -> SetPx ( P4.Px() );
187  fProbeP4 -> SetPy ( P4.Py() );
188  fProbeP4 -> SetPz ( P4.Pz() );
189 }
190 //___________________________________________________________________________
191 void InitialState::SetTgtP4(const TLorentzVector & P4)
192 {
193  fTgtP4 -> SetE ( P4.E() );
194  fTgtP4 -> SetPx ( P4.Px() );
195  fTgtP4 -> SetPy ( P4.Py() );
196  fTgtP4 -> SetPz ( P4.Pz() );
197 }
198 //___________________________________________________________________________
199 bool InitialState::IsNuP(void) const
200 {
201  int prob = fProbePdg;
202  int nucl = fTgt->HitNucPdg();
203  bool isvp = pdg::IsNeutrino(prob) && pdg::IsProton(nucl);
204 
205  return isvp;
206 }
207 //___________________________________________________________________________
208 bool InitialState::IsNuN(void) const
209 {
210  int prob = fProbePdg;
211  int nucl = fTgt->HitNucPdg();
212  bool isvn = pdg::IsNeutrino(prob) && pdg::IsNeutron(nucl);
213 
214  return isvn;
215 }
216 //___________________________________________________________________________
217 bool InitialState::IsNuBarP(void) const
218 {
219  int prob = fProbePdg;
220  int nucl = fTgt->HitNucPdg();
221  bool isvbp = pdg::IsAntiNeutrino(prob) && pdg::IsProton(nucl);
222 
223  return isvbp;
224 }
225 //___________________________________________________________________________
226 bool InitialState::IsNuBarN(void) const
227 {
228  int prob = fProbePdg;
229  int nucl = fTgt->HitNucPdg();
230  bool isvbn = pdg::IsAntiNeutrino(prob) && pdg::IsNeutron(nucl);
231 
232  return isvbn;
233 }
234 //___________________________________________________________________________
235 bool InitialState::IsDMP(void) const
236 {
237 // Check if DM - proton interaction
238  int prob = fProbePdg;
239  int nucl = fTgt->HitNucPdg();
240  bool isdp = pdg::IsDarkMatter(prob) && pdg::IsProton(nucl);
241 
242  return isdp;
243 }
244 //___________________________________________________________________________
245 bool InitialState::IsDMN(void) const
246 {
247 // Check if DM - neutron interaction
248  int prob = fProbePdg;
249  int nucl = fTgt->HitNucPdg();
250  bool isdn = pdg::IsDarkMatter(prob) && pdg::IsNeutron(nucl);
251 
252  return isdn;
253 }
254 //___________________________________________________________________________
255 bool InitialState::IsDMBP(void) const
256 {
257 // Check if DM - proton interaction
258  int prob = fProbePdg;
259  int nucl = fTgt->HitNucPdg();
260  bool isdp = pdg::IsAntiDarkMatter(prob) && pdg::IsProton(nucl);
261 
262  return isdp;
263 }
264 //___________________________________________________________________________
265 bool InitialState::IsDMBN(void) const
266 {
267 // Check if DM - neutron interaction
268  int prob = fProbePdg;
269  int nucl = fTgt->HitNucPdg();
270  bool isdn = pdg::IsAntiDarkMatter(prob) && pdg::IsNeutron(nucl);
271 
272  return isdn;
273 }
274 //___________________________________________________________________________
275 TLorentzVector * InitialState::GetTgtP4(RefFrame_t ref_frame) const
276 {
277 // Return the target 4-momentum in the specified reference frame
278 // Note: the caller adopts the TLorentzVector object
279 
280  switch (ref_frame) {
281 
282  //------------------ NUCLEAR TARGET REST FRAME:
283  case (kRfTgtRest) :
284  {
285  // for now make sure that the target nucleus is always at
286  // rest and it is only the struck nucleons that can move:
287  // so the [target rest frame] = [LAB frame]
288 
289  return this->GetTgtP4(kRfLab);
290  }
291 
292  //------------------ STRUCK NUCLEON REST FRAME:
293  case (kRfHitNucRest) :
294  {
295  // make sure that 'struck nucleon' properties were set in
296  // the nuclear target object
297  assert(fTgt->HitNucIsSet());
298  TLorentzVector * pnuc4 = fTgt->HitNucP4Ptr();
299 
300  // compute velocity vector (px/E, py/E, pz/E)
301  double bx = pnuc4->Px() / pnuc4->Energy();
302  double by = pnuc4->Py() / pnuc4->Energy();
303  double bz = pnuc4->Pz() / pnuc4->Energy();
304 
305  // BOOST
306  TLorentzVector * p4 = new TLorentzVector(*fTgtP4);
307  p4->Boost(-bx,-by,-bz);
308 
309  return p4;
310  break;
311  }
312  //------------------ LAB:
313  case (kRfLab) :
314  {
315  TLorentzVector * p4 = new TLorentzVector(*fTgtP4);
316  return p4;
317  break;
318  }
319  default:
320  LOG("Interaction", pERROR) << "Uknown reference frame";
321  }
322  return 0;
323 }
324 //___________________________________________________________________________
325 TLorentzVector * InitialState::GetProbeP4(RefFrame_t ref_frame) const
326 {
327 // Return the probe 4-momentum in the specified reference frame
328 // Note: the caller adopts the TLorentzVector object
329 
330  switch (ref_frame) {
331 
332  //------------------ NUCLEAR TARGET REST FRAME:
333  case (kRfTgtRest) :
334  {
335  // for now assure that the target nucleus is always at rest
336  // and it is only the struck nucleons that can move:
337  // so the [target rest frame] = [LAB frame]
338 
339  TLorentzVector * p4 = new TLorentzVector(*fProbeP4);
340  return p4;
341  }
342 
343  //------------------ STRUCK NUCLEON REST FRAME:
344  case (kRfHitNucRest) :
345  {
346  // make sure that 'struck nucleon' properties were set in
347  // the nuclear target object
348 
349  assert( fTgt->HitNucP4Ptr() != 0 );
350 
351  TLorentzVector * pnuc4 = fTgt->HitNucP4Ptr();
352 
353  // compute velocity vector (px/E, py/E, pz/E)
354 
355  double bx = pnuc4->Px() / pnuc4->Energy();
356  double by = pnuc4->Py() / pnuc4->Energy();
357  double bz = pnuc4->Pz() / pnuc4->Energy();
358 
359  // BOOST
360 
361  TLorentzVector * p4 = new TLorentzVector(*fProbeP4);
362 
363  p4->Boost(-bx,-by,-bz);
364 
365  return p4;
366 
367  break;
368  }
369  //------------------ LAB:
370  case (kRfLab) :
371  {
372  TLorentzVector * p4 = new TLorentzVector(*fProbeP4);
373  return p4;
374 
375  break;
376  }
377  default:
378 
379  LOG("Interaction", pERROR) << "Uknown reference frame";
380  }
381  return 0;
382 }
383 //___________________________________________________________________________
384 double InitialState::ProbeE(RefFrame_t ref_frame) const
385 {
386  TLorentzVector * p4 = this->GetProbeP4(ref_frame);
387  double E = p4->Energy();
388 
389  delete p4;
390  return E;
391 }
392 
393 //___________________________________________________________________________
395 {
396  TLorentzVector * k4 = this->GetProbeP4(kRfLab);
397  TLorentzVector * p4 = fTgt->HitNucP4Ptr();
398 
399  *k4 += *p4; // now k4 represents centre-of-mass 4-momentum
400  double s = k4->Dot(*k4); // dot-product with itself
401  double E = TMath::Sqrt(s);
402 
403  delete k4;
404 
405  return E;
406 }
407 //___________________________________________________________________________
408 string InitialState::AsString(void) const
409 {
410 // Code-ify the interaction in a string to be used as (part of a) keys
411 // Template:
412 // nu_pdg:code;tgt-pdg:code;
413 
414  ostringstream init_state;
415 
416  if (this->Probe()->Mass() > 0) {
417  init_state << "dm_mass:" << this->Probe()->Mass() << ";";
418  }
419  else {
420  init_state << "nu-pdg:" << this->ProbePdg() << ";";
421  }
422  init_state << "tgt-pdg:" << this->Tgt().Pdg() << ";";
423 
424  return init_state.str();
425 }
426 //___________________________________________________________________________
427 void InitialState::Print(ostream & stream) const
428 {
429  stream << "[-] [Init-State] " << endl;
430 
431  stream << " |--> probe : "
432  << "PDG-code = " << fProbePdg
433  << " (" << this->Probe()->GetName() << ")" << endl;
434 
435  stream << " |--> nucl. target : "
436  << "Z = " << fTgt->Z()
437  << ", A = " << fTgt->A()
438  << ", PDG-Code = " << fTgt->Pdg();
439 
440  TParticlePDG * tgt = PDGLibrary::Instance()->Find( fTgt->Pdg() );
441  if(tgt) {
442  stream << " (" << tgt->GetName() << ")";
443  }
444  stream << endl;
445 
446  stream << " |--> hit nucleon : ";
447  int nuc_pdgc = fTgt->HitNucPdg();
448 
449  if ( pdg::IsNeutronOrProton(nuc_pdgc) ) {
450  TParticlePDG * p = PDGLibrary::Instance()->Find(nuc_pdgc);
451  stream << "PDC-Code = " << nuc_pdgc << " (" << p->GetName() << ")";
452  } else {
453  stream << "no set";
454  }
455  stream << endl;
456 
457  stream << " |--> hit quark : ";
458  int qrk_pdgc = fTgt->HitQrkPdg();
459 
460  if ( pdg::IsQuark(qrk_pdgc) || pdg::IsAntiQuark(qrk_pdgc)) {
461  TParticlePDG * p = PDGLibrary::Instance()->Find(qrk_pdgc);
462  stream << "PDC-Code = " << qrk_pdgc << " (" << p->GetName() << ") ";
463  stream << (fTgt->HitSeaQrk() ? "[sea]" : "[valence]");
464  } else {
465  stream << "no set";
466  }
467  stream << endl;
468 
469  stream << " |--> probe 4P : "
470  << "(E = " << setw(12) << setprecision(6) << fProbeP4->E()
471  << ", Px = " << setw(12) << setprecision(6) << fProbeP4->Px()
472  << ", Py = " << setw(12) << setprecision(6) << fProbeP4->Py()
473  << ", Pz = " << setw(12) << setprecision(6) << fProbeP4->Pz()
474  << ")"
475  << endl;
476  stream << " |--> target 4P : "
477  << "(E = " << setw(12) << setprecision(6) << fTgtP4->E()
478  << ", Px = " << setw(12) << setprecision(6) << fTgtP4->Px()
479  << ", Py = " << setw(12) << setprecision(6) << fTgtP4->Py()
480  << ", Pz = " << setw(12) << setprecision(6) << fTgtP4->Pz()
481  << ")"
482  << endl;
483 
484  if ( pdg::IsNeutronOrProton(nuc_pdgc) ) {
485 
486  TLorentzVector * nuc_p4 = fTgt->HitNucP4Ptr();
487 
488  stream << " |--> nucleon 4P : "
489  << "(E = " << setw(12) << setprecision(6) << nuc_p4->E()
490  << ", Px = " << setw(12) << setprecision(6) << nuc_p4->Px()
491  << ", Py = " << setw(12) << setprecision(6) << nuc_p4->Py()
492  << ", Pz = " << setw(12) << setprecision(6) << nuc_p4->Pz()
493  << ")";
494  }
495 }
496 //___________________________________________________________________________
497 bool InitialState::Compare(const InitialState & init_state) const
498 {
499  int probe = init_state.ProbePdg();
500  const Target & target = init_state.Tgt();
501 
502  bool equal = (fProbePdg == probe) && (*fTgt == target);
503 
504  return equal;
505 }
506 //___________________________________________________________________________
507 bool InitialState::operator == (const InitialState & init_state) const
508 {
509  return this->Compare(init_state);
510 }
511 //___________________________________________________________________________
513 {
514  this->Copy(init_state);
515  return (*this);
516 }
517 //___________________________________________________________________________
void SetTgtP4(const TLorentzVector &P4)
bool HitSeaQrk(void) const
Definition: Target.cxx:299
bool IsNuBarN(void) const
is anti-neutrino + neutron?
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
void SetProbeP4(const TLorentzVector &P4)
#define pERROR
Definition: Messenger.h:59
void SetTgtPdg(int pdg_code)
int HitNucPdg(void) const
Definition: Target.cxx:304
int A(void) const
Definition: Target.h:70
int HitQrkPdg(void) const
Definition: Target.cxx:242
bool IsNuN(void) const
is neutrino + neutron?
int Pdg(void) const
Definition: Target.h:71
static constexpr double s
Definition: Units.h:95
bool IsAntiQuark(int pdgc)
Definition: PDGUtils.cxx:258
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:127
TParticlePDG * Probe(void) const
double Mass(Resonance_t res)
resonance mass (GeV)
void SetProbeE(double E)
bool Compare(const InitialState &init_state) const
enum genie::ERefFrame RefFrame_t
bool operator==(const InitialState &i) const
equal?
bool IsDMN(void) const
is dark matter + neutron?
bool IsDMBP(void) const
is anti-dark matter + proton?
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
bool IsNuP(void) const
is neutrino + proton?
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool IsAntiDarkMatter(int pdgc)
Definition: PDGUtils.cxx:133
void Copy(const InitialState &init_state)
static constexpr double A
Definition: Units.h:74
Target * fTgt
nuclear target
Definition: InitialState.h:110
TLorentzVector * fProbeP4
probe 4-momentum in LAB-frame
Definition: InitialState.h:111
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int ProbePdg(void) const
Definition: InitialState.h:64
string AsString(void) const
TLorentzVector * fTgtP4
nuclear target 4-momentum in LAB-frame
Definition: InitialState.h:112
void SetPdgs(int tgt_pdgc, int probe_pdgc)
int Z(void) const
Definition: Target.h:68
bool IsDMP(void) const
is dark matter + proton?
int fTgt
bool IsDMBN(void) const
is anti-dark matter + neutron?
void Print(ostream &stream) const
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
bool HitNucIsSet(void) const
Definition: Target.cxx:283
ostream & operator<<(ostream &stream, const AlgConfigPool &config_pool)
bool IsQuark(int pdgc)
Definition: PDGUtils.cxx:250
TLorentzVector * GetTgtP4(RefFrame_t rf=kRfLab) const
ClassImp(CacheBranchFx)
double CMEnergy() const
centre-of-mass energy (sqrt s)
int TgtPdg(void) const
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:351
bool IsNuBarP(void) const
is anti-neutrino + proton?
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
InitialState & operator=(const InitialState &i)
copy
void Copy(const Target &t)
Definition: Target.cxx:116
int fProbePdg
probe PDG code
Definition: InitialState.h:109
void SetProbePdg(int pdg_code)
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
const Target & Tgt(void) const
Definition: InitialState.h:66
double ProbeE(RefFrame_t rf) const
static constexpr double m
Definition: Units.h:71
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
Initial State information.
Definition: InitialState.h:48