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