include/apfel/mkpdf.h

0001 //
0002 // APFEL++ 2017
0003 //
0004 // Author: Valerio Bertone: valerio.bertone@cern.ch
0005 //
0006
0007 #pragma once
0008
0009 #include "apfel/config.h"
0010
0011 #ifdef WITH_LHAPDF
0012
0013 #include "apfel/initialiseevolution.h"
0014
0015 #include <LHAPDF/LHAPDF.h>
0016 #include <LHAPDF/GridPDF.h>
0017
0018 /**
0019  * @brief Function that constructs a pointer to an LHAPDF::PDF object
0020  * using an apfel::InitialiseEvolution object as an input in the place
0021  * of a standard LHAPDF grid.
0022  * @param ev: APFEL++ InitialiseEvolution object
0023  * @return a pointer to an LHAPDF::PDF object
0024  */
0025 LHAPDF::PDF* mkPDF(apfel::InitialiseEvolution const& ev)
0026 {
0027   // Get knot array from APFEL
0028   const std::map<double, std::map<int, apfel::LHKnotArray>> ka = ev.KnotArray();
0029
0030   // Knot array to be fed to LHAPDF
0031   LHAPDF::KnotArray data;
0032
0033   // Loop over Q subgrids to accumulate x-grid, q2-grid, and flavour
0034   // IDs.
0035   for (auto const& sg : ka)
0036     {
0037       // Get x-grid from first flavour of the first subgrid
0038       if (data.xs().empty())
0039         data.setxknots() = sg.second.begin()->second.xs;
0040
0041       // Accumulate q2 grid using the first flavour
0042       data.setq2knots().insert(data.setq2knots().end(), sg.second.begin()->second.q2s.begin(), sg.second.begin()->second.q2s.end());
0043
0044       // Fill in the flavour ID vector (use 21 for the gluon)
0045       if (data.setPids().empty())
0046         for (auto const& id : sg.second)
0047           data.setPids().push_back((id.first == 0 ? 21 : id.first));
0048     }
0049
0050   // Set up the knots of the Knotarray
0051   data.setShape() = std::vector<size_t> {data.xs().size(), data.q2s().size(), data.setPids().size()};
0052   data.fillLogKnots();
0053   data.initPidLookup();
0054
0055   // Set size of data vector
0056   data.setGrid().resize(data.shape(0) * data.shape(1) * data.shape(2), 0.);
0057
0058   // Loop over Q subgrids to accumulate data
0059   int nqprev = 0;
0060   for (auto const& sg : ka)
0061     {
0062       // Initialise flavour index
0063       int ifl = 0;
0064
0065       // Get q2-subgrid size
0066       const int q2size = (int) sg.second.begin()->second.xfs.size() / data.shape(0);
0067
0068       // Loop over flavours
0069       for (auto const& id : sg.second)
0070         {
0071           // Get vector of distributions
0072           const std::vector<double> f = id.second.xfs;
0073
0074           // Loop over x-grid
0075           for (int ix = 0; ix < (int) data.shape(0); ix++)
0076             {
0077               int iq = nqprev;
0078               // Loop over q2-sugrid
0079               for (int iqs = 0; iqs < q2size; iqs++)
0080                 data.setGrid()[ix * data.shape(2) * data.shape(1) + iq++ * data.shape(2) + ifl] = f[ix * q2size + iqs];
0081             }
0082           ifl++;
0083         }
0084       nqprev += q2size;
0085     }
0086
0087   // Fill in alpha_s vector
0088   std::vector<double> als;
0089   for (auto const& q2 : data.q2s())
0090     als.push_back(ev.Alphas(sqrt(q2)));
0091
0092   // Construct alpha_s object
0093   LHAPDF::AlphaS_Ipol* as = new LHAPDF::AlphaS_Ipol{};
0094   as->setQ2Values(data.q2s());
0095   as->setAlphaSValues(als);
0096
0097   // Define an LHAPDF GridPDF object
0098   LHAPDF::GridPDF* dist = new LHAPDF::GridPDF{};
0099
0100   // Pass knot arrays to LHAPDF
0101   dist->Data() = data;
0102
0103   // Set interpolator
0104   dist->setInterpolator((std::string) "logcubic");
0105
0106   // Set extrapolator
0107   dist->setExtrapolator((std::string) "continuation");
0108
0109   // Set flavours
0110   dist->setFlavors(data.setPids());
0111
0112   // Set alpha_s
0113   dist->setAlphaS(as);
0114
0115   // Set scale bounds
0116   dist->info().set_entry("QMin", sqrt(data.q2s().front()));
0117   dist->info().set_entry("QMax", sqrt(data.q2s().back()));
0118
0119   // Set x bounds
0120   dist->info().set_entry("XMin", data.xs().front());
0121   dist->info().set_entry("XMax", data.xs().back());
0122
0123   // Set quark masses and thresholds
0124   const std::vector<std::string> Qnames{"Down", "Up", "Strange", "Charm", "Bottom", "Top"};
0125   const std::vector<double> trhs = ev.GetEvolutionSetup().Thresholds;
0126   for (int iq = 0; iq < (int) Qnames.size(); iq++)
0127     {
0128       dist->info().set_entry("M" + Qnames[iq], (iq < (int) trhs.size() ? trhs[iq] : 1e8 + iq));
0129       dist->info().set_entry("Threshold" + Qnames[iq], (iq < (int) trhs.size() ? trhs[iq] : 1e8 + iq));
0130     }
0131
0132   // Return object
0133   return dist;
0134 }
0135
0136 #endif