GRV98LO.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3 
4  Copyright (c) 2003-2019, The GENIE Collaboration
5  For the full text of the license visit http://copyright.genie-mc.org
6  or see $GENIE/LICENSE
7 
8  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
9  University of Liverpool STFC, Rutherford Appleton Laboratory
10 
11  For the class documentation see the corresponding header file.
12 
13  Important revisions:
14 
15 */
16 //____________________________________________________________________________
17 
18 #include <sstream>
19 #include <fstream>
20 #include <cstdlib>
21 
22 #include <TSystem.h>
23 #include <TMath.h>
24 
27 
28 using namespace std;
29 using namespace genie;
30 
31 //____________________________________________________________________________
32 GRV98LO::GRV98LO() :
33 PDFModelI("genie::GRV98LO"),
34  fXUVF(NULL),
35  fXDVF(NULL),
36  fXDEF(NULL),
37  fXUDF(NULL),
38  fXSF (NULL),
39  fXGF (NULL)
40 {
41  this->Initialize();
42 }
43 //____________________________________________________________________________
45 PDFModelI("genie::GRV98LO", config),
46  fXUVF(NULL),
47  fXDVF(NULL),
48  fXDEF(NULL),
49  fXUDF(NULL),
50  fXSF (NULL),
51  fXGF (NULL)
52 {
53  LOG("GRV98LO", pDEBUG) << "GRV98LO configuration:\n " << GetConfig() ;
54 
55  this->Initialize();
56 }
57 //____________________________________________________________________________
59 {
60  if (fXUVF) {delete fXUVF; fXUVF = NULL;}
61  if (fXDVF) {delete fXDVF; fXDVF = NULL;}
62  if (fXDEF) {delete fXDEF; fXDEF = NULL;}
63  if (fXUDF) {delete fXUDF; fXUDF = NULL;}
64  if (fXSF ) {delete fXSF ; fXSF = NULL;}
65  if (fXGF ) {delete fXGF ; fXGF = NULL;}
66 }
67 //____________________________________________________________________________
68 double GRV98LO::UpValence(double x, double Q2) const
69 {
70  return AllPDFs(x,Q2).uval;
71 }
72 //____________________________________________________________________________
73 double GRV98LO::DownValence(double x, double Q2) const
74 {
75  return AllPDFs(x,Q2).dval;
76 }
77 //____________________________________________________________________________
78 double GRV98LO::UpSea(double x, double Q2) const
79 {
80  return AllPDFs(x,Q2).usea;
81 }
82 //____________________________________________________________________________
83 double GRV98LO::DownSea(double x, double Q2) const
84 {
85  return AllPDFs(x,Q2).dsea;
86 }
87 //____________________________________________________________________________
88 double GRV98LO::Strange(double x, double Q2) const
89 {
90  return AllPDFs(x,Q2).str;
91 }
92 //____________________________________________________________________________
93 double GRV98LO::Charm(double x, double Q2) const
94 {
95  return AllPDFs(x,Q2).chm;
96 }
97 //____________________________________________________________________________
98 double GRV98LO::Bottom(double x, double Q2) const
99 {
100  return AllPDFs(x,Q2).bot;
101 }
102 //____________________________________________________________________________
103 double GRV98LO::Top(double x, double Q2) const
104 {
105  return AllPDFs(x,Q2).top;
106 }
107 //____________________________________________________________________________
108 double GRV98LO::Gluon(double x, double Q2) const
109 {
110  return AllPDFs(x,Q2).gl;
111 }
112 //____________________________________________________________________________
113 PDF_t GRV98LO::AllPDFs(double x, double Q2) const
114 {
115  PDF_t pdf;
116 
117  if(!fInitialized) {
118  LOG("GRV98LO", pWARN)
119  << "GRV98LO algorithm was not initialized succesfully";
120  pdf.uval = 0.;
121  pdf.dval = 0.;
122  pdf.usea = 0.;
123  pdf.dsea = 0.;
124  pdf.str = 0.;
125  pdf.chm = 0.;
126  pdf.bot = 0.;
127  pdf.top = 0.;
128  pdf.gl = 0.;
129  return pdf;
130  }
131 
132  LOG("GRV98LO", pDEBUG)
133  << "Inputs x = " << x << ", Q2 = " << Q2;
134 
135  // apply kinematical limits
136 //Q2 = TMath::Max(Q2, fGridQ2[0]);
137  if(Q2 <= 0.8) Q2 = 0.80001;
138  Q2 = std::min(Q2, fGridQ2[kNQ2-1]);
139  x = std::max(x, fGridXbj[0]);
140  x = std::min(x, fGridXbj[kNXbj-1]);
141 
142  double logx = std::log(x);
143  double logQ2 = std::log(Q2);
144  double x1 = 1-x;
145  double xv = std::sqrt(x);
146  double xs = std::pow(x, -0.2);
147  double x1p3 = x1*x1*x1;
148  double x1p4 = x1*x1p3;
149  double x1p5 = x1*x1p4;
150  double x1p7 = x1p3*x1p4;
151 
152  double uv = fXUVF->Eval(logx,logQ2) * x1p3 * xv;
153  double dv = fXDVF->Eval(logx,logQ2) * x1p4 * xv;
154  double de = fXDEF->Eval(logx,logQ2) * x1p7 * xv;
155  double ud = fXUDF->Eval(logx,logQ2) * x1p7 * xs;
156  double us = 0.5 * (ud - de);
157  double ds = 0.5 * (ud + de);
158  double ss = fXSF->Eval(logx,logQ2) * x1p7 * xs;
159  double gl = fXGF->Eval(logx,logQ2) * x1p5 * xs;
160 
161  pdf.uval = uv;
162  pdf.dval = dv;
163  pdf.usea = us;
164  pdf.dsea = ds;
165  pdf.str = ss;
166  pdf.chm = 0.;
167  pdf.bot = 0.;
168  pdf.top = 0.;
169  pdf.gl = gl;
170 
171  return pdf;
172 }
173 //____________________________________________________________________________
175 {
176  Algorithm::Configure(config);
177  this->Initialize();
178 }
179 //____________________________________________________________________________
181 {
182  Algorithm::Configure(config);
183  this->Initialize();
184 }
185 //____________________________________________________________________________
187 {
188  fInitialized = false;
189 
190  const char * genie_dir = gSystem->Getenv("GENIE");
191  if(!genie_dir) return;
192 
193  string grid_file_name =
194  string(gSystem->Getenv("GENIE")) + string("/data/evgen/pdfs/GRV98lo_patched.LHgrid");
195 
196  LOG("GRV98LO", pNOTICE)
197  << "Reading grid file from:\n " << grid_file_name;
198 
199  ifstream grid_file;
200  grid_file.open (grid_file_name.c_str());
201 
202  char rubbish[1000];
203 
204  const int nskip = 21;
205  for(int i=0; i<nskip; i++) {
206  grid_file.getline(rubbish,1000);
207  LOG("GRV98LO", pDEBUG) << "Skipping: " << rubbish;
208  }
209 
210  // x's
211  //
212 
213  LOG("GRV98LO", pDEBUG) << "Reading x_bj grid values";
214  for(int j=0; j < kNXbj; j++) {
215  double xbj = -1;
216  grid_file >> xbj;
217  // check against known limits
218  // ...
219  fGridXbj[j] = xbj;
220  fGridLogXbj[j] = TMath::Log(xbj);
221  }
222  ostringstream grid_values;
223  grid_values << "(";
224  for(int j=0; j < kNXbj; j++) {
225  grid_values << fGridXbj[j];
226  if(j == kNXbj - 1) { grid_values << ")"; }
227  else { grid_values << ", "; }
228  }
229  LOG("GRV98LO", pDEBUG)
230  << "x_bj grid values: " << grid_values.str();
231 
232  // Q^2s
233  //
234 
235  LOG("GRV98LO", pDEBUG) << "Reading Q^2 grid values.";
236  for(int i=0; i < kNQ2; i++) {
237  double Q2 = -1;
238  grid_file >> Q2;
239  // check against known limits
240  // ...
241  fGridQ2[i] = Q2;
242  fGridLogQ2[i] = TMath::Log(Q2);
243  }
244  grid_values.str("");
245  grid_values << "(";
246  for(int i=0; i < kNQ2; i++) {
247  grid_values << fGridQ2[i];
248  if(i == kNQ2 - 1) { grid_values << ")"; }
249  else { grid_values << ", "; }
250  }
251  LOG("GRV98LO", pDEBUG)
252  << "Q^2 grid values: " << grid_values.str() << "GeV^2";
253 
254  // skip again
255  grid_file.getline(rubbish,1000);
256  LOG("GRV98LO", pDEBUG) << "Skipping: " << rubbish;
257  grid_file.getline(rubbish,1000);
258  LOG("GRV98LO", pDEBUG) << "Skipping: " << rubbish;
259 
260  // pdf values on grid points
261  //
262 
263  LOG("GRV98LO", pDEBUG) << "Reading PDF values on grid points";
264 
265  int k=0;
266  for(int j=0; j < kNXbj-1; j++) {
267  for(int i=0; i < kNQ2; i++) {
268  double p0 = 0;
269  double p1 = 0;
270  double p2 = 0;
271  double p3 = 0;
272  double p4 = 0;
273  double p5 = 0;
274  grid_file >> p0 >> p1 >> p2 >> p3 >> p4 >> p5;
275  LOG("GRV98LO", pDEBUG)
276  << "Row: " << k << ", grid point: (" << i << ", " << j << ") ->"
277  << " p0 = " << p0
278  << ", p1 = " << p1
279  << ", p2 = " << p2
280  << ", p3 = " << p3
281  << ", p4 = " << p4
282  << ", p5 = " << p5;
283  fParton [0][i][j] = p0;
284  fParton [1][i][j] = p1;
285  fParton [2][i][j] = p2;
286  fParton [3][i][j] = p3;
287  fParton [4][i][j] = p4;
288  fParton [5][i][j] = p5;
289  k++;
290  }
291  }
292 
293  grid_file.close();
294 
295  // arrays for interpolation routines
296  //
297 
298  vector<double> gridLogQ2 (kNQ2);
299  vector<double> gridLogXbj(kNXbj);
300  vector<double> knotsXUVF(kNQ2*kNXbj);
301  vector<double> knotsXDVF(kNQ2*kNXbj);
302  vector<double> knotsXDEF(kNQ2*kNXbj);
303  vector<double> knotsXUDF(kNQ2*kNXbj);
304  vector<double> knotsXSF (kNQ2*kNXbj);
305  vector<double> knotsXGF (kNQ2*kNXbj);
306 
307  k=0;
308  for(int i=0; i < kNQ2; i++) {
309  double logQ2 = std::log(fGridQ2[i]);
310  gridLogQ2[i] = logQ2;
311  for(int j=0; j < kNXbj - 1; j++) {
312  double logx = std::log(fGridXbj[j]);
313  gridLogXbj[j] = logx;
314  double xb0v = std::sqrt(fGridXbj[j]);
315  double xb0s = std::pow(fGridXbj[j], -0.2);
316  double xb1 = 1 - fGridXbj[j];
317  double xb1p3 = std::pow(xb1, 3.);
318  double xb1p4 = std::pow(xb1, 4.);
319  double xb1p5 = std::pow(xb1, 5.);
320  double xb1p7 = std::pow(xb1, 7.);
321  knotsXUVF[k] = fParton[0][i][j] / (xb1p3 * xb0v);
322  knotsXDVF[k] = fParton[1][i][j] / (xb1p4 * xb0v);
323  knotsXDEF[k] = fParton[2][i][j] / (xb1p7 * xb0v);
324  knotsXUDF[k] = fParton[3][i][j] / (xb1p7 * xb0s);
325  knotsXSF [k] = fParton[4][i][j] / (xb1p7 * xb0s);
326  knotsXGF [k] = fParton[5][i][j] / (xb1p5 * xb0s);
327  k++;
328  }
329  double logxmax = TMath::Log(fGridXbj[kNXbj-1]);
330  gridLogXbj[kNXbj-1] = logxmax;
331  knotsXUVF[k] = 0;
332  knotsXDVF[k] = 0;
333  knotsXDEF[k] = 0;
334  knotsXUDF[k] = 0;
335  knotsXSF [k] = 0;
336  knotsXGF [k] = 0;
337  k++;
338  }
339 
340  fXUVF = new Interpolator2D(gridLogXbj.size(),&gridLogXbj[0],gridLogQ2.size(),&gridLogQ2[0],&knotsXUVF[0]);
341  fXDVF = new Interpolator2D(gridLogXbj.size(),&gridLogXbj[0],gridLogQ2.size(),&gridLogQ2[0],&knotsXDVF[0]);
342  fXDEF = new Interpolator2D(gridLogXbj.size(),&gridLogXbj[0],gridLogQ2.size(),&gridLogQ2[0],&knotsXDEF[0]);
343  fXUDF = new Interpolator2D(gridLogXbj.size(),&gridLogXbj[0],gridLogQ2.size(),&gridLogQ2[0],&knotsXUDF[0]);
344  fXSF = new Interpolator2D(gridLogXbj.size(),&gridLogXbj[0],gridLogQ2.size(),&gridLogQ2[0],&knotsXSF [0]);
345  fXGF = new Interpolator2D(gridLogXbj.size(),&gridLogXbj[0],gridLogQ2.size(),&gridLogQ2[0],&knotsXGF [0]);
346 
347  fInitialized = true;
348 }
349 //____________________________________________________________________________
T max(const caf::Proxy< T > &a, T b)
double Top(double x, double Q2) const
Definition: GRV98LO.cxx:103
double Charm(double x, double Q2) const
Definition: GRV98LO.cxx:93
double fGridLogQ2[kNQ2]
Definition: GRV98LO.h:81
PDF_t AllPDFs(double x, double Q2) const
Definition: GRV98LO.cxx:113
Interpolator2D * fXDVF
Definition: GRV98LO.h:89
double Bottom(double x, double Q2) const
Definition: GRV98LO.cxx:98
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
double Eval(const double &x, const double &y) const
void Initialize(void)
Definition: GRV98LO.cxx:186
double DownSea(double x, double Q2) const
Definition: GRV98LO.cxx:83
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
Float_t ss
Definition: plot.C:24
Definition: config.py:1
double fGridLogXbj[kNXbj]
Definition: GRV98LO.h:83
Interpolator2D * fXUVF
Definition: GRV98LO.h:88
Pure abstract base class. Defines the PDFModelI interface to be implemented by wrapper classes to exi...
Definition: PDFModelI.h:29
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
void Configure(const Registry &config)
Definition: GRV98LO.cxx:174
double fGridQ2[kNQ2]
Definition: GRV98LO.h:80
Interpolator2D * fXGF
Definition: GRV98LO.h:93
double DownValence(double x, double Q2) const
Definition: GRV98LO.cxx:73
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
static const int kNXbj
Definition: GRV98LO.h:73
Interpolator2D * fXDEF
Definition: GRV98LO.h:90
Interpolator2D * fXUDF
Definition: GRV98LO.h:91
double Strange(double x, double Q2) const
Definition: GRV98LO.cxx:88
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
const double j
Definition: BetheBloch.cxx:29
double fGridXbj[kNXbj]
Definition: GRV98LO.h:82
bool fInitialized
Definition: GRV98LO.h:68
A very simple service to remember what detector we&#39;re working in.
double UpValence(double x, double Q2) const
Definition: GRV98LO.cxx:68
#define pWARN
Definition: Messenger.h:61
bool logx
static constexpr double us
A struct to hold PDF set data.
Interpolator2D * fXSF
Definition: GRV98LO.h:92
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
A 2D interpolator using the GSL spline type If GSL version is not sufficient, does an inefficient ver...
string genie_dir(std::getenv("GENIE"))
#define pNOTICE
Definition: Messenger.h:62
double fParton[kNParton][kNQ2][kNXbj-1]
Definition: GRV98LO.h:84
double UpSea(double x, double Q2) const
Definition: GRV98LO.cxx:78
T min(const caf::Proxy< T > &a, T b)
static const int kNQ2
Definition: GRV98LO.h:72
double Gluon(double x, double Q2) const
Definition: GRV98LO.cxx:108
#define pDEBUG
Definition: Messenger.h:64
enum BeamMode string