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