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