GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HEDISPXSec.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  or see $GENIE/LICENSE
6 
7  Author: Alfonso Garcia <alfonsog \at nikhef.nl>
8  NIKHEF
9 
10  For the class documentation see the corresponding header file.
11 
12 */
13 //____________________________________________________________________________
14 
22 
23 #include <TMath.h>
24 
25 using namespace genie;
26 using namespace genie::constants;
27 
28 //____________________________________________________________________________
30 XSecAlgorithmI("genie::HEDISPXSec")
31 {
32 
33 }
34 //____________________________________________________________________________
35 HEDISPXSec::HEDISPXSec(string config) :
36 XSecAlgorithmI("genie::HEDISPXSec", config)
37 {
38 
39 }
40 //____________________________________________________________________________
42 {
43 
44 }
45 //____________________________________________________________________________
47  const Interaction * interaction, KinePhaseSpace_t kps) const
48 {
49 
50  if(! this -> ValidKinematics (interaction) ) return 0.;
51 
52  // Load SF tables
54 
55  // W limits are computed using kinematics assumption.
56  // The lower limit is tuneable because hadronization might have issues at low W (as in PYTHIA6).
57  // To be consistent the cross section must be computed in the W range where the events are generated.
58  const Kinematics & kinematics = interaction -> Kine();
59  const KPhaseSpace & ps = interaction -> PhaseSpace();
60  double W = kinematics.W();
61  Range1D_t Wl = ps.WLim();
62  Wl.min = TMath::Max(Wl.min,fWmin);
63  if ( W<Wl.min ) return 0.;
64  else if ( W>Wl.max ) return 0.;
65 
66  const InitialState & init_state = interaction -> InitState();
67 
68  double y = kinematics.y();
69  double Q2 = kinematics.Q2();
70  double x = kinematics.x();
71  double E = init_state.ProbeE(kRfLab);
72  double Mnuc = init_state.Tgt().HitNucMass();
73  double Mlep2 = TMath::Power(interaction->FSPrimLepton()->Mass(),2);
74 
75  // Get F1,F2,F3 for particular quark channel and compute differential xsec
76  SF_xQ2 sf = sf_tbl->EvalQrkSFLO( interaction, x, Q2 );
77  double xsec = (fMassTerms) ? ds_dxdy_mass( sf, x, y, E, Mnuc, Mlep2 ) : ds_dxdy( sf, x, y );
78 
79  // If NLO is enable we compute sigma_NLO/sigma_LO. Then the quark xsec
80  // is multiplied by this ratio.
81  // This is done because at NLO we can only compute the nucleon xsec. But
82  // for the hadronization we need the different quark contributions.
83  // This could be avoid if a NLO parton showering is introduced.
84  if (fSFinfo.IsNLO && xsec>0.) {
85  SF_xQ2 sflo = sf_tbl->EvalNucSFLO(interaction,x,Q2);
86  SF_xQ2 sfnlo = sf_tbl->EvalNucSFNLO(interaction,x,Q2);
87  double lo = (fMassTerms) ? ds_dxdy_mass( sflo, x, y, E, Mnuc, Mlep2 ) : ds_dxdy( sflo, x, y );
88  if (lo>0.) {
89  double nlo = (fMassTerms) ? ds_dxdy_mass( sfnlo, x, y, E, Mnuc, Mlep2 ) : ds_dxdy( sfnlo, x, y );
90  xsec *= nlo / lo;
91  }
92  }
93 
94  // Compute the front factor
95  double propagator = 0;
96  if (interaction -> ProcInfo().IsWeakCC()) propagator = TMath::Power( fSFinfo.MassW*fSFinfo.MassW/(Q2+fSFinfo.MassW*fSFinfo.MassW), 2);
97  else propagator = TMath::Power( fSFinfo.MassZ*fSFinfo.MassZ/(Q2+fSFinfo.MassZ*fSFinfo.MassZ)/(1.-fSFinfo.Rho), 2);
98 
99  xsec *= kGF2/(2*kPi*x) * propagator;
100 
101  LOG("HEDISPXSec", pINFO) << "d2xsec/dxdy[FreeN] (x= " << x << ", y= " << y << ", Q2= " << Q2 << ") = " << xsec;
102 
103  // The algorithm computes d^2xsec/dxdy. Check whether variable tranformation is needed
104  if( kps!=kPSxQ2fE ) xsec *= utils::kinematics::Jacobian(interaction,kPSxQ2fE,kps);
105 
106  // If requested return the free nucleon xsec even for input nuclear tgt
107  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
108 
109  // Compute nuclear cross section (simple scaling here, corrections must have been included in the structure functions)
110  int NNucl = (pdg::IsProton(init_state.Tgt().HitNucPdg())) ? init_state.Tgt().Z() : init_state.Tgt().N();
111  xsec *= NNucl;
112 
113  return xsec;
114 
115 }
116 //____________________________________________________________________________
117 double HEDISPXSec::ds_dxdy(SF_xQ2 sf, double x, double y ) const
118 {
119 
120  // We neglect F4 and F5 and higher order terms.
121  double term1 = y * ( x*y );
122  double term2 = ( 1 - y );
123  double term3 = ( x*y*(1-y/2) );
124 
125  LOG("HEDISPXSec", pDEBUG) << sf.F1 << " " << sf.F2 << " " << sf.F3;
126  LOG("HEDISPXSec", pDEBUG) << term1*sf.F1 + term2*sf.F2 + term3*sf.F3;
127 
128  return fmax( term1*sf.F1 + term2*sf.F2 + term3*sf.F3 , 0.);
129 
130 }
131 //____________________________________________________________________________
132 double HEDISPXSec::ds_dxdy_mass(SF_xQ2 sf, double x, double y, double e, double mt, double ml2 ) const
133 {
134 
135  double term1 = y * ( x*y + ml2/2/e/mt );
136  double term2 = ( 1 - y - mt*x*y/2/e - ml2/4/e/e );
137  double term3 = (x*y*(1-y/2) - y*ml2/4/mt/e);
138  double term4 = x*y*ml2/2/mt/e + ml2*ml2/4/mt/mt/e/e;
139  double term5 = -1.*ml2/2/mt/e;
140 
141  double F4 = 0.;
142  double F5 = sf.F2/x;
143 
144  LOG("HEDISPXSec", pDEBUG) << sf.F1 << " " << sf.F2 << " " << sf.F3;
145  LOG("HEDISPXSec", pDEBUG) << term1*sf.F1 + term2*sf.F2 + term3*sf.F3;
146 
147  return fmax( term1*sf.F1 + term2*sf.F2 + term3*sf.F3 + term4*F4 + term5*F5 , 0.);
148 
149 }
150 //____________________________________________________________________________
151 double HEDISPXSec::Integral(const Interaction * interaction) const
152 {
153 
154  return fXSecIntegrator->Integrate(this,interaction);
155 
156 }
157 //____________________________________________________________________________
158 bool HEDISPXSec::ValidProcess(const Interaction * interaction) const
159 {
160 
161  if(interaction->TestBit(kISkipProcessChk)) return true;
162 
163  const ProcessInfo & proc_info = interaction->ProcInfo();
164  if(!proc_info.IsDeepInelastic()) return false;
165 
166  const InitialState & init_state = interaction -> InitState();
167  int probe_pdg = init_state.ProbePdg();
168  if(!pdg::IsLepton(probe_pdg)) return false;
169 
170  if(! init_state.Tgt().HitNucIsSet()) return false;
171 
172  int hitnuc_pdg = init_state.Tgt().HitNucPdg();
173  if(!pdg::IsNeutronOrProton(hitnuc_pdg)) return false;
174 
175  return true;
176 }
177 //____________________________________________________________________________
178 void HEDISPXSec::Configure(const Registry & config)
179 {
180  Algorithm::Configure(config);
181  this->LoadConfig();
182 }
183 //____________________________________________________________________________
184 void HEDISPXSec::Configure(string config)
185 {
186  Algorithm::Configure(config);
187  this->LoadConfig();
188 }
189 //____________________________________________________________________________
191 {
192 
193  //-- load the differential cross section integrator
194  fXSecIntegrator = dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
195  assert(fXSecIntegrator);
196 
197  // Minimum value of W (typically driven by hadronization limitation)
198  GetParam("Xsec-Wmin", fWmin);
199  GetParam("Mass-Terms", fMassTerms);
200 
201  // Information about Structure Functions
202  GetParam("LHAPDF-set", fSFinfo.LHAPDFset );
203  GetParam("LHAPDF-member", fSFinfo.LHAPDFmember );
204  GetParam("Is-NLO", fSFinfo.IsNLO );
205  GetParam("Scheme", fSFinfo.Scheme );
206  GetParam("Quark-Threshold", fSFinfo.QrkThrs );
207  GetParam("NGridX", fSFinfo.NGridX );
208  GetParam("NGridQ2", fSFinfo.NGridQ2 );
209  GetParam("XGrid-Min", fSFinfo.XGridMin );
210  GetParam("Q2Grid-Min", fSFinfo.Q2GridMin );
211  GetParam("Q2Grid-Max", fSFinfo.Q2GridMax );
212  GetParam("MassW", fSFinfo.MassW );
213  GetParam("MassZ", fSFinfo.MassZ );
214  GetParam("Rho", fSFinfo.Rho );
215  GetParam("Sin2ThW", fSFinfo.Sin2ThW );
216  GetParam("Mud", fSFinfo.Vud );
217  GetParam("Mus", fSFinfo.Vus );
218  GetParam("Mub", fSFinfo.Vub );
219  GetParam("Mcd", fSFinfo.Vcd );
220  GetParam("Mcs", fSFinfo.Vcs );
221  GetParam("Mcb", fSFinfo.Vcb );
222  GetParam("Mtd", fSFinfo.Vtd );
223  GetParam("Mts", fSFinfo.Vts );
224  GetParam("Mtb", fSFinfo.Vtb );
225 
226 }
Cross Section Calculation Interface.
virtual ~HEDISPXSec()
Definition: HEDISPXSec.cxx:41
double W(bool selected=false) const
Definition: Kinematics.cxx:157
double ds_dxdy(SF_xQ2 sf, double x, double y) const
Definition: HEDISPXSec.cxx:117
double Integral(const Interaction *i) const
Definition: HEDISPXSec.cxx:151
void Configure(const Registry &config)
Definition: HEDISPXSec.cxx:178
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
int HitNucPdg(void) const
Definition: Target.cxx:304
const XSecIntegratorI * fXSecIntegrator
diff. xsec integrator
Definition: HEDISPXSec.h:52
A simple [min,max] interval for doubles.
Definition: Range1.h:42
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
Definition: HEDISPXSec.cxx:46
double HitNucMass(void) const
Definition: Target.cxx:233
SF_xQ2 EvalNucSFLO(const Interaction *in, double x, double Q2)
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double x(bool selected=false) const
Definition: Kinematics.cxx:99
enum genie::EKinePhaseSpace KinePhaseSpace_t
double y(bool selected=false) const
Definition: Kinematics.cxx:112
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
static HEDISStrucFunc * Instance(SF_info sfinfo)
Summary information for an interaction.
Definition: Interaction.h:56
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Definition: HEDISPXSec.cxx:158
void LoadConfig(void)
Definition: HEDISPXSec.cxx:190
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
const double e
double fWmin
Minimum value of W.
Definition: HEDISPXSec.h:54
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
double ds_dxdy_mass(SF_xQ2 sf, double x, double y, double e, double mt, double ml2) const
Definition: HEDISPXSec.cxx:132
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
int Z(void) const
Definition: Target.h:68
#define pINFO
Definition: Messenger.h:62
static constexpr double ps
Definition: Units.h:99
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
double max
Definition: Range1.h:53
int N(void) const
Definition: Target.h:69
bool HitNucIsSet(void) const
Definition: Target.cxx:283
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:351
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
SF_xQ2 EvalQrkSFLO(const Interaction *in, double x, double Q2)
double min
Definition: Range1.h:52
SF_info fSFinfo
Information used to computed SF.
Definition: HEDISPXSec.h:56
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
bool fMassTerms
Account for second order effects in DDxsec.
Definition: HEDISPXSec.h:55
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double ProbeE(RefFrame_t rf) const
SF_xQ2 EvalNucSFNLO(const Interaction *in, double x, double Q2)
bool IsLepton(int pdgc)
Definition: PDGUtils.cxx:86
Range1D_t WLim(void) const
W limits.
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345