GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
EvtLibPXSec.cxx
Go to the documentation of this file.
2 #include "Tools/EvtLib/Utils.h"
3 
5 
9 
10 #include "TFile.h"
11 #include "TGraph.h"
12 
13 using namespace genie;
14 using namespace genie::evtlib;
15 
16 //____________________________________________________________________________
17 EvtLibPXSec::EvtLibPXSec() :
18 XSecAlgorithmI("genie::evtlib::EvtLibPXSec")
19 {
20 
21 }
22 
23 //____________________________________________________________________________
24 EvtLibPXSec::EvtLibPXSec(string config) :
25 XSecAlgorithmI("genie::evtlib::EvtLibPXSec", config)
26 {
27 
28 }
29 
30 //____________________________________________________________________________
32 {
33  ClearXSecs();
34 }
35 
36 //____________________________________________________________________________
37 double EvtLibPXSec::XSec(const Interaction* /*in*/,
38  KinePhaseSpace_t /*kps*/) const
39 {
40  LOG("ELI", pFATAL) << "This point should not be reached" ;
41  abort();
42 }
43 
44 //____________________________________________________________________________
45 double EvtLibPXSec::Integral(const Interaction* in) const
46 {
47  TGraph* g = GetXSec(in);
48  if(!g) return 0; // Reason already printed
49 
50  const InitialState& init_state = in->InitState();
51  const double E = init_state.ProbeE(kRfLab);
52 
53  // Units of the cross-section graph are expected to be 10^-38 cm^2 / nucleus
54  return g->Eval(E) * 1e-38 * genie::units::cm2;
55 }
56 
57 //____________________________________________________________________________
58 bool EvtLibPXSec::ValidProcess(const Interaction* /*in*/) const
59 {
60  return true;
61 }
62 
63 //____________________________________________________________________________
64 void EvtLibPXSec::Configure(const Registry& config)
65 {
66  Algorithm::Configure(config);
67  LoadXSecs();
68 }
69 
70 //____________________________________________________________________________
71 void EvtLibPXSec::Configure(string config)
72 {
73  Algorithm::Configure(config);
74  LoadXSecs();
75 }
76 
77 //____________________________________________________________________________
79 {
80  for(auto it: fXSecs) delete it.second;
81  fXSecs.clear();
82 }
83 
84 //____________________________________________________________________________
86 {
87  ClearXSecs();
88 
89  std::string libPath;
90  GetParam("EventLibraryPath", libPath);
91  Expand(libPath);
92 
93  PDGLibrary* pdglib = PDGLibrary::Instance();
94 
95  TFile fin(libPath.c_str());
96  if(fin.IsZombie()) exit(1);
97 
98  TIter next(fin.GetListOfKeys());
99  while(TObject* dir = next()){
100  const std::string& tgtName = dir->GetName();
101  const TParticlePDG* tgtPart = pdglib->DBase()->GetParticle(tgtName.c_str());
102  if(!tgtPart){
103  LOG("ELI", pWARN) << "Unknown nucleus " << tgtName
104  << " found in " << libPath
105  << " -- skipping";
106  continue;
107  }
108 
109  for(int pdg: {kPdgNuE, kPdgAntiNuE,
112 
113  for(bool iscc: {true, false}){
114  // NCs should be the same for all flavours. Use nu_mu as a convention
115  // internal to this code to index into the xsecs map.
116  if(!iscc && abs(pdg) != kPdgNuMu) continue;
117 
118  std::string nuName = pdglib->Find(pdg)->GetName();
119  if(!iscc) nuName = pdg::IsAntiNeutrino(pdg) ? "nu_bar" : "nu";
120 
121  const std::string graphName =
122  TString::Format("%s/%s/%s/xsec",
123  tgtName.c_str(),
124  iscc ? "cc" : "nc",
125  nuName.c_str()).Data();
126 
127  const Key key(tgtPart->PdgCode(), pdg, iscc);
128 
129  TGraph* g = (TGraph*)fin.Get(graphName.c_str());
130  if(!g){
131  LOG("ELI", pINFO) << graphName << " not found in "
132  << libPath << " -- skipping";
133  continue;
134  }
135 
136  fXSecs[key] = new TGraph(*g);
137  } // end for iscc
138  } // end for pdg
139  } // end for dir
140 }
141 
142 //____________________________________________________________________________
143 TGraph* EvtLibPXSec::GetXSec(const Interaction* in) const
144 {
145  const InitialState& init_state = in->InitState();
146 
147  const ProcessInfo& proc = in->ProcInfo();
148  if(!proc.IsWeakCC() && !proc.IsWeakNC()){
149  LOG("ELI", pINFO) << "Skipping unknown process " << proc;
150  return 0;
151  }
152 
153  int pdg = init_state.ProbePdg();
154  // Use nu_mu for NC as a convention internal to this code to index into the
155  // xsec graphs map.
156  if(proc.IsWeakNC()){
157  if(pdg > 0) pdg = kPdgNuMu; else pdg = kPdgAntiNuMu;
158  }
159 
160  const Key key(init_state.TgtPdg(), pdg, proc.IsWeakCC());
161 
162  auto it = fXSecs.find(key);
163  if(it == fXSecs.end()){
164  LOG("ELI", pINFO) << "Skipping process without xsec " << key;
165  return 0;
166  }
167 
168  return it->second;
169 }
Cross Section Calculation Interface.
bool IsWeakCC(void) const
const int kPdgNuE
Definition: PDGCodes.h:28
static constexpr double g
Definition: Units.h:144
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Definition: EvtLibPXSec.cxx:58
const int kPdgAntiNuE
Definition: PDGCodes.h:29
void Expand(std::string &s)
Expand env vars/homedirs/wildcards in s.
Definition: Utils.cxx:8
#define pFATAL
Definition: Messenger.h:56
TDatabasePDG * DBase(void)
Definition: PDGLibrary.cxx:81
const int kPdgNuMu
Definition: PDGCodes.h:30
string dir
enum genie::EKinePhaseSpace KinePhaseSpace_t
void Configure(const Registry &config)
Definition: EvtLibPXSec.cxx:64
Summary information for an interaction.
Definition: Interaction.h:56
bool IsWeakNC(void) const
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
static constexpr double cm2
Definition: Units.h:69
TGraph * GetXSec(const Interaction *in) const
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
std::map< Key, TGraph * > fXSecs
Definition: EvtLibPXSec.h:35
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int ProbePdg(void) const
Definition: InitialState.h:64
#define pINFO
Definition: Messenger.h:62
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
#define pWARN
Definition: Messenger.h:60
double Integral(const Interaction *i) const
Definition: EvtLibPXSec.cxx:45
const int kPdgNuTau
Definition: PDGCodes.h:32
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
Definition: EvtLibPXSec.cxx:37
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
int TgtPdg(void) const
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:146
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Initial State information.
Definition: InitialState.h:48