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