NewQELXSec.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: Steven Gardiner <gardiner \at fnal.gov>
8  Fermi National Accelerator Laboratory - Feb 26, 2019
9 
10  For the class documentation see the corresponding header file.
11 
12 */
13 //____________________________________________________________________________
14 
24 
32 #include "Framework/Utils/Range1.h"
37 
38 using namespace genie;
39 using namespace genie::constants;
40 using namespace genie::utils::gsl;
41 
42 //____________________________________________________________________________
43 NewQELXSec::NewQELXSec() : XSecIntegratorI("genie::NewQELXSec")
44 {
45 
46 }
47 //____________________________________________________________________________
49 {
50 
51 }
52 //____________________________________________________________________________
54 {
55  LOG("NewQELXSec",pDEBUG) << "Beginning integrate";
56  if ( !model->ValidProcess(in) ) return 0.;
57 
58  Interaction* interaction = new Interaction( *in );
59  interaction->SetBit( kISkipProcessChk );
60  //interaction->SetBit( kISkipKinematicChk );
61 
62  const NuclearModelI* nucl_model = dynamic_cast<const NuclearModelI*>(
63  model->SubAlg("IntegralNuclearModel") );
64  assert( nucl_model );
65 
67  const VertexGenerator* vtx_gen = dynamic_cast<const VertexGenerator*>(
68  algf->GetAlgorithm(fVertexGenID) );
69  assert( vtx_gen );
70 
71  // Determine the appropriate binding energy mode to use.
72  // The default given here is for the case of a free nucleon.
73  QELEvGen_BindingMode_t bind_mode = kOnShell;
74  Target* tgt = interaction->InitState().TgtPtr();
75  if ( tgt->IsNucleus() ) {
76  std::string bind_mode_str = model->GetConfig()
77  .GetString( "IntegralNucleonBindingMode" );
78  bind_mode = genie::utils::StringToQELBindingMode( bind_mode_str );
79  }
80 
82  interaction, bind_mode, fMinAngleEM);
85 
86  // Switch to using the copy of the interaction in the integrator rather than
87  // the copy that we made in this function
88  delete interaction;
89  interaction = func->GetInteractionPtr();
90 
91  // Also update the pointer to the Target
92  tgt = interaction->InitState().TgtPtr();
93 
94  double abstol = 1e-16; // We mostly care about relative tolerance
95  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
96 
97  // Integration ranges for the lepton COM frame scattering angles (in the
98  // kPSQELEvGen phase space, these are measured with respect to the COM
99  // velocity as observed in the lab frame)
100  Range1D_t cos_theta_0_lim( -1., 1. );
101  Range1D_t phi_0_lim( 0., 2.*kPi );
102 
103  double kine_min[2] = { cos_theta_0_lim.min, phi_0_lim.min };
104  double kine_max[2] = { cos_theta_0_lim.max, phi_0_lim.max };
105 
106  // For a free nucleon target (hit nucleon is at rest in the lab frame), we
107  // don't need to do an MC integration over the initial state variables. In
108  // this case, just set up the nucleon at the origin, on-shell, and at rest,
109  // then integrate over the angles and return the result.
110 
111  // Also use this approach if we're over the "nuclear influence" cutoff
112  // energy for the probe. Beyond the cutoff, the effects of Fermi motion
113  // and the removal energy are assumed to be small enough to be neglected
114  double E_lab_cutoff = model->GetConfig()
115  .GetDouble("IntegralNuclearInfluenceCutoffEnergy");
116 
117  double probeE = interaction->InitState().ProbeE( kRfLab );
118  if ( !tgt->IsNucleus() || probeE > E_lab_cutoff ) {
119  tgt->SetHitNucPosition(0.);
120 
121  if ( tgt->IsNucleus() ) nucl_model->GenerateNucleon(*tgt, 0.);
122  else {
123  nucl_model->SetRemovalEnergy(0.);
124  interaction->SetBit( kIAssumeFreeNucleon );
125  }
126 
127  nucl_model->SetMomentum3( TVector3(0., 0., 0.) );
128  double xsec_total = ig.Integral(kine_min, kine_max);
129  delete func;
130  return xsec_total;
131  }
132 
133  // For a nuclear target, we need to loop over a bunch of nucleons sampled
134  // from the nuclear model (with positions sampled from the vertex generator
135  // to allow for using the local Fermi gas model). The MC estimator for the
136  // total cross section is simply the mean of ig.Integral() for all of the
137  // sampled nucleons.
138  double xsec_sum = 0.;
139  for (int n = 0; n < fNumNucleonThrows; ++n) {
140 
141  // Select a new position for the initial hit nucleon (needed for the local
142  // Fermi gas model, but other than slowing things down a bit, it doesn't
143  // hurt to do this for other models)
144  TVector3 vertex_pos = vtx_gen->GenerateVertex( interaction, tgt->A() );
145  double radius = vertex_pos.Mag();
146  tgt->SetHitNucPosition( radius );
147 
148  // Sample a new nucleon 3-momentum and removal energy (this will be applied
149  // to the nucleon via a call to genie::utils::ComputeFullQELPXSec(), so
150  // there's no need to mess with its 4-momentum here)
151  nucl_model->GenerateNucleon(*tgt, radius);
152 
153  // The initial state variables have all been defined, so integrate over
154  // the final lepton angles.
155  double xsec = ig.Integral(kine_min, kine_max);
156 
157  xsec_sum += xsec;
158  }
159 
160  delete func;
161 
162  // MC estimator of the total cross section is the mean of the xsec values
163  double xsec_mean = xsec_sum / fNumNucleonThrows;
164 
165  return xsec_mean;
166 }
167 //____________________________________________________________________________
169 {
170  Algorithm::Configure(config);
171  this->LoadConfig();
172 }
173 //____________________________________________________________________________
174 void NewQELXSec::Configure(string config)
175 {
176  Algorithm::Configure(config);
177  this->LoadConfig();
178 }
179 //____________________________________________________________________________
181 {
182  // Get GSL integration type & relative tolerance
183  GetParamDef( "gsl-integration-type", fGSLIntgType, std::string("adaptive") ) ;
184  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1e-2 ) ;
185  int max;
186  GetParamDef( "gsl-max-eval", max, 500000 ) ;
187  fGSLMaxEval = static_cast<unsigned int>( max );
188 
189  RgAlg vertexGenID;
190  GetParamDef( "VertexGenAlg", vertexGenID, RgAlg("genie::VertexGenerator", "Default") );
191  fVertexGenID = AlgId( vertexGenID );
192 
193  GetParamDef( "NumNucleonThrows", fNumNucleonThrows, 5000 );
194 
195  // TODO: This is a parameter that may also be specified in the XML
196  // configuration for QELEventGenerator. Avoid duplication here to ensure
197  // consistency.
198  GetParamDef( "SF-MinAngleEMscattering", fMinAngleEM, 0. ) ;
199 }
200 
202  const Interaction* interaction, QELEvGen_BindingMode_t binding_mode, double min_angle_EM)
203  : fXSecModel( xsec_model ), fInteraction( new Interaction(*interaction) ),
204  fHitNucleonBindingMode( binding_mode ), fMinAngleEM( min_angle_EM )
205 {
206  fNuclModel = dynamic_cast<const NuclearModelI*>( fXSecModel->SubAlg("IntegralNuclearModel") );
207  assert( fNuclModel );
208 }
209 
211 {
212  delete fInteraction;
213 }
214 
216 {
217  return fInteraction;
218 }
219 
221 {
222  return *fInteraction;
223 }
224 
225 ROOT::Math::IBaseFunctionMultiDim* genie::utils::gsl::FullQELdXSec::Clone(void) const
226 {
228 }
229 
231 {
232  return 2;
233 }
234 
235 double genie::utils::gsl::FullQELdXSec::DoEval(const double* xin) const
236 {
237  // Elements of "xin"
238  //
239  // element 0: "cos_theta0" = Cosine of theta0, the angle between the COM frame
240  // 3-momentum of the outgoing lepton and the COM frame velocity
241  // as measured in the laboratory frame
242  // element 1: "phi_theta0" = Azimuthal angle of the COM frame 3-momentum of the
243  // outgoing lepton measured with respect to the COM frame
244  // velocity as measured in the laboratory frame
245 
246  double cos_theta0 = xin[0];
247  double phi0 = xin[1];
248 
249  // Dummy storage for the binding energy of the hit nucleon
250  double dummy_Eb = 0.;
251 
252  // Compute the full differential cross section
254  fXSecModel, cos_theta0, phi0, dummy_Eb, fHitNucleonBindingMode, fMinAngleEM, true);
255 
256  return xsec;
257 }
Cross Section Calculation Interface.
const double kPi
TVector3 GenerateVertex(const Interaction *in, double A) const
Basic constants.
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:67
int A(void) const
Definition: Target.h:71
A simple [min,max] interval for doubles.
Definition: Range1.h:43
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
Definition: NewQELXSec.cxx:53
bool IsNucleus(void) const
Definition: Target.cxx:289
std::string fGSLIntgType
Definition: NewQELXSec.h:88
Simple utilities for integrating GSL in the GENIE framework.
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
void SetHitNucPosition(double r)
Definition: Target.cxx:227
Definition: config.py:1
RgDbl GetDouble(RgKey key) const
Definition: Registry.cxx:489
double DoEval(const double *xin) const
Definition: NewQELXSec.cxx:235
double ComputeFullQELPXSec(Interaction *interaction, const NuclearModelI *nucl_model, const XSecAlgorithmI *xsec_model, double cos_theta_0, double phi_0, double &Eb, QELEvGen_BindingMode_t hitNucleonBindingMode, double min_angle_EM=0., bool bind_nucleon=true)
Definition: QELUtils.cxx:94
enum genie::EQELEvGenBindingMode QELEvGen_BindingMode_t
void Configure(const Registry &config)
Definition: NewQELXSec.cxx:168
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
unsigned int fGSLMaxEval
Definition: NewQELXSec.h:90
unsigned int NDim(void) const
Definition: NewQELXSec.cxx:230
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
std::void_t< T > n
const Interaction & GetInteraction() const
Definition: NewQELXSec.cxx:220
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:86
double func(double x, double y)
Double_t radius
void LoadConfig(void)
Definition: NewQELXSec.cxx:180
Double_t xsec[nknots]
Definition: testXsec.C:47
const NuclearModelI * fNuclModel
Definition: NewQELXSec.h:55
QELEvGen_BindingMode_t fHitNucleonBindingMode
Definition: NewQELXSec.h:57
Algorithm ID (algorithm name + configuration set name)
Definition: AlgId.h:35
const XSecAlgorithmI * fXSecModel
Definition: NewQELXSec.h:54
void SetRemovalEnergy(double E) const
Definition: NuclearModelI.h:82
FullQELdXSec(const XSecAlgorithmI *xsec_model, const Interaction *interaction, QELEvGen_BindingMode_t binding_mode, double min_angle_EM)
Definition: NewQELXSec.cxx:201
double max
Definition: Range1.h:54
RgStr GetString(RgKey key) const
Definition: Registry.cxx:496
ifstream in
Definition: comparison.C:7
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
void SetMomentum3(const TVector3 &mom) const
Definition: NuclearModelI.h:78
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
Definition: QELUtils.cxx:195
Target * TgtPtr(void) const
Definition: InitialState.h:68
virtual bool ValidProcess(const Interaction *i) const =0
Can this cross section algorithm handle the input process?
assert(nhit_max >=nhit_nbins)
virtual bool GenerateNucleon(const Target &) const =0
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:53
bool GetParamDef(const RgKey &name, T &p, const T &def) const
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
Float_t e
Definition: plot.C:35
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:66
const XML_Char XML_Content * model
Definition: expat.h:151
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
Definition: NewQELXSec.cxx:225
double fMinAngleEM
Definition: NewQELXSec.h:93
double ProbeE(RefFrame_t rf) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
#define pDEBUG
Definition: Messenger.h:64
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353
enum BeamMode string