/examples/CCTest2.cc

An example of a program using the C++ interface to LHAPDF to calculate PDF errors.

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

Generated on 8 Feb 2016 for LHAPDF C++ wrapper by  doxygen 1.4.7