GENIEGenerator
Loading...
Searching...
No Matches
DMELXSec.cxx
Go to the documentation of this file.
1//____________________________________________________________________________
2/*
3 Copyright (c) 2003-2025, The GENIE Collaboration
4 For the full text of the license visit http://copyright.genie-mc.org
5
6
7 Author: Joshua Berger <jberger \at physics.wisc.edu>
8 University of Wisconsin-Madison
9
10 Costas Andreopoulos <c.andreopoulos \at cern.ch>
11 University of Liverpool
12*/
13//____________________________________________________________________________
14
15#include <TMath.h>
16#include <Math/IFunction.h>
17#include <Math/Integrator.h>
18
21#include "Framework/Conventions/GBuild.h"
28
38
39using namespace genie;
40using namespace genie::constants;
41using namespace genie::controls;
42
43//____________________________________________________________________________
45XSecIntegratorI("genie::DMELXSec")
46{
47
48}
49//____________________________________________________________________________
50DMELXSec::DMELXSec(string config) :
51XSecIntegratorI("genie::DMELXSec", config)
52{
53
54}
55//____________________________________________________________________________
60//____________________________________________________________________________
62 const XSecAlgorithmI * model, const Interaction * in) const
63{
64 LOG("DMELXSec",pDEBUG) << "Beginning integrate";
65 if(! model->ValidProcess(in)) return 0.;
66
67 const KPhaseSpace & kps = in->PhaseSpace();
68 if(!kps.IsAboveThreshold()) {
69 LOG("DMELXSec", pDEBUG) << "*** Below energy threshold";
70 return 0;
71 }
72 Range1D_t rQ2 = kps.Limits(kKVQ2);
73
74 if(rQ2.min<0 || rQ2.max<0) return 0;
75 LOG("DMELXSec", pDEBUG)
76 << "Q2 integration range = (" << rQ2.min << ", " << rQ2.max << ")";
77
78 Interaction * interaction = new Interaction(*in);
79 interaction->SetBit(kISkipProcessChk);
80 interaction->SetBit(kISkipKinematicChk);
81
82 ROOT::Math::IBaseFunctionOneDim * func = new
83 utils::gsl::dXSec_dQ2_E(model, interaction);
84 ROOT::Math::IntegrationOneDim::Type ig_type =
86
87 double abstol = 1; //We mostly care about relative tolerance
88 ROOT::Math::Integrator ig(*func,ig_type,abstol,fGSLRelTol,fGSLMaxEval);
89 double xsec = ig.Integral(rQ2.min, rQ2.max) * (1E-38 * units::cm2);
90
91 //LOG("DMELXSec", pDEBUG) << "XSec[DMEL] (E = " << E << ") = " << xsec;
92
93 delete func;
94 delete interaction;
95
96 return xsec;
97}
98//____________________________________________________________________________
99void DMELXSec::Configure(const Registry & config)
100{
101 Algorithm::Configure(config);
102 this->LoadConfig();
103}
104//____________________________________________________________________________
105void DMELXSec::Configure(string config)
106{
107 Algorithm::Configure(config);
108 this->LoadConfig();
109}
110//____________________________________________________________________________
112{
113 // Get GSL integration type & relative tolerance
114 GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive"));
115 GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 0.01 ) ;
116 int max;
117 GetParamDef( "gsl-max-eval", max, 100000) ;
118 fGSLMaxEval = (unsigned int) max ;
119}
120//____________________________________________________________________________
#define pDEBUG
Definition Messenger.h:63
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
virtual void Configure(const Registry &config)
Definition Algorithm.cxx:62
bool GetParamDef(const RgKey &name, T &p, const T &def) const
void LoadConfig(void)
Definition DMELXSec.cxx:111
virtual ~DMELXSec()
Definition DMELXSec.cxx:56
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
Definition DMELXSec.cxx:61
void Configure(const Registry &config)
Definition DMELXSec.cxx:99
Summary information for an interaction.
Definition Interaction.h:56
const KPhaseSpace & PhaseSpace(void) const
Definition Interaction.h:73
Kinematical phase space.
Definition KPhaseSpace.h:33
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
A simple [min,max] interval for doubles.
Definition Range1.h:43
A registry. Provides the container for algorithm configuration parameters.
Definition Registry.h:65
Cross Section Calculation Interface.
virtual bool ValidProcess(const Interaction *i) const =0
Can this cross section algorithm handle the input process?
string fGSLIntgType
name of GSL numerical integrator
int fGSLMaxEval
GSL max evaluations.
double fGSLRelTol
required relative tolerance (error)
double func(double x, double y)
Basic constants.
Misc GENIE control constants.
static constexpr double cm2
Definition Units.h:69
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition GSLUtils.cxx:23
THE MAIN GENIE PROJECT NAMESPACE
Definition AlgCmp.h:25
@ kKVQ2
Definition KineVar.h:33
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition Interaction.h:48
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition Interaction.h:47