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