GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BardinIMDRadCorPXSec.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 #include <Math/Integrator.h>
13 
15 #include "Framework/Conventions/GBuild.h"
22 
23 using namespace genie;
24 using namespace genie::constants;
25 
26 //____________________________________________________________________________
28 XSecAlgorithmI("genie::BardinIMDRadCorPXSec")
29 {
30 
31 }
32 //____________________________________________________________________________
34 XSecAlgorithmI("genie::BardinIMDRadCorPXSec", config)
35 {
36 
37 }
38 //____________________________________________________________________________
40 {
41 
42 }
43 //____________________________________________________________________________
45  const Interaction * interaction, KinePhaseSpace_t kps) const
46 {
47  if(! this -> ValidProcess (interaction) ) return 0.;
48  if(! this -> ValidKinematics (interaction) ) return 0.;
49 
50  const InitialState & init_state = interaction -> InitState();
51 
52  double E = init_state.ProbeE(kRfLab);
53  double sig0 = kGF2 * kElectronMass * E / kPi;
54  double re = 0.5 * kElectronMass / E;
55  double r = (kMuonMass2 / kElectronMass2) * re;
56  double y = interaction->Kine().y();
57 
58  y = 1-y; //Note: y = (Ev-El)/Ev but in Bardin's paper y=El/Ev.
59 
60  double ymin = r + re;
61  double ymax = 1 + re + r*re / (1+re);
62 
63  double e = 1E-5;
64  ymax = TMath::Min(ymax,1-e); // avoid ymax=1, due to a log(1-y)
65 
66 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
67  LOG("BardinIMD", pDEBUG)
68  << "sig0 = " << sig0 << ", r = " << r << ", re = " << re;
69  LOG("BardinIMD", pDEBUG)
70  << "allowed y: [" << ymin << ", " << ymax << "]";
71 #endif
72 
73  if(y<ymin || y>ymax) return 0;
74 
75  double xsec = 2 * sig0 * ( 1 - r + (kAem/kPi) * Fa(re,r,y) );
76 
77 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
78  LOG("BardinIMD", pINFO)
79  << "dxsec[1-loop]/dy (Ev = " << E << ", y = " << y << ") = " << xsec;
80 #endif
81 
82  // The algorithm computes dxsec/dy
83  // Check whether variable tranformation is needed
84  if(kps!=kPSyfE) {
85  double J = utils::kinematics::Jacobian(interaction,kPSyfE,kps);
86  xsec *= J;
87  }
88 
89  // If requested return the free electron xsec even for nuclear target
90  if( interaction->TestBit(kIAssumeFreeElectron) ) return xsec;
91 
92  // Scale for the number of scattering centers at the target
93  int Ne = init_state.Tgt().Z(); // num of scattering centers
94  xsec *= Ne;
95 
96  return xsec;
97 }
98 //____________________________________________________________________________
99 double BardinIMDRadCorPXSec::Integral(const Interaction * interaction) const
100 {
101  double xsec = fXSecIntegrator->Integrate(this,interaction);
102  return xsec;
103 }
104 //____________________________________________________________________________
105 bool BardinIMDRadCorPXSec::ValidProcess(const Interaction * interaction) const
106 {
107  if(interaction->TestBit(kISkipProcessChk)) return true;
108  return true;
109 }
110 //____________________________________________________________________________
111 double BardinIMDRadCorPXSec::Fa(double re, double r, double y) const
112 {
113  double y2 = y * y;
114  double rre = r * re;
115  double r_y = r/y;
116  double y_r = y/r;
117 
118  double fa = 0;
119 
120  fa = (1-r) * ( TMath::Log(y2/rre) * TMath::Log(1-r_y) +
121  TMath::Log(y_r) * TMath::Log(1-y) -
122  this->Li2( r ) +
123  this->Li2( y ) +
124  this->Li2( (r-y) / (1-y) ) +
125  1.5 * (1-r) * TMath::Log(1-r)
126  )
127  +
128 
129  0.5*(1+3*r) * ( this->Li2( (1-r_y) / (1-r) ) -
130  this->Li2( (y-r) / (1-r) ) -
131  TMath::Log(y_r) * TMath::Log( (y-r) / (1-r) )
132  )
133  +
134 
135  this->P(1,r,y) -
136  this->P(2,r,y) * TMath::Log(r) -
137  this->P(3,r,y) * TMath::Log(re) +
138  this->P(4,r,y) * TMath::Log(y) +
139  this->P(5,r,y) * TMath::Log(1-y) +
140  this->P(6,r,y) * (1 - r_y) * TMath::Log(1-r_y);
141 
142  return fa;
143 }
144 //____________________________________________________________________________
145 double BardinIMDRadCorPXSec::P(int i, double r, double y) const
146 {
147  int kmin = -3;
148  int kmax = 2;
149  double p = 0;
150  for(int k = kmin; k <= kmax; k++) {
151  double c = this->C(i,k,r);
152  double yk = TMath::Power(y,k);
153  p += (c*yk);
154  }
155  return p;
156 }
157 //____________________________________________________________________________
158 double BardinIMDRadCorPXSec::Li2(double z) const
159 {
160  double epsilon = 1e-2;
161  double tmin = epsilon;
162  double tmax = 1. - epsilon;
163 
164 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
165  LOG("BardinIMD", pDEBUG)
166  << "Summing BardinIMDRadCorIntegrand in [" << tmin<< ", " << tmax<< "]";
167 #endif
168 
169  ROOT::Math::IBaseFunctionOneDim * integrand = new
171  ROOT::Math::IntegrationOneDim::Type ig_type =
173 
174  double abstol = 1; // We mostly care about relative tolerance
175  double reltol = 1E-4;
176  int nmaxeval = 100000;
177  ROOT::Math::Integrator ig(*integrand,ig_type,abstol,reltol,nmaxeval);
178  double li2 = ig.Integral(tmin, tmax);
179 
180 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
181  LOG("BardinIMD", pDEBUG) << "Li2(z = " << z << ")" << li2;
182 #endif
183 
184  delete integrand;
185 
186  return li2;
187 }
188 //____________________________________________________________________________
189 double BardinIMDRadCorPXSec::C(int i, int k, double r) const
190 {
191  if ( i == 1 ) {
192 
193  if (k == -3) return -0.19444444*TMath::Power(r,3.);
194  else if (k == -2) return (0.083333333+0.29166667*r)*TMath::Power(r,2.);
195  else if (k == -1) return -0.58333333*r - 0.5*TMath::Power(r,2.) - TMath::Power(r,3.)/6.;
196  else if (k == 0) return -1.30555560 + 3.125*r + 0.375*TMath::Power(r,2.);
197  else if (k == 1) return -0.91666667 - 0.25*r;
198  else if (k == 2) return 0.041666667;
199  else return 0.;
200 
201  } else if ( i == 2 ) {
202 
203  if (k == -3) return 0.;
204  else if (k == -2) return 0.5*TMath::Power(r,2.);
205  else if (k == -1) return 0.5*r - 2*TMath::Power(r,2.);
206  else if (k == 0) return 0.25 - 0.75*r + 1.5*TMath::Power(r,2);
207  else if (k == 1) return 0.5;
208  else if (k == 2) return 0.;
209  else return 0.;
210 
211  } else if ( i == 3 ) {
212 
213  if (k == -3) return 0.16666667*TMath::Power(r,3.);
214  else if (k == -2) return 0.25*TMath::Power(r,2.)*(1-r);
215  else if (k == -1) return r-0.5*TMath::Power(r,2.);
216  else if (k == 0) return 0.66666667;
217  else if (k == 1) return 0.;
218  else if (k == 2) return 0.;
219  else return 0.;
220 
221  } else if ( i == 4 ) {
222 
223  if (k == -3) return 0.;
224  else if (k == -2) return TMath::Power(r,2.);
225  else if (k == -1) return r*(1-4.*r);
226  else if (k == 0) return 1.5*TMath::Power(r,2.);
227  else if (k == 1) return 1.;
228  else if (k == 2) return 0.;
229  else return 0.;
230 
231  } else if ( i == 5 ) {
232 
233  if (k == -3) return 0.16666667*TMath::Power(r,3.);
234  else if (k == -2) return -0.25*TMath::Power(r,2.)*(1+r);
235  else if (k == -1) return 0.5*r*(1+3*r);
236  else if (k == 0) return -1.9166667+2.25*r-1.5*TMath::Power(r,2);
237  else if (k == 1) return -0.5;
238  else if (k == 2) return 0.;
239  else return 0.;
240 
241  } else if ( i == 6 ) {
242 
243  if (k == -3) return 0.;
244  else if (k == -2) return 0.16666667*TMath::Power(r,2.);
245  else if (k == -1) return -0.25*r*(r+0.33333333);
246  else if (k == 0) return 1.25*(r+0.33333333);
247  else if (k == 1) return 0.5;
248  else if (k == 2) return 0.;
249  else return 0.;
250 
251  } else return 0.;
252 }
253 //____________________________________________________________________________
255 {
256  Algorithm::Configure(config);
257  this->LoadConfig();
258 }
259 //____________________________________________________________________________
260 void BardinIMDRadCorPXSec::Configure(string param_set)
261 {
262  Algorithm::Configure(param_set);
263  this->LoadConfig();
264 }
265 //____________________________________________________________________________
267 {
268  ////fIntegrator =
269 //// dynamic_cast<const IntegratorI *> (this->SubAlg("Integrator"));
270 ///// assert(fIntegrator);
271 
273  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
274  assert(fXSecIntegrator);
275 }
276 //____________________________________________________________________________
277 // Auxiliary scalar function for internal integration
278 //____________________________________________________________________________
280 ROOT::Math::IBaseFunctionOneDim()
281 {
282  fZ = z;
283 }
284 //____________________________________________________________________________
286 {
287 
288 }
289 //____________________________________________________________________________
291 {
292  return 1;
293 }
294 //____________________________________________________________________________
296 {
297  if(xin<=0) return 0.;
298  if(xin*fZ >= 1.) return 0.;
299  double f = TMath::Log(1.-fZ*xin)/xin;
300  return f;
301 }
302 //____________________________________________________________________________
303 ROOT::Math::IBaseFunctionOneDim *
305 {
307 }
308 //____________________________________________________________________________
Cross Section Calculation Interface.
void Configure(const Registry &config)
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition: GSLUtils.cxx:23
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
Cross Section Integrator Interface.
enum genie::EKinePhaseSpace KinePhaseSpace_t
static const double kElectronMass
const double epsilon
double y(bool selected=false) const
Definition: Kinematics.cxx:112
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?
double P(int i, double r, double y) const
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
double C(int i, int k, double r) const
static const double kElectronMass2
const Kinematics & Kine(void) const
Definition: Interaction.h:71
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int Z(void) const
Definition: Target.h:68
#define pINFO
Definition: Messenger.h:62
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
const XSecIntegratorI * fXSecIntegrator
differential x-sec integrator
double Fa(double re, double r, double y) const
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
const Target & Tgt(void) const
Definition: InitialState.h:66
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double ProbeE(RefFrame_t rf) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
double Integral(const Interaction *i) const
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
const UInt_t kIAssumeFreeElectron
Definition: Interaction.h:50
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345