diff --git a/Makefile b/Makefile index 38f354d..36595d3 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ include Makefile.inc -CXXFLAGS = $(PRJCXXFLAGS) +CXXFLAGS = $(PRJCXXFLAGS) LDLIBS = ./src/libac_core.a $(PRJLDFLAGS) VPATH=./src diff --git a/Makefile.inc b/Makefile.inc index e74a3fe..e5798f5 100644 --- a/Makefile.inc +++ b/Makefile.inc @@ -10,17 +10,18 @@ SHELL=/bin/bash # APPLCOMB paths RESULTSDIR= -D RESULTS_PATH="./results/" +PINEAPPLGRID= -D PINEAPPL_PATH="../pineapplgrids/" APPLGRIDDIR= -D APPL_PATH="../applgrids/" DATABASEDIR= -D DB_PATH="./db/" -ALLDIR= $(RESULTSDIR) $(APPLGRIDDIR) $(DATABASEDIR) +ALLDIR= $(RESULTSDIR) $(APPLGRIDDIR) $(DATABASEDIR) $(PINEAPPLGRID) # root -ROOTINCS = $(shell root-config --cflags) -ROOTLIBS = $(shell root-config --libs) +ROOTINCS = $(shell root-config --cflags) +ROOTLIBS = $(shell root-config --libs) # APFEL -APFELINCS = $(shell apfel-config --cppflags) -APFELLIBS = $(shell apfel-config --ldflags) +APFELINCS = $(shell apfel-config --cppflags) +APFELLIBS = $(shell apfel-config --ldflags) #LHAPDF LHAPDFINCS = -I$(shell lhapdf-config --prefix)/include @@ -29,7 +30,7 @@ LHAPDFLIBS = -L$(LHAPDFDIR) -lLHAPDF # applgrid APPLINCS = -I$(shell applgrid-config --prefix)/include -APPLCLIBS = -L$(shell applgrid-config --prefix)/lib -lAPPLgrid +APPLCLIBS = -L$(shell applgrid-config --prefix)/lib -lAPPLgrid # libnnpdf NNPDFINCLUDE=$(shell pkg-config nnpdf --cflags) @@ -39,8 +40,12 @@ NNPDFLIBS=$(shell pkg-config nnpdf --libs) GSLINCLUDE=$(shell gsl-config --cflags) GSLLIBS=$(shell gsl-config --libs) -# additional libraries to be included -PRJLDFLAGS = $(LHAPDFLIBS) $(APPLCLIBS) $(ROOTLIBS) $(APFELLIBS) $(NNPDFLIBS) $(GSLLIBS) -lsqlite3 +# pineappl +PINEAPPLINCLUDE=$(shell pkg-config pineappl_capi --cflags) +PINEAPPLLIBS=$(shell pkg-config pineappl_capi --libs) + +# additional libraries to be included +PRJLDFLAGS = $(LHAPDFLIBS) $(APPLCLIBS) $(ROOTLIBS) $(APFELLIBS) $(NNPDFLIBS) $(GSLLIBS) $(PINEAPPLLIBS) -lsqlite3 # scheduling and optimization options (such as -DSSE -DSSE2 -DP4) -PRJCXXFLAGS = -Wall -O3 $(ALLDIR) $(LHAPDFINCS) $(APPLINCS) $(ROOTINCS) $(APFELINCS) $(NNPDFINCLUDE) $(GSLINCLUDE) +PRJCXXFLAGS = -Wall -O3 $(ALLDIR) $(LHAPDFINCS) $(APPLINCS) $(ROOTINCS) $(APFELINCS) $(NNPDFINCLUDE) $(GSLINCLUDE) $(PINEAPPLINCLUDE) diff --git a/db/apfelcomb.dat b/db/apfelcomb.dat index a4dc3ff..4a43cd4 100644 --- a/db/apfelcomb.dat +++ b/db/apfelcomb.dat @@ -1,5 +1,16 @@ PRAGMA foreign_keys=OFF; BEGIN TRANSACTION; + +CREATE TABLE IF NOT EXISTS "pineappl_subgrids" ( + `id` INTEGER, + `fktarget` TEXT, + `pineapplgrid` TEXT, + `mask` TEXT, + `operators` TEXT, + PRIMARY KEY(id) +); +INSERT INTO pineappl_subgrids VALUES(170,'ATLASZHIGHMASS49FB','ATLASZHIGHMASS49FB.pineappl.lz4','1 1 1 1 1 1 1 1 1 1 1 1 1',NULL); + CREATE TABLE IF NOT EXISTS "app_subgrids" ( `id` INTEGER, `fktarget` TEXT, @@ -590,6 +601,7 @@ INSERT INTO dyp_subgrids VALUES(20,'POSDYUD',NULL); INSERT INTO dyp_subgrids VALUES(21,'POSDYUDB',NULL); INSERT INTO dyp_subgrids VALUES(22,'POSDYUS',NULL); INSERT INTO dyp_subgrids VALUES(23,'POSDYUSB',NULL); + CREATE TABLE IF NOT EXISTS "grids" ( `id` INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT UNIQUE, `setname` TEXT NOT NULL, @@ -602,7 +614,7 @@ CREATE TABLE IF NOT EXISTS "grids" ( INSERT INTO grids VALUES(0,'ATLASWZRAP36PB','ATLASWZRAP36PB',' ATLAS 36PB W/Z rapidity',25,0,'APP'); INSERT INTO grids VALUES(1,'ATLASR04JETS2P76TEV','ATLASR04JETS2P76TEV','ATLAS 2.76 TeV R=0.4 inclusive jets',30,0,'APP'); INSERT INTO grids VALUES(2,'ATLASR04JETS36PB','ATLASR04JETS36PB','ATLAS 7 TeV R=0.4 inclusive jets',25,0,'APP'); -INSERT INTO grids VALUES(3,'ATLASZHIGHMASS49FB','ATLASZHIGHMASS49FB','ATLAS 49FB High mass Z',20,0,'APP'); +INSERT INTO grids VALUES(3,'ATLASZHIGHMASS49FB','ATLASZHIGHMASS49FB','ATLAS 49FB High mass Z',50,0,'APP'); INSERT INTO grids VALUES(4,'ATLASWPT31PB','ATLASWPT31PB_WP','ATLAS 31PB W+ pT',25,0,'APP'); INSERT INTO grids VALUES(5,'ATLASWPT31PB','ATLASWPT31PB_WM','ATLAS 31PB W- pT',35,0,'APP'); INSERT INTO grids VALUES(6,'ATLASWPT31PB','ATLASWPT31PB_WP_TOT','ATLAS 31PB W+ total',20,0,'APP'); @@ -894,6 +906,18 @@ INSERT INTO grids VALUES(294,'INTEGXV8','INTEGXV8','xV8 integrability',70,1,'DIS INSERT INTO grids VALUES(295,'INTEGXV','INTEGXV','xV integrability',70,1,'DIS'); INSERT INTO grids VALUES(296,'ATLAS_TTBARTOT_13TEV_FULLLUMI','ATLAS_TTBARTOT_13TEV_FULLLUMI','ATLAS 13 TeV inclusive ttbar cross-section',30,0,'APP'); +/* TODO: FIX THIS MESS */ +CREATE TABLE IF NOT EXISTS "fixmegrids" ( + `id` INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT UNIQUE, + `setname` TEXT NOT NULL, + `name` TEXT NOT NULL UNIQUE, + `description` TEXT NOT NULL, + `nx` INTEGER NOT NULL, + `positivity` INTEGER NOT NULL, + `source` TEXT NOT NULL +); +INSERT INTO fixmegrids VALUES(0,'ATLASZHIGHMASS49FB','ATLASZHIGHMASS49FB','ATLAS 49FB High mass Z',50,0,'PINE'); + DELETE FROM sqlite_sequence; INSERT INTO sqlite_sequence VALUES('dis_subgrids',71); INSERT INTO sqlite_sequence VALUES('dyp_subgrids',23); diff --git a/src/apfel_comb.cc b/src/apfel_comb.cc index f9e2473..6211c8f 100644 --- a/src/apfel_comb.cc +++ b/src/apfel_comb.cc @@ -35,7 +35,7 @@ int main(int argc, char* argv[]) { if (argc!=4) { - cout << "Usage: "< "< "< +#include +#include +#include + +// NNPDF +#include "NNPDF/fkgenerator.h" +#include "NNPDF/fastkernel.h" +#include "NNPDF/nnpdfdb.h" + +// PineAPPL includes +#include "pineappl_capi.h" + +#include "fk_utils.h" +#include "fk_qcd.h" +#include "fk_grids.h" + +using QCD::qcd_param; + +namespace PINE +{ + class SubGrid: public FKSubGrid + { + public: + void Compute(qcd_param const&, vector&) const; //!< Compute APPLgrid results mapped to Commondata + void Combine(QCD::qcd_param const&, NNPDF::FKGenerator*) const; //!< Perform the FK combination on a subgrid + private: + friend class ::FKTarget; + SubGrid(FKTarget const& parent, NNPDF::IndexDB const& db, int const& iDB); + + void Splash(ostream&) const; //!< Print metadata to stream + size_t GetNdat() const {return ndata;}; //!< Return number of datapoints in subgrid + double GetQ2max() const; //!< Return maximum scale used in a subgrid + double GetXmin() const; //!< Return minimum x-value used in this sub grid + double GetComputeXmin() const; //!< Return minimal x-value used in computation of this observable + + // *********************************************************** + + const string pineapplfile; //!< Path for the pineapplg file + const string readme; //!< APPLgrid README + + const vector maskmap; //!< Map of masked applgrid points to datapoints + const size_t ndata; //!< Number of selected datapoints + + const pineappl_grid* grid; //!< PineAPPL class + + static vector parse_maskmap(string const&); + }; + +} + diff --git a/src/apfelcomb/fk_qcd.h b/src/apfelcomb/fk_qcd.h index 62ceaa4..bb3efb4 100644 --- a/src/apfelcomb/fk_qcd.h +++ b/src/apfelcomb/fk_qcd.h @@ -109,6 +109,7 @@ namespace QCD // APFEL PDF access functions void evolpdf(const double& x, const double& Q, double* pdf); + double flvpdf(int32_t pdgid, double x, double q2); double alphas(const double& Q); double beta0(); diff --git a/src/apfelcomb/fk_utils.h b/src/apfelcomb/fk_utils.h index 707d4b5..a577649 100644 --- a/src/apfelcomb/fk_utils.h +++ b/src/apfelcomb/fk_utils.h @@ -22,6 +22,10 @@ typedef std::chrono::system_clock::duration time_span; #define RESULTS_PATH run/ #endif +#ifndef PINEAPPL_PATH +#define PINEAPPL_PATH ./ +#endif + #ifndef APPL_PATH #define APPL_PATH ./ #endif @@ -60,6 +64,7 @@ namespace Colour { void DisplayHR(); void Splash(); +std::string pineapplPath(); std::string applPath(); std::string dataPath(); std::string resultsPath(); diff --git a/src/core/Makefile b/src/core/Makefile index bd4198c..1792962 100644 --- a/src/core/Makefile +++ b/src/core/Makefile @@ -6,7 +6,7 @@ CXXFLAGS = $(PRJCXXFLAGS) -I .. LDFLAGS = $(PRJLDFLAGS) CORELIB = ../libac_core.a -COREOBJ = fk_appl.o fk_dis.o fk_qcd.o fk_utils.o fk_grids.o fk_ftdy.o +COREOBJ = fk_pine.o fk_appl.o fk_dis.o fk_qcd.o fk_utils.o fk_grids.o fk_ftdy.o all : $(CORELIB) $(CORELIB) : $(COREOBJ) diff --git a/src/core/fk_grids.cc b/src/core/fk_grids.cc index 7dbdcb2..90bbf59 100644 --- a/src/core/fk_grids.cc +++ b/src/core/fk_grids.cc @@ -10,6 +10,7 @@ #include "NNPDF/nnpdfdb.h" #include "NNPDF/commondata.h" +#include "apfelcomb/fk_pine.h" #include "apfelcomb/fk_appl.h" #include "apfelcomb/fk_dis.h" #include "apfelcomb/fk_ftdy.h" @@ -61,9 +62,10 @@ void FKTarget::Combine(QCD::qcd_param const& par, NNPDF::FKGenerator* FK) const FKTarget::source FKTarget::parse_source(string const& ss) { - if (ss.compare("APP") == 0) return FKTarget::APP; - if (ss.compare("DIS") == 0) return FKTarget::DIS; - if (ss.compare("DYP") == 0) return FKTarget::DYP; + if (ss.compare("PINE") == 0) return FKTarget::PINE; + if (ss.compare("APP") == 0) return FKTarget::APP; + if (ss.compare("DIS") == 0) return FKTarget::DIS; + if (ss.compare("DYP") == 0) return FKTarget::DYP; return FKTarget::NSR; } @@ -73,6 +75,9 @@ void FKTarget::ReadSubGrids(NNPDF::IndexDB const& db) for (int i : subgridIDs) switch(subgrid_source) { + case PINE: + components.insert(pair(i, new PINE::SubGrid(*this, db, i))); + break; case APP: components.insert(pair(i, new APP::SubGrid(*this, db, i))); break; diff --git a/src/core/fk_pine.cc b/src/core/fk_pine.cc new file mode 100644 index 0000000..4913fca --- /dev/null +++ b/src/core/fk_pine.cc @@ -0,0 +1,462 @@ +// Classes for A matrices (combined evolution-rotation) +// and Sigma matrices (combined applgrid - A matrices) +// n.p.hartland@ed.ac.uk - 03/12 + +#include "apfelcomb/fk_pine.h" +#include "apfelcomb/fk_qcd.h" +#include "apfelcomb/fk_utils.h" + +#include "NNPDF/common.h" +#include "NNPDF/nnpdfdb.h" +#include "NNPDF/commondata.h" + +#include +#include +#include +#include +#include +#include + +using namespace std; +using NNPDF::FKHeader; + +namespace PINE +{ + +// // ********************* Evolution factors ****************************** + + + class EvolutionFactors + { + public: + EvolutionFactors(const int nxin, const int nxout): + b1(14), + b2(14*14), + b3(14*14*nxin), + data(nxout*b3) + {} + + double* operator()(size_t ox, size_t ix, size_t fi ) {return &data.at(b3*ox+b2*ix+b1*fi);} + const double* operator()(size_t ox, size_t ix, size_t fi) const {return &data.at(b3*ox+b2*ix+b1*fi);} + private: + const size_t b1; + const size_t b2; + const size_t b3; + std::vector data; + }; + +// ********************************* Basis rotation helpers ************************************* + + // Rotates APFEL flavour basis into APPLgrid flavour basis (photon moves from 0 to 13) + double evolpdf_pineapplgrid(int32_t pdgid, double x, double q2, void*) + { + if (x < APFEL::xGrid(0)) // A nice trick of APPLgrid is to request PDF x-values smaller than are actually used + return 0; + return QCD::flvpdf(pdgid, x, q2); + } + + double alphas(double q2, void*) + { + return QCD::alphas(sqrt(q2)); + } + + // convert pineappl pdg to applgrid + int32_t pdgid_to_apfel_array_id(int32_t pdgid) + { + switch (pdgid) + { + case 22: + return 13; + + case 21: + return 6; + + case -6: + case -5: + case -4: + case -3: + case -2: + case -1: + case 0: + case 6: + case 5: + case 4: + case 3: + case 2: + case 1: + return pdgid + 6; + + default: + throw std::runtime_error("pdgid `" + std::to_string(pdgid) + "` not available."); + } + } + + // ********************************* Kinematics Helpers ************************************* +/* + // Returns the minimum and maximum x-grid points for a specified subgrid slice. + // igrid is the requested subgrid, nsubproc the number of subprocesses held within igrid. + // tau specified the bin in scale to be investigated, and alpha specifies the bin in x1. + // return.first and return.second return the minumum and maxiumum bins in x2 respectively. + std::pair get_slice_limits(appl::igrid const* igrid, int const& nsubproc, int const& tau, int const& alpha) + { + std::pair limits(igrid->Ny2(), 0); + appl::igrid* igrid_nc = const_cast(igrid); + for (int tsp=0; tspweightgrid(tsp))[tau] != NULL) + for (int ix2=0; ix2Ny2(); ix2++) + if ( (*(const SparseMatrix3d*) igrid_nc->weightgrid(tsp))(tau,alpha,ix2) != 0 ) + { + limits.first = std::min(ix2, limits.first); + limits.second = std::max(ix2, limits.second); + } + + return limits; + } + + // Returns the minimum and maximum used values of x1 across all grid slices + std::pair get_igrid_limits_x1(appl::igrid const* igrid, int const& nsubproc, int const& tau) + { + std::pair limits(igrid->Ny1(), 0); + for (int alpha = 0; alpha < igrid->Ny1(); alpha++) + { + const std::pair sl = get_slice_limits(igrid, nsubproc, tau, alpha); + if (sl.first <= sl.second) // This alpha is not trimmed + { + limits.first = std::min(alpha, limits.first); + limits.second = std::max(alpha, limits.second); + } + } + return limits; + } + + // Returns the minimum and maximum used values of x2 across all grid slices + std::pair get_igrid_limits_x2(appl::igrid const* igrid, int const& nsubproc, int const& tau) + { + std::pair limits(igrid->Ny2(), 0); + for (int alpha = 0; alpha < igrid->Ny1(); alpha++) + { + const std::pair sl = get_slice_limits(igrid, nsubproc, tau, alpha); + limits.first = std::min(sl.first, limits.first); + limits.second = std::max(sl.second, limits.second); + } + return limits; + } + + + // ********************************* Combination Helpers ************************************* + + // Translates 'loop' order to appl::grid index + // This is specifically in order to translate aMC@NLO-like four-part grids + // into LO and NLO components. + // aMC@NLO-like convolution uses Born = grid 3 + // NLO = grid 0 + int get_grid_idx( const appl::grid *g, int const& pto ) + { + if (g->calculation() == appl::grid::AMCATNLO) + return (pto==0) ? 3:0; + return pto; + } + + // Returns the APPLgrid PDF object associated with the ith subgrid of g + appl::appl_pdf* get_appl_pdf( const appl::grid *g, int const& i ) + { + // Split PDF string + const std::string pdfnames = g->getGenpdf(); + std::vector pdfvec; + std::stringstream ss(pdfnames); std::string s; + while (getline(ss, s, ':')) pdfvec.push_back(s); + + const size_t isubproc = pdfvec.size() == 1 ? 0:i; + return appl::appl_pdf::getpdf( pdfvec[isubproc] ); + } + + // Computes the normalisation required for translating APPLgrid weights to FK ones + // e.g including factors of alpha_S, bin width etc. + // g is the APPLgrid being combined with evolution factors. + // pto specified the perturbative order being combined, as the value of alpha_S in the current bin, + // and x1/x2 specify the numerical values of the PDF x-values for the first and second PDF respectively. + double compute_wgt_norm( const appl::grid *g, int const& d, double const& pto, double const& as, double const& x1, double const& x2) + { + // PDF x and bin width normalisation + double norm = 1.0/(x1*x2*g->deltaobs(d)); + + // Normalisation by number of runs + appl::grid* g_nc = const_cast(g); + if ( !g->getNormalised() && g_nc->run() ) + norm*=1.0/double(g_nc->run()); + + // Factor of alpha_S + const double LO = g->leadingOrder(); + if (g->calculation() == appl::grid::AMCATNLO) + { norm*=pow(as*(4.0*M_PI), LO+pto ); } + else + { norm*=pow(as/(2.0*M_PI), LO+pto ); } + + return norm; + } + +// // Progress update **************************************************** + + int countElements(vector const& maskmap, int const& min_pto, int const& max_pto, const appl::grid* g) + { + // Counter + int nElm = 0; + for (auto bin : maskmap ) + for (int pto=min_pto; pto<=max_pto; pto++) + { + const int gidx = get_grid_idx(g, pto); // APPLgrid grid index + const appl::igrid *igrid = g->weightgrid(gidx, bin); + const size_t nsubproc = g->subProcesses(gidx); // Number of subprocesses + for (int t=0; tNtau(); t++) // Loop over scale bins + for (int a=0; aNy1(); a++ ) // Loop over x1 bins + { + const std::pair sl = get_slice_limits(igrid, nsubproc, t, a); + nElm += std::max(0, sl.second - sl.first + 1); + } + } + return nElm; + } +*/ + + // *********************************** APPLGrid convolutions ************************************** + + SubGrid::SubGrid(FKTarget const& parent, NNPDF::IndexDB const& db, int const& iDB): + FKSubGrid(parent, iDB, NNPDF::dbquery(db, iDB, "operators")), + pineapplfile(NNPDF::dbquery(db,iDB,"pineapplgrid")), + readme(), + maskmap(parse_maskmap(NNPDF::dbquery(db,iDB,"mask"))), + ndata(maskmap.size()), + grid(pineappl_grid_read((pineapplPath() + parent.GetSetName() + "/" + pineapplfile).c_str())) + { + const size_t bin_count = pineappl_grid_bin_count(grid); + if (ndata > bin_count) + { + cerr <<"Error: number of datapoints in settings: "< appl grid Nobs: "<& results) const + { + if (par.evol_pto == 2) + std::cout << "WARNING: APPLgrid does not currently support NNLO convolutions, fixing convolution to NLO" < xsec(pineappl_grid_bin_count(grid)); + pineappl_grid_convolute(grid, evolpdf_pineapplgrid, evolpdf_pineapplgrid, alphas, + nullptr, nullptr, nullptr, par.xiR, par.xiF, xsec.data()); + + // apply database normalization rule + for (double& obs : xsec) + obs *= nrmdat; + + // Relate back to results + for (size_t i=0; i afl = QCD::active_flavours(par); + + vector orders(4 * pineappl_grid_order_count(grid)); + pineappl_grid_order_params(grid, orders.data()); + + auto *grid_lumi = pineappl_grid_lumi(grid); + const size_t lumi_count = pineappl_lumi_count(grid_lumi); + + vector q2values(pineappl_subgrid_q2_grid_count(grid)); + pineappl_subgrid_q2_grid(grid, q2values.data()); + + vector xvalues(pineappl_subgrid_x_grid_count(grid)); + pineappl_subgrid_x_grid(grid, xvalues.data()); + const size_t nxpi = xvalues.size(); + + vector bin_sizes(pineappl_grid_bin_count(grid)); + pineappl_grid_bin_sizes(grid, bin_sizes.data()); + + vector weight_matrix(nxpi * nxpi); + + vector pdgids; + vector factors; + + int completedElements = 0; + int nXelements = 0; + for (size_t d=0; d 0 || orders.at(4 * pto + 3) > 0) + continue; + + for (size_t lumi = 0; lumi < lumi_count; lumi++) + { + std::array slice_indices; + pineappl_subgrid_filled_q2_slices(grid, pto, bin, lumi, slice_indices.data()); + nXelements += (slice_indices.at(1) - slice_indices.at(0)) * nxpi * nxpi; + } + } + } + + const time_point t1 = std::chrono::system_clock::now(); + + for (size_t d=0; d 0 || orders.at(4 * pto + 3) > 0) + continue; + + for (size_t lumi = 0; lumi < lumi_count; lumi++) + { + // prepare luminosity combinations + const size_t combinations = pineappl_lumi_combinations(grid_lumi, lumi); + pdgids.resize(2 * combinations); + factors.resize(combinations); + pineappl_lumi_entry(grid_lumi, lumi, pdgids.data(), factors.data()); + + // Fetch grid pointer and loop over Q + std::array slice_indices; + pineappl_subgrid_filled_q2_slices(grid, pto, bin, lumi, slice_indices.data()); + for (size_t t = slice_indices.at(0); t < slice_indices.at(1); t++) + { + // Scales and strong coupling + const double Q2 = q2values.at(t); + const double QF = sqrt(Q2)*par.xiF; + const double QR = sqrt(Q2)*par.xiR; + const double as = QCD::alphas(QR); + + // Renormalisation and factorisation scale variation terms + //const bool vary_ren = pto == 0 && par.evol_pto == 1 && par.xiR != 1.0; + //const bool vary_fac = pto == 0 && par.evol_pto == 1 && par.xiF != 1.0; + //const double renscale = (as/(2.0*M_PI))*2.0*M_PI*QCD::beta0()*g->leadingOrder()*log(par.xiR*par.xiR); + //const double facscale = -(as/(2.0*M_PI))*log(par.xiF*par.xiF); + + // define evolution factor arrays + const int nxin = fk->GetNx(); + EvolutionFactors fA(nxin, nxpi); // PDF 1 and 2 Evolution factors + //EvolutionFactors fdA1(nxin, igrid->Ny1()); // PDF 1 Splitting function x Evolution factors + //EvolutionFactors fdA2(nxin, igrid->Ny2()); // PDF 2 Splitting function x Evolution factors + + // Compute nonzero evolution factors + for (int ix = 0; ix < nxin; ix++) + for (size_t fl : afl) + { + for (size_t ox = 0; ox < nxpi; ox++) // Loop over applgrid x1 + { + QCD::EvolutionOperator(false, true, ix, xvalues.at(ox), fl, QF, fA(ox, ix, fl)); + //if (vary_fac) QCD::DerivativeOperator(false,applgrid_nfl==14,ix,igrid->fx(igrid->gety1(ox)),fl,QF,fdA1(ox,ix,fl)); + } + } + + // take the grid from pineappl + pineappl_subgrid_q2_slice(grid, pto, bin, lumi, t, weight_matrix.data()); + + for (size_t a = 0; a < nxpi; a++) + for (size_t b = 0; b < nxpi; b++) // Loop over applgrid x2 + { + // fetch weight values + const double W = weight_matrix.at(a * nxpi + b); + + // Calculate normalisation factors + const double norm = pow(as, orders.at(4 * pto + 0))/bin_sizes.at(bin)*nrmdat; + + for (int i = 0; i < nxin; i++) // Loop over input pdf x1 + for (int j = 0; j < nxin; j++) // Loop over input pdf x2 + for (size_t k : afl) // loop over flavour 1 + for (size_t l : afl) // loop over flavour 2 + { + double H = 0.0; + for (std::size_t c = 0; c != combinations; ++c) + { + auto const ai = pdgid_to_apfel_array_id(pdgids.at(2 * c + 0)); + auto const bi = pdgid_to_apfel_array_id(pdgids.at(2 * c + 1)); + H += factors.at(c) * (fA(a, i, k)[ai]) * (fA(b, j, l)[bi]); + } + + // if (vary_fac) + // { + // genpdf->evaluate(fdA1(a,i,k),fA2(b,j,l),H1); + // genpdf->evaluate(fA1(a,i,k),fdA2(b,j,l),H2); + // } + + double fill = H; // Basic fill + //if (vary_ren) fill += renscale*H[ip]; // Ren. scale variation + //if (vary_fac) fill += facscale*(H1[ip] + H2[ip]); // Fac. scale variation + if ( fill != 0.0 ) + for (int const& td : datamap[d]) + fk->Fill(td, i, j, k, l, norm*fill*W); + } + + // Update progress + completedElements++; + StatusUpdate(t1, (double)completedElements/(double)nXelements, cout); + } + } + } + } // /pto + } // /data + + pineappl_lumi_delete(grid_lumi); + + return; + } + + // *************************************** Metadata methods ******************************************** + + + void SubGrid::Splash(ostream& o) const + { + FKSubGrid::Splash(o); + o << "- PineAPPL: " << pineapplfile << endl; + } + + // Get the maximum scale of an applgrid + double SubGrid::GetQ2max() const + { + vector q2values(pineappl_subgrid_q2_grid_count(grid)); + pineappl_subgrid_q2_grid(grid, q2values.data()); + return std::max(q2values.front(), q2values.back()); + } + + // Return the minimum x used in the subgrid + double SubGrid::GetXmin() const + { + // TODO: check this implementation. + return GetComputeXmin(); + }; + + double SubGrid::GetComputeXmin() const + { + vector xvalues(pineappl_subgrid_x_grid_count(grid)); + pineappl_subgrid_x_grid(grid, xvalues.data()); + return std::min(xvalues.front(), xvalues.back()); + }; + + vector SubGrid::parse_maskmap(string const& mask) + { + const vector masksplit = ssplit(mask); + vector _maskmap; + for (size_t i=0; i