GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
KLVOxygenIBDPXSec.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  Corey Reed <cjreed \at nikhef.nl>
7  Nikhef
8  */
9 //____________________________________________________________________________
10 
11 #include <TSpline.h>
12 #include <TGraph.h>
13 
21 
22 using namespace genie;
23 using namespace genie::constants;
24 
25 const double KLVOxygenIBDPXSec::kO16NubarThr = 0.0114; // GeV
26 const double KLVOxygenIBDPXSec::kO16NuMinE = 0.0150; // GeV
27 const double KLVOxygenIBDPXSec::kMaxE = 0.1000; // GeV
28 
30 
31 //____________________________________________________________________________
34  fXsplNue(0),
35  fXsplNuebar(0)
36 {
37 
38 }
39 //____________________________________________________________________________
41  XSecAlgorithmI("genie::KLVOxygenIBDPXSec", config),
42  fXsplNue(0),
43  fXsplNuebar(0)
44 {
45 
46 }
47 //____________________________________________________________________________
49 {
50  delete fXsplNue;
51  delete fXsplNuebar;
52 }
53 //____________________________________________________________________________
54 double KLVOxygenIBDPXSec::XSec(const Interaction * interaction,
55  KinePhaseSpace_t /* kps */) const
56 {
57  // compute the differential cross section ds/dt
58  // currently not implemented (only total)
59 
60  if(! this -> ValidProcess (interaction) ) return 0.;
61  if(! this -> ValidKinematics (interaction) ) return 0.;
62 
63  LOG("KLVOxygen", pWARN)
64  << "*** No differential cross section calculation is implemented yet";
65 
66  return 1;
67 }
68 //____________________________________________________________________________
70 {
71  // make a new xsec spline from the KLV paper's calculation
72  // for the 16O + nu_e_bar reaction
73  // remove any spline that might already exist
74 
75  delete fXsplNuebar;
76 
77  static const Int_t npts_nuebar = 21;
78  static const Double_t Evnuebar[npts_nuebar] = {
80  15.0e-3, 17.5e-3, 20.0e-3, 22.5e-3, 25.0e-3,
81  27.5e-3, 30.0e-3, 32.5e-3, 35.0e-3, 37.5e-3,
82  40.0e-3, 45.0e-3, 50.0e-3, 55.0e-3, 60.0e-3,
83  65.0e-3, 70.0e-3, 80.0e-3, 90.0e-3, 100.0e-3
84  };
85  static const Double_t xsunit = 1e-42*units::cm2;
86  static const Double_t Onuebar[npts_nuebar] = {
87  0,
88  2.53e-2*xsunit, 7.27e-2*xsunit, 1.81e-1*xsunit, 4.21e-1*xsunit, 8.90e-1*xsunit,
89  1.69*xsunit, 2.94*xsunit, 4.76*xsunit, 7.26*xsunit, 1.06e1*xsunit,
90  1.48e1*xsunit, 2.64e1*xsunit, 4.29e1*xsunit, 6.46e1*xsunit, 9.17e1*xsunit,
91  1.25e2*xsunit, 1.63e2*xsunit, 2.57e2*xsunit, 3.77e2*xsunit, 5.18e2*xsunit
92  };
93  // make spline via dummy TGraph because TSpline3's ctor isn't const correct
94  const TGraph dummy(npts_nuebar,Evnuebar,Onuebar);
95  fXsplNuebar = new TSpline3("16O_nu_e_bar_xsec",&dummy);
96  fXsplNuebar->SetNpx(500);
97 }
98 //____________________________________________________________________________
100 {
101  // make a new xsec spline from the KLV paper's calculation
102  // for the 16O + nu_e reaction
103  // remove any spline that might already exist
104 
105  delete fXsplNue;
106 
107  static const Int_t npts_nue = 20;
108  static const Double_t Evnue[npts_nue] = {
109  15.0e-3, 17.5e-3, 20.0e-3, 22.5e-3, 25.0e-3,
110  27.5e-3, 30.0e-3, 32.5e-3, 35.0e-3, 37.5e-3,
111  40.0e-3, 45.0e-3, 50.0e-3, 55.0e-3, 60.0e-3,
112  65.0e-3, 70.0e-3, 80.0e-3, 90.0e-3, 100.0e-3
113  };
114  static const Double_t xsunit = 1e-42*units::cm2;
115  static const Double_t Onue[npts_nue] = {
116  1.56e-6*xsunit, 8.42e-4*xsunit, 7.26e-3*xsunit, 3.99e-2*xsunit, 1.77e-1*xsunit,
117  5.23e-1*xsunit, 1.25*xsunit, 2.58*xsunit, 4.76*xsunit, 8.05*xsunit,
118  1.28e1*xsunit, 2.76e1*xsunit, 5.21e1*xsunit, 8.89e1*xsunit, 1.41e2*xsunit,
119  2.12e2*xsunit, 3.02e2*xsunit, 5.52e2*xsunit, 8.92e2*xsunit, 1.32e3*xsunit
120  };
121  const TGraph dummy(npts_nue,Evnue,Onue);
122  fXsplNue = new TSpline3("16O_nu_e_xsec",&dummy);
123  fXsplNue->SetNpx(500);
124 }
125 //____________________________________________________________________________
126 double KLVOxygenIBDPXSec::Integral(const Interaction * interaction) const
127 {
128  // Compute the total cross section for a free nucleon target
129 
130  assert(interaction!=0);
131  if(! this -> ValidProcess (interaction) ) return 0.;
132  if(! this -> ValidKinematics (interaction) ) return 0.;
133 
134  const InitialState & init_state = interaction -> InitState();
135  const double Ev = init_state.ProbeE(kRfHitNucRest);
136  const int prbpdg = init_state.ProbePdg();
137 
138  double xsec = 0;
139 
140  if (pdg::IsNuE(prbpdg)) {
141  assert(fXsplNue!=0);
142  xsec = fXsplNue->Eval(Ev);
143  } else if (pdg::IsAntiNuE(prbpdg)) {
144  assert(fXsplNuebar!=0);
145  xsec = fXsplNuebar->Eval(Ev);
146  } else {
147  LOG("KLVOxygen", pERROR) << "*** <Integral> Probe has invalid pdg ["
148  << init_state.ProbePdg() << "]";
149  }
150 
151  return xsec;
152 }
153 //____________________________________________________________________________
154 bool KLVOxygenIBDPXSec::ValidProcess(const Interaction * interaction) const
155 {
156  if(interaction->TestBit(kISkipProcessChk)) return true;
157 
158  // should be IBD and either nu_e + O16 or anu_e + O16
159  if (interaction->ProcInfo().IsInverseBetaDecay()) {
160 
161  const InitialState & init_state = interaction -> InitState();
162  if (init_state.TgtPdg() == kPdgTgtO16) {
163 
164  if ( (pdg::IsNuE(init_state.ProbePdg())) ||
165  (pdg::IsAntiNuE(init_state.ProbePdg())) ) {
166 
167  return true;
168 
169  } else {
170  LOG("KLVOxygen", pERROR) << "*** Probe has invalid pdg ["
171  << init_state.ProbePdg() << "]";
172  }
173 
174  } else {
175  LOG("KLVOxygen", pERROR) << "*** Target has pdg ["
176  << init_state.TgtPdg()
177  << "], not 16O ("
178  << kPdgTgtO16 << ")!";
179  }
180 
181  }
182 
183  return false;
184 }
185 //____________________________________________________________________________
186 bool KLVOxygenIBDPXSec::ValidKinematics(const Interaction* interaction) const
187 {
188  // check energy range
189 
190  if(interaction->TestBit(kISkipKinematicChk)) return true;
191 
192  const InitialState & init_state = interaction -> InitState();
193  const double Ev = init_state.ProbeE(kRfHitNucRest);
194  if ( pdg::IsNuE(init_state.ProbePdg()) ) {
195  if ( (Ev < kO16NuMinE) || (Ev > kMaxE) ) {
196  LOG("KLVOxygen", pERROR) << "*** Ev=" << Ev
197  << " outside range ("
198  << kO16NuMinE << ", " << kMaxE << ")!";
199  return false;
200  }
201  } else if ( pdg::IsAntiNuE(init_state.ProbePdg()) ) {
202  if ( (Ev < kO16NubarThr) || (Ev > kMaxE) ) {
203  LOG("KLVOxygen", pERROR) << "*** Ev=" << Ev
204  << " outside range ("
205  << kO16NubarThr << ", " << kMaxE << ")!";
206  return false;
207  }
208  }
209 
210  const KPhaseSpace & kps = interaction->PhaseSpace();
211  return kps.IsAboveThreshold();
212 }
213 //____________________________________________________________________________
215 {
216  Algorithm::Configure(config);
217  this->LoadConfig();
218 }
219 //____________________________________________________________________________
220 void KLVOxygenIBDPXSec::Configure(string config)
221 {
222  Algorithm::Configure(config);
223  this->LoadConfig();
224 }
225 //____________________________________________________________________________
227 {
228  // make splines
230  MakeNuESpline();
231 }
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
#define pERROR
Definition: Messenger.h:59
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
TSpline3 * fXsplNuebar
a spline around the 16O+nu_e xsec points listed in the reference paper
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
static const double kMaxE
bool IsNuE(int pdgc)
Definition: PDGUtils.cxx:158
bool IsInverseBetaDecay(void) const
enum genie::EKinePhaseSpace KinePhaseSpace_t
An implementation of the neutrino - Oxygen16 cross section.
static const double kO16NuMinE
const int kPdgTgtO16
Definition: PDGCodes.h:203
Summary information for an interaction.
Definition: Interaction.h:56
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
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int ProbePdg(void) const
Definition: InitialState.h:64
Kinematical phase space.
Definition: KPhaseSpace.h:33
static const double kO16NubarThr
#define pWARN
Definition: Messenger.h:60
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
ClassImp(CacheBranchFx)
int TgtPdg(void) const
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
void Configure(const Registry &config)
double Integral(const Interaction *i) const
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
bool IsAntiNuE(int pdgc)
Definition: PDGUtils.cxx:173
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48