#include "LHAPDF/LHAPDF.h"
#include <iostream>
#include <cmath>
#include <cstdio>
#include <cstdlib>
using namespace std;
using namespace LHAPDF;
double logdist_x(double xmin, double xmax, size_t ix, size_t nx) {
const double log10xmin = log10(xmin);
const double log10xmax = log10(xmax);
const double log10x = log10xmin + (ix/static_cast<double>(nx-1))*(log10xmax-log10xmin);
const double x = pow(10.0, log10x);
return x;
}
int main() {
const string NAME = "MRST2006nnlo";
initPDFSetM(1, NAME, LHGRID);
initPDFSetM(2, NAME, LHGRID);
initPDFSetM(3, NAME, LHGRID);
const int neigen = numberPDFM(1)/2;
cout << "Number of eigensets in this fit = " << neigen << endl;
cout << "Valid x-range = [" << xmin << ", " << xmax << "]" << endl;
const int nx = 10;
double q = 10;
int flav = 4;
initPDFM(1, 0);
vector<double> fc(nx), x(nx);
for (int ix = 0; ix < nx; ++ix) {
x[ix] = logdist_x(xmin, 0.9*xmax, ix, nx);
fc[ix] = xfxM(1, x[ix], q, flav);
}
vector<double> summax(nx), summin(nx), sum(nx);
#ifndef LHAPDF_LOWMEM
cout << "Using efficient set looping" << endl;
for (int ieigen = 1; ieigen <= neigen; ++ieigen) {
initPDFM(2, 2*ieigen-1);
initPDFM(3, 2*ieigen);
for (int ix = 0; ix < nx; ++ix) {
const double fp = xfxM(2, x[ix], q, flav);
const double fm = xfxM(3, x[ix], q, flav);
const double plus = max(max(fp-fc[ix], fm-fc[ix]),0.0);
const double minus = min(min(fp-fc[ix], fm-fc[ix]),0.0);
const double diff = fp-fm;
summax[ix] += plus*plus;
summin[ix] += minus*minus;
sum[ix] += diff*diff;
}
}
#else
cout << "Using low-mem mode set looping" << endl;
for (int ieigen = 1; ieigen <= neigen; ++ieigen) {
vector<double> fp(nx), fm(nx);
initPDFM(2, 2*ieigen-1);
for (int ix = 0; ix < nx; ++ix) {
fp[ix] = xfxM(2, x[ix], q, flav);
}
initPDFM(3, 2*ieigen);
for (int ix = 0; ix < nx; ++ix) {
fm[ix] = xfxM(3, x[ix], q, flav);
}
for (int ix = 0; ix < nx; ++ix) {
const double plus = max(max(fp[ix]-fc[ix], fm[ix]-fc[ix]), 0.0);
const double minus = min(min(fp[ix]-fc[ix], fm[ix]-fc[ix]), 0.0);
const double diff = fp[ix]-fm[ix];
summax[ix] += plus*plus;
summin[ix] += minus*minus;
sum[ix] += diff*diff;
}
}
#endif
cout << "flavour = " << flav << " Asymmetric (%) Symmetric (%)" << endl;
cout << " x Q**2 xf(x) plus minus +- " << endl;
for (int ix = 0; ix < nx; ++ix) {
printf("%0.7f %.0f %10.2E %8.2f %8.2f %8.2f \n",
x[ix], q*q, fc[ix],
sqrt(summax[ix])*100/fc[ix],
sqrt(summin[ix])*100/fc[ix],
0.5*sqrt(sum[ix])*100/fc[ix]);
}
return EXIT_SUCCESS;
}
#include "FortranWrappers.h"
#ifdef FC_DUMMY_MAIN
int FC_DUMMY_MAIN() { return 1; }
#endif