00001 00002 // Program to demonstrate usage of the MRST 2006 NNLO PDFs. // 00003 // to calculate errors. // 00005 #include "LHAPDF/LHAPDF.h" 00006 #include <iostream> 00007 #include <cmath> 00008 #include <cstdio> 00009 #include <cstdlib> 00010 using namespace std; 00011 00012 using namespace LHAPDF; 00013 00014 00015 double logdist_x(double xmin, double xmax, size_t ix, size_t nx) { 00016 const double log10xmin = log10(xmin); 00017 const double log10xmax = log10(xmax); 00018 const double log10x = log10xmin + (ix/static_cast<double>(nx-1))*(log10xmax-log10xmin); 00019 const double x = pow(10.0, log10x); 00020 return x; 00021 } 00022 00023 00024 int main() { 00025 // Show initialisation banner only once 00026 setVerbosity(LOWKEY); // or SILENT, for no banner at all 00027 00028 // You could explicitly set the path to the PDFsets directory 00029 // setPDFPath("/home/whalley/local/share/lhapdf/PDFsets"); 00030 00031 // Initialize PDF sets 00032 const string NAME = "MRST2006nnlo"; 00033 initPDFSetM(1, NAME, LHGRID); 00034 initPDFSetM(2, NAME, LHGRID); 00035 initPDFSetM(3, NAME, LHGRID); 00036 00037 // Find the number of eigensets from numberPDF() 00038 const int neigen = numberPDFM(1)/2; 00039 cout << "Number of eigensets in this fit = " << neigen << endl; 00040 // Find the min and max values of x and Q2 00041 const double xmin = getXmin(0); 00042 const double xmax = getXmax(0); 00043 cout << "Valid x-range = [" << xmin << ", " << xmax << "]" << endl; 00044 // Number of x values to sample 00045 const int nx = 10; 00046 // Set the Q scale and flavour 00047 double q = 10; 00048 int flav = 4; 00049 00050 // Get x's and central PDF values 00051 initPDFM(1, 0); 00052 vector<double> fc(nx), x(nx); 00053 for (int ix = 0; ix < nx; ++ix) { 00054 x[ix] = logdist_x(xmin, 0.9*xmax, ix, nx); 00055 fc[ix] = xfxM(1, x[ix], q, flav); 00056 } 00057 00058 // Sum over error contributions (two ways, depending on how LHDPAF was compiled) 00059 vector<double> summax(nx), summin(nx), sum(nx); 00060 #ifndef LHAPDF_LOWMEM 00061 // This is the normal, efficient, way to do this, with the error 00062 // sets being initialised the minimum number of times 00063 cout << "Using efficient set looping" << endl; 00064 for (int ieigen = 1; ieigen <= neigen; ++ieigen) { 00065 initPDFM(2, 2*ieigen-1); 00066 initPDFM(3, 2*ieigen); 00067 for (int ix = 0; ix < nx; ++ix) { 00068 // Find central and plus/minus values 00069 const double fp = xfxM(2, x[ix], q, flav); 00070 const double fm = xfxM(3, x[ix], q, flav); 00071 // Construct shifts 00072 const double plus = max(max(fp-fc[ix], fm-fc[ix]),0.0); 00073 const double minus = min(min(fp-fc[ix], fm-fc[ix]),0.0); 00074 const double diff = fp-fm; 00075 // Add it together 00076 summax[ix] += plus*plus; 00077 summin[ix] += minus*minus; 00078 sum[ix] += diff*diff; 00079 } 00080 } 00081 #else 00082 // In low memory mode, the sets need to be re-initialised with every 00083 // change of member. Using the approach above gives wrong answers, and 00084 // reinitialising in all the nested loops is sloooooow! The best way is 00085 // to calculate the values, plus and minus errors separately. 00086 cout << "Using low-mem mode set looping" << endl; 00087 for (int ieigen = 1; ieigen <= neigen; ++ieigen) { 00088 vector<double> fp(nx), fm(nx); 00089 initPDFM(2, 2*ieigen-1); 00090 for (int ix = 0; ix < nx; ++ix) { 00091 fp[ix] = xfxM(2, x[ix], q, flav); 00092 } 00093 initPDFM(3, 2*ieigen); 00094 for (int ix = 0; ix < nx; ++ix) { 00095 fm[ix] = xfxM(3, x[ix], q, flav); 00096 } 00097 for (int ix = 0; ix < nx; ++ix) { 00098 // Construct shifts 00099 const double plus = max(max(fp[ix]-fc[ix], fm[ix]-fc[ix]), 0.0); 00100 const double minus = min(min(fp[ix]-fc[ix], fm[ix]-fc[ix]), 0.0); 00101 const double diff = fp[ix]-fm[ix]; 00102 // Add it together 00103 summax[ix] += plus*plus; 00104 summin[ix] += minus*minus; 00105 sum[ix] += diff*diff; 00106 } 00107 } 00108 #endif 00109 00110 // Print out results 00111 cout << "flavour = " << flav << " Asymmetric (%) Symmetric (%)" << endl; 00112 cout << " x Q**2 xf(x) plus minus +- " << endl; 00113 for (int ix = 0; ix < nx; ++ix) { 00114 printf("%0.7f %.0f %10.2E %8.2f %8.2f %8.2f \n", 00115 x[ix], q*q, fc[ix], 00116 sqrt(summax[ix])*100/fc[ix], 00117 sqrt(summin[ix])*100/fc[ix], 00118 0.5*sqrt(sum[ix])*100/fc[ix]); 00119 } 00120 00121 return EXIT_SUCCESS; 00122 } 00123 00124 00125 00126 #include "FortranWrappers.h" 00127 #ifdef FC_DUMMY_MAIN 00128 int FC_DUMMY_MAIN() { return 1; } 00129 #endif