From b04ebd0d6467a408d675be0ca2a677bfc90b39d3 Mon Sep 17 00:00:00 2001
From: Stefano Carrazza <stefano.carrazza@cern.ch>
Date: Wed, 26 Aug 2020 10:58:02 +0200
Subject: [PATCH 01/14] adding pineappl related structure

---
 Makefile.inc             |  20 +-
 src/apfel_comb.cc        |   2 +-
 src/apfelcomb/fk_grids.h |   2 +-
 src/apfelcomb/fk_pine.h  |  74 ++++++
 src/core/Makefile        |   2 +-
 src/core/fk_grids.cc     |  11 +-
 src/core/fk_pine.cc      | 490 +++++++++++++++++++++++++++++++++++++++
 7 files changed, 587 insertions(+), 14 deletions(-)
 create mode 100644 src/apfelcomb/fk_pine.h
 create mode 100644 src/core/fk_pine.cc

diff --git a/Makefile.inc b/Makefile.inc
index e74a3fe..1937867 100644
--- a/Makefile.inc
+++ b/Makefile.inc
@@ -15,12 +15,12 @@ DATABASEDIR=  -D  DB_PATH="./db/"
 ALLDIR= $(RESULTSDIR) $(APPLGRIDDIR) $(DATABASEDIR)
 
 # 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 +29,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 +39,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/src/apfel_comb.cc b/src/apfel_comb.cc
index f9e2473..0b5bf11 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: "<<argv[0]<<" <source=app/dis/dyp> <database id> <theory id>"<<endl;
+    cout << "Usage: "<<argv[0]<<" <source=pine/app/dis/dyp> <database id> <theory id>"<<endl;
     exit(1);
   }
 
diff --git a/src/apfelcomb/fk_grids.h b/src/apfelcomb/fk_grids.h
index 128ea55..30e3ff3 100644
--- a/src/apfelcomb/fk_grids.h
+++ b/src/apfelcomb/fk_grids.h
@@ -27,7 +27,7 @@ class FKTarget
 public:
 	FKTarget(NNPDF::IndexDB const& db, int const& targetID, int const& global_nx);
 	~FKTarget(){};
-	enum source { APP, DIS, DYP, NSR };
+	enum source { PINE, APP, DIS, DYP, NSR };
 
 	virtual void Splash( ostream& ) const; 				//!< Write table information to stream
 	void ReadSubGrids(NNPDF::IndexDB const& db);		//!< Read information on subgrids from database
diff --git a/src/apfelcomb/fk_pine.h b/src/apfelcomb/fk_pine.h
new file mode 100644
index 0000000..e86fe3e
--- /dev/null
+++ b/src/apfelcomb/fk_pine.h
@@ -0,0 +1,74 @@
+#pragma once
+
+#include <stdlib.h>
+#include <fstream>
+#include <vector>
+#include <sys/stat.h>
+
+// NNPDF
+#include "NNPDF/fkgenerator.h"
+#include "NNPDF/fastkernel.h"
+#include "NNPDF/nnpdfdb.h"
+
+// APPLGrid main includes
+#include "appl_grid/appl_grid.h"
+#include "appl_grid/appl_igrid.h"
+#include "appl_grid/appl_pdf.h"
+#include "appl_grid/fastnlo.h"
+
+#include "fk_utils.h"
+#include "fk_qcd.h"
+#include "fk_grids.h"
+
+using QCD::qcd_param;
+
+namespace PINE
+{
+
+	// // Wrapper for appl::grids (handles fastNLO based grids)
+	class grid
+	{
+	public:
+		grid(std::string const& filename, int const& fnlobin);
+		~grid();
+
+		fastnlo* 	fg;
+		appl::grid* g;
+	};
+
+
+	class SubGrid: public FKSubGrid
+	{
+	public:
+	void Compute(qcd_param const&, vector<double>&) 			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 		applfile;	//!< Path for the applgrid file
+		const string 		readme;		//!< APPLgrid README
+
+		const vector<int> 	maskmap;	//!< Map of masked applgrid points to datapoints
+		const size_t 		ndata;     	//!< Number of selected datapoints
+		const size_t 		ptmin;     	//!< Minimum perturbative order to contribute
+
+		const int    		fnlobin;   	//!< Which bin in the fastNLO grid? (yes, this needs to be signed)
+		const bool 			pdfwgt;     //!< Using a PDF weight? Typically used for fastNLO grids
+		const bool 			ppbar;		//!< Does the grid need a pp -> ppbar transform?
+
+		const grid 			applgrid;	//!< APPLgrid class
+
+		static vector<int>  parse_maskmap(string const&);
+	};
+
+}
+
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<int, FKSubGrid*>(i, new PINE::SubGrid(*this, db, i)));
+				break;
 			case APP:
 				components.insert(pair<int, FKSubGrid*>(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..aa5cb13
--- /dev/null
+++ b/src/core/fk_pine.cc
@@ -0,0 +1,490 @@
+// 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 <chrono>
+#include <ctime>
+#include <limits>
+
+using namespace std;
+using NNPDF::FKHeader;
+
+namespace PINE
+{
+  #ifdef __APPL_PHOTON
+  const int applgrid_nfl = 14;
+  #else
+  const int applgrid_nfl = 13;
+  #endif
+
+
+// // ********************* Evolution factors ******************************
+
+  class EvolutionFactors
+  {
+  public:
+    EvolutionFactors(const int nxin, const int nxout):
+    b1(applgrid_nfl),
+    b2(applgrid_nfl*14),
+    b3(applgrid_nfl*14*nxin),
+    data(new double[nxout*b3])
+    {
+      for (int i=0; i<nxout*b3; i++)
+        data[i] = 0;
+    };
+    ~EvolutionFactors() {delete[] data;};
+
+    double* operator()(int const& ox, int const& ix, int const& fi ) {return data+b3*ox+b2*ix+b1*fi;};
+    const double* operator()(int const& ox, int const& ix, int const& fi ) const {return data+b3*ox+b2*ix+b1*fi;};
+  private:
+    const int b1;
+    const int b2;
+    const int b3;
+    double* data;
+  };
+
+// ********************************* Basis rotation helpers *************************************
+
+  // Rotates APFEL flavour basis into APPLgrid flavour basis (photon moves from 0 to 13)
+  void evolpdf_applgrid(const double& x, const double& Q, double* pdf)
+  {
+    for (int i=0; i<applgrid_nfl; i++) pdf[i]=0;
+    if (x<APFEL::xGrid(0))  // A nice trick of APPLgrid is to request PDF x-values smaller than are actually used
+      return;
+
+    double *APFEL_basis = new double[14]();
+    QCD::evolpdf(x, Q, APFEL_basis);
+
+    for (int i=0; i<13; i++)
+      pdf[i]=APFEL_basis[i+1];
+    if (applgrid_nfl == 14)
+      pdf[13] = APFEL_basis[0];
+
+    delete[] APFEL_basis;
+  }
+
+  // antiproton version of evolpdf_applgrid, swaps q<->qbar
+  void evolpdf_applgrid_pbar(const double& x, const double& Q, double* pdf)
+  {
+    evolpdf_applgrid(x,Q,pdf);
+    for (int i=1; i<7; i++)
+      std::swap(pdf[6+i], pdf[6-i]);
+  }
+
+  // ********************************* 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<int,int> get_slice_limits(appl::igrid const* igrid, int const& nsubproc, int const& tau, int const& alpha)
+  {
+    std::pair<int,int> limits(igrid->Ny2(), 0);
+    appl::igrid* igrid_nc = const_cast<appl::igrid*>(igrid);
+    for (int tsp=0; tsp<nsubproc; tsp++)
+      if ((*(const SparseMatrix3d*) igrid_nc->weightgrid(tsp))[tau] != NULL)
+        for (int ix2=0; ix2<igrid->Ny2(); 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<int,int> get_igrid_limits_x1(appl::igrid const* igrid, int const& nsubproc, int const& tau)
+  {
+    std::pair<int,int> limits(igrid->Ny1(), 0);
+    for (int alpha = 0; alpha < igrid->Ny1(); alpha++)
+    {
+      const std::pair<int,int> 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<int,int> get_igrid_limits_x2(appl::igrid const* igrid, int const& nsubproc, int const& tau)
+  {
+    std::pair<int,int> limits(igrid->Ny2(), 0);
+    for (int alpha = 0; alpha < igrid->Ny1(); alpha++)
+    {
+      const std::pair<int,int> 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<std::string> 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<appl::grid*>(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<int> 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; t<igrid->Ntau(); t++)     // Loop over scale bins
+          for (int a=0; a<igrid->Ny1(); a++  )  // Loop over x1 bins
+          {
+            const std::pair<int,int> sl = get_slice_limits(igrid, nsubproc, t, a);
+            nElm += std::max(0, sl.second - sl.first + 1);
+          }
+      }
+    return nElm;
+  }
+
+ // *********************************** Grid object **************************************
+
+  grid::grid(std::string const& filename, int const& fnlobin):
+    fg(fnlobin >= 0 ? new fastnlo(filename):0),
+    g(fg ? fg->grids()[fnlobin]:new appl::grid(filename))
+    {}
+
+  grid::~grid() {
+    fg ? delete fg : delete g;
+  }
+
+
+ // *********************************** APPLGrid convolutions **************************************
+
+  SubGrid::SubGrid(FKTarget const& parent, NNPDF::IndexDB const& db, int const& iDB):
+    FKSubGrid(parent, iDB, NNPDF::dbquery<string>(db, iDB, "operators")),
+    applfile(NNPDF::dbquery<string>(db,iDB,"applgrid")),
+    readme(),
+    maskmap(parse_maskmap(NNPDF::dbquery<string>(db,iDB,"mask"))),
+    ndata(	maskmap.size() ),
+    ptmin(	NNPDF::dbquery<size_t>(db,iDB,"ptmin")),
+    fnlobin(NNPDF::dbquery<int>(db,iDB,"fnlobin")),
+    pdfwgt(	NNPDF::dbquery<bool>(db,iDB,"pdfwgt")),
+    ppbar(	NNPDF::dbquery<bool>(db,iDB,"ppbar")),
+    applgrid(applPath() + parent.GetSetName() + "/" + applfile , fnlobin)
+  {
+    if (ndata>static_cast<size_t>(applgrid.g->Nobs()))
+    {
+      cerr <<"Error: number of datapoints in settings: "<<ndata<<" > appl grid Nobs: "<<applgrid.g->Nobs()<<endl;
+      cerr <<"Please check the provided mask" <<endl;
+      exit(-1);
+    }
+  }
+
+  void SubGrid::Compute(qcd_param const& par, vector<double>& results) const
+  {
+    if (ppbar == true && par.xiF != 1)
+      std::cout << "WARNING: ppbar ROTATION NOT TESTED - APPLgrid does not support fac. scale variation with ppbar so I cannot cross-check" <<std::endl;
+    if (par.evol_pto == 2)
+      std::cout << "WARNING: APPLgrid does not currently support NNLO convolutions, fixing convolution to NLO" <<std::endl;
+
+    // Compute with applgrid interface
+    const int pto = (ptmin == 1) ? -1: min(par.evol_pto,(size_t)1);
+    vector<double> xsec;
+    if ( ppbar == true && par.xiF == 1)
+      xsec = applgrid.g->vconvolute( evolpdf_applgrid, evolpdf_applgrid_pbar, QCD::alphas, pto, par.xiR, par.xiF );
+    else
+      xsec = applgrid.g->vconvolute( evolpdf_applgrid, QCD::alphas, pto, par.xiR, par.xiF);
+    for (double& obs : xsec)
+      obs *= nrmdat;
+
+    // Relate back to results
+    for (size_t i=0; i<maskmap.size(); i++)
+      for (int const& j : datamap[i])
+        results[j] += xsec[maskmap[i]];
+  }
+
+
+//  *************************************** FK Table computation ********************************************
+
+  void SubGrid::Combine(QCD::qcd_param const& par, NNPDF::FKGenerator* fk) const
+  {
+    if (par.evol_pto == 2)
+      std::cout << "WARNING: APPLgrid does not currently support NNLO convolutions, fixing convolution to NLO" <<std::endl;
+
+    if (applgrid_nfl == 14)
+      std::cout << "WARNING: Combining with photon channel in APPLgrid" <<std::endl;
+
+    // APPLgrid pointer
+    const appl::grid* g = applgrid.g;
+    const int ptmax = min(par.evol_pto,(size_t)1); // Maximum perturbative order limited to NLO
+
+    // Progress monitoring
+    int completedElements = 0;
+    const int nXelements = countElements(maskmap, ptmin, ptmax, g);
+    const time_point t1 = std::chrono::system_clock::now();
+    const vector<size_t> afl = QCD::active_flavours(par);
+
+    for (size_t d=0; d<maskmap.size(); d++)
+    {
+      // Fetch associated applgrid info
+      const size_t bin = maskmap[d];
+      for (size_t pto=0; pto<=(ptmax-ptmin); pto++) // Loop over perturbative order
+      {
+        // Determine grid index, allocate subprocess arrays
+        const int gidx = get_grid_idx(g, pto+ptmin);
+        appl::appl_pdf *genpdf = get_appl_pdf( g, gidx );
+        const size_t nsubproc = g->subProcesses(gidx);
+        double *W = new double[nsubproc];   // Weight array
+        double *H = new double[nsubproc];   // Evolution factor array
+        double *H1 = new double[nsubproc];  // Split evolution factor array 1
+        double *H2 = new double[nsubproc];  // Split evolution factor array 2
+
+        // Fetch grid pointer and loop over Q
+        appl::igrid const *igrid = g->weightgrid(gidx, bin);
+        for (int t=0; t<igrid->Ntau(); t++)
+        {
+          // Scales and strong coupling
+          const double Q2  = igrid->fQ2( igrid->gettau(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 fA1(nxin, igrid->Ny1());   // PDF 1 Evolution factors
+          EvolutionFactors fA2(nxin, igrid->Ny2());   // PDF 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
+          const std::pair<int,int> l1 = get_igrid_limits_x1(igrid, nsubproc, t);
+          const std::pair<int,int> l2 = get_igrid_limits_x2(igrid, nsubproc, t);
+
+          for (int ix = 0; ix < nxin; ix++)
+          for (size_t fl : afl)
+          {
+            for (int ox=l1.first; ox<=l1.second; ox++) // Loop over applgrid x1
+            {
+              QCD::EvolutionOperator(false,applgrid_nfl==14,ix,igrid->fx(igrid->gety1(ox)),fl,QF,fA1(ox,ix,fl));
+              if (vary_fac) QCD::DerivativeOperator(false,applgrid_nfl==14,ix,igrid->fx(igrid->gety1(ox)),fl,QF,fdA1(ox,ix,fl));
+            }
+            for (int ox=l2.first; ox<=l2.second; ox++) // Loop over applgrid x2
+            {
+              QCD::EvolutionOperator(ppbar,applgrid_nfl==14,ix,igrid->fx(igrid->gety2(ox)),fl,QF,fA2(ox,ix,fl));
+              if (vary_fac) QCD::DerivativeOperator(ppbar,applgrid_nfl==14,ix,igrid->fx(igrid->gety2(ox)),fl,QF,fdA2(ox,ix,fl));
+            }
+          }
+
+          for (int a=l1.first; a<=l1.second; a++ )
+          {
+            const double x1 = igrid->fx(igrid->gety1(a));
+            const std::pair<int,int> limits = get_slice_limits(igrid, nsubproc, t, a);
+            for (int b=limits.first; b<=limits.second; b++) // Loop over applgrid x2
+            {
+              // fetch weight values
+              for (size_t ip=0; ip<nsubproc; ip++)
+                W[ip] = (*(const SparseMatrix3d*) const_cast<appl::igrid*>(igrid)->weightgrid(ip))(t,a,b);
+
+              // Calculate normalisation factors
+              const double x2 = igrid->fx(igrid->gety2(b));
+              const double pdfnrm = pdfwgt ? igrid->weightfun(x1)*igrid->weightfun(x2) : 1.0;
+              const double norm = pdfnrm*compute_wgt_norm(g, bin, pto+ptmin, as, x1, x2)*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
+                    {
+                      // Rotate to subprocess basis
+                      genpdf->evaluate(fA1(a,i,k),fA2(b,j,l),H);
+
+                      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);
+                      }
+
+                      for (size_t ip=0; ip<nsubproc; ip++)
+                        if (W[ip] != 0)
+                        {
+                          double fill = H[ip];                              // 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[ip]);
+                        }
+                    }
+
+              // Update progress
+              completedElements++;
+              StatusUpdate(t1, (double)completedElements/(double)nXelements, cout);
+            }
+          }
+        }
+
+        // Free subprocess arrays
+        delete[] W;
+        delete[] H;
+        delete[] H1;
+        delete[] H2;
+
+      } // /pto
+    } // /data
+    return;
+  }
+
+  //  *************************************** Metadata methods ********************************************
+
+
+  void SubGrid::Splash(ostream& o) const
+  {
+    FKSubGrid::Splash(o);
+    o << "- APPLgrid: " << applfile << endl
+      << "- PTMin: "    << ptmin    << endl
+      << "- fnlobin: "  << fnlobin  << endl
+      << "- PDFWeight: "<< pdfwgt   << endl
+      << "- ppbar: "    << ppbar    << endl
+      << endl;
+  }
+
+    // Get the maximum scale of an applgrid
+  double SubGrid::GetQ2max() const
+  {
+    // Find maximum required scale
+    double Q2max = 0;
+    const double nloops = applgrid.g->calculation() == appl::grid::AMCATNLO ? 4 : 2;
+    for(int i=0; i<nloops; i++)  // pto
+      for (int j=0; j<applgrid.g->Nobs(); j++) // bin
+      {
+        appl::igrid const *igrid = applgrid.g->weightgrid(i, j);
+        Q2max = max(Q2max, igrid->getQ2max());
+      }
+    return Q2max;
+  }
+
+
+  // Return the minimum x used in the subgrid
+  double SubGrid::GetXmin() const
+  {
+    const double nloops = applgrid.g->calculation() == appl::grid::AMCATNLO ? 4 : 2;
+    double xmin = 1.0;
+    for(int i=0; i<nloops; i++)
+      for (int j=0; j<applgrid.g->Nobs(); j++)
+      {
+        appl::igrid const *igrid = applgrid.g->weightgrid(i, j);
+        const size_t nsubproc = applgrid.g->subProcesses(i);
+
+          for (int t=0; t<igrid->Ntau(); t++) // Loop over Q^2 integral
+          {
+            const std::pair<int,int> l1 = get_igrid_limits_x1(igrid, nsubproc, t);
+            if (l1.first <= l1.second)
+            {
+              xmin = min(xmin, igrid->fx(igrid->gety1(l1.first)));
+              xmin = min(xmin, igrid->fx(igrid->gety1(l1.second)));
+            }
+            const std::pair<int,int> l2 = get_igrid_limits_x2(igrid, nsubproc, t);
+            if (l2.first <= l2.second)
+            {
+              xmin = min(xmin, igrid->fx(igrid->gety2(l2.first)));
+              xmin = min(xmin, igrid->fx(igrid->gety2(l2.second)));
+            }
+          }
+        }
+    return xmin;
+  };
+
+  double SubGrid::GetComputeXmin() const
+  {
+    const double nloops = applgrid.g->calculation() == appl::grid::AMCATNLO ? 4 : 2;
+    double xmin = 1.0;
+    for(int i=0; i<nloops; i++)
+      for (int j=0; j<applgrid.g->Nobs(); j++)
+      {
+        appl::igrid const *igrid = applgrid.g->weightgrid(i, j);
+        xmin = min(xmin, igrid->fx(igrid->gety1(0)));
+        xmin = min(xmin, igrid->fx(igrid->gety1(igrid->Ny1()-1)));
+        xmin = min(xmin, igrid->fx(igrid->gety2(0)));
+        xmin = min(xmin, igrid->fx(igrid->gety2(igrid->Ny2()-1)));
+      }
+
+    return xmin;
+  };
+
+
+  vector<int> SubGrid::parse_maskmap(string const& mask)
+  {
+    const vector<string> masksplit = ssplit(mask);
+    vector<int> _maskmap;
+    for (size_t i=0; i<masksplit.size(); i++)
+      if ((bool) atoi(masksplit[i].c_str()))
+        _maskmap.push_back(i);
+    return _maskmap;
+  }
+
+}

From 855e20cd4179fd0334995cd9bf85b326de56f0b5 Mon Sep 17 00:00:00 2001
From: Stefano Carrazza <stefano.carrazza@cern.ch>
Date: Wed, 26 Aug 2020 11:53:02 +0200
Subject: [PATCH 02/14] adding more interfacing

---
 Makefile.inc             |  3 ++-
 db/apfelcomb.dat         | 24 ++++++++++++++++++++
 src/apfel_comb.cc        |  5 +++--
 src/apfelcomb/fk_pine.h  | 31 +++++---------------------
 src/apfelcomb/fk_utils.h |  5 +++++
 src/core/fk_pine.cc      | 47 ++++++++++++++++------------------------
 src/core/fk_utils.cc     |  6 +++++
 7 files changed, 64 insertions(+), 57 deletions(-)

diff --git a/Makefile.inc b/Makefile.inc
index 1937867..e5798f5 100644
--- a/Makefile.inc
+++ b/Makefile.inc
@@ -10,9 +10,10 @@ 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)
diff --git a/db/apfelcomb.dat b/db/apfelcomb.dat
index a4dc3ff..c73dafd 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,
@@ -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',20,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 0b5bf11..345efad 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: "<<argv[0]<<" <source=pine/app/dis/dyp> <database id> <theory id>"<<endl;
+    cout << "Usage: "<<argv[0]<<" <source=pineappl/app/dis/dyp> <database id> <theory id>"<<endl;
     exit(1);
   }
 
@@ -49,7 +49,8 @@ int main(int argc, char* argv[]) {
   Splash(); QCD::qcd_param par; QCD::parse_input(iTh, par);
 
   NNPDF::SetVerbosity(0);
-  NNPDF::IndexDB grid_db(databasePath()+"apfelcomb.db", "grids");
+  std::string gridname = (source == "pineappl") ? "fixmegrids" : "grids";
+  NNPDF::IndexDB grid_db(databasePath()+"apfelcomb.db", gridname);
   NNPDF::IndexDB subgrid_db(databasePath()+"apfelcomb.db", source+"_subgrids");
 
   // Read grid information
diff --git a/src/apfelcomb/fk_pine.h b/src/apfelcomb/fk_pine.h
index e86fe3e..920382b 100644
--- a/src/apfelcomb/fk_pine.h
+++ b/src/apfelcomb/fk_pine.h
@@ -10,11 +10,8 @@
 #include "NNPDF/fastkernel.h"
 #include "NNPDF/nnpdfdb.h"
 
-// APPLGrid main includes
-#include "appl_grid/appl_grid.h"
-#include "appl_grid/appl_igrid.h"
-#include "appl_grid/appl_pdf.h"
-#include "appl_grid/fastnlo.h"
+// PineAPPL includes
+#include "pineappl_capi.h"
 
 #include "fk_utils.h"
 #include "fk_qcd.h"
@@ -24,19 +21,6 @@ using QCD::qcd_param;
 
 namespace PINE
 {
-
-	// // Wrapper for appl::grids (handles fastNLO based grids)
-	class grid
-	{
-	public:
-		grid(std::string const& filename, int const& fnlobin);
-		~grid();
-
-		fastnlo* 	fg;
-		appl::grid* g;
-	};
-
-
 	class SubGrid: public FKSubGrid
 	{
 	public:
@@ -54,18 +38,13 @@ namespace PINE
 
 		// ***********************************************************
 
-		const string 		applfile;	//!< Path for the applgrid file
-		const string 		readme;		//!< APPLgrid README
+		const string 		pineapplfile; //!< Path for the pineapplg file
+		const string 		readme;		  //!< APPLgrid README
 
 		const vector<int> 	maskmap;	//!< Map of masked applgrid points to datapoints
 		const size_t 		ndata;     	//!< Number of selected datapoints
-		const size_t 		ptmin;     	//!< Minimum perturbative order to contribute
-
-		const int    		fnlobin;   	//!< Which bin in the fastNLO grid? (yes, this needs to be signed)
-		const bool 			pdfwgt;     //!< Using a PDF weight? Typically used for fastNLO grids
-		const bool 			ppbar;		//!< Does the grid need a pp -> ppbar transform?
 
-		const grid 			applgrid;	//!< APPLgrid class
+		const pineappl_grid* grid;	//!< PineAPPL class
 
 		static vector<int>  parse_maskmap(string const&);
 	};
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/fk_pine.cc b/src/core/fk_pine.cc
index aa5cb13..9b56f28 100644
--- a/src/core/fk_pine.cc
+++ b/src/core/fk_pine.cc
@@ -80,7 +80,7 @@ namespace PINE
   }
 
   // ********************************* 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.
@@ -204,36 +204,22 @@ namespace PINE
       }
     return nElm;
   }
-
- // *********************************** Grid object **************************************
-
-  grid::grid(std::string const& filename, int const& fnlobin):
-    fg(fnlobin >= 0 ? new fastnlo(filename):0),
-    g(fg ? fg->grids()[fnlobin]:new appl::grid(filename))
-    {}
-
-  grid::~grid() {
-    fg ? delete fg : delete g;
-  }
-
+*/
 
  // *********************************** APPLGrid convolutions **************************************
 
   SubGrid::SubGrid(FKTarget const& parent, NNPDF::IndexDB const& db, int const& iDB):
     FKSubGrid(parent, iDB, NNPDF::dbquery<string>(db, iDB, "operators")),
-    applfile(NNPDF::dbquery<string>(db,iDB,"applgrid")),
+    pineapplfile(NNPDF::dbquery<string>(db,iDB,"pineapplgrid")),
     readme(),
     maskmap(parse_maskmap(NNPDF::dbquery<string>(db,iDB,"mask"))),
-    ndata(	maskmap.size() ),
-    ptmin(	NNPDF::dbquery<size_t>(db,iDB,"ptmin")),
-    fnlobin(NNPDF::dbquery<int>(db,iDB,"fnlobin")),
-    pdfwgt(	NNPDF::dbquery<bool>(db,iDB,"pdfwgt")),
-    ppbar(	NNPDF::dbquery<bool>(db,iDB,"ppbar")),
-    applgrid(applPath() + parent.GetSetName() + "/" + applfile , fnlobin)
+    ndata(maskmap.size()),
+    grid(pineappl_grid_read((pineapplPath() + parent.GetSetName() + "/" + pineapplfile).c_str()))
   {
-    if (ndata>static_cast<size_t>(applgrid.g->Nobs()))
+    const size_t bin_count = pineappl_grid_bin_count(grid);
+    if (ndata > bin_count)
     {
-      cerr <<"Error: number of datapoints in settings: "<<ndata<<" > appl grid Nobs: "<<applgrid.g->Nobs()<<endl;
+      cerr <<"Error: number of datapoints in settings: "<<ndata<<" > appl grid Nobs: "<<bin_count<<endl;
       cerr <<"Please check the provided mask" <<endl;
       exit(-1);
     }
@@ -241,6 +227,7 @@ namespace PINE
 
   void SubGrid::Compute(qcd_param const& par, vector<double>& results) const
   {
+    /*
     if (ppbar == true && par.xiF != 1)
       std::cout << "WARNING: ppbar ROTATION NOT TESTED - APPLgrid does not support fac. scale variation with ppbar so I cannot cross-check" <<std::endl;
     if (par.evol_pto == 2)
@@ -260,6 +247,7 @@ namespace PINE
     for (size_t i=0; i<maskmap.size(); i++)
       for (int const& j : datamap[i])
         results[j] += xsec[maskmap[i]];
+        */
   }
 
 
@@ -267,6 +255,7 @@ namespace PINE
 
   void SubGrid::Combine(QCD::qcd_param const& par, NNPDF::FKGenerator* fk) const
   {
+    /*
     if (par.evol_pto == 2)
       std::cout << "WARNING: APPLgrid does not currently support NNLO convolutions, fixing convolution to NLO" <<std::endl;
 
@@ -397,6 +386,7 @@ namespace PINE
       } // /pto
     } // /data
     return;
+    */
   }
 
   //  *************************************** Metadata methods ********************************************
@@ -405,17 +395,13 @@ namespace PINE
   void SubGrid::Splash(ostream& o) const
   {
     FKSubGrid::Splash(o);
-    o << "- APPLgrid: " << applfile << endl
-      << "- PTMin: "    << ptmin    << endl
-      << "- fnlobin: "  << fnlobin  << endl
-      << "- PDFWeight: "<< pdfwgt   << endl
-      << "- ppbar: "    << ppbar    << endl
-      << endl;
+    o << "- APPLgrid: " << pineapplfile << endl;
   }
 
     // Get the maximum scale of an applgrid
   double SubGrid::GetQ2max() const
   {
+    /*
     // Find maximum required scale
     double Q2max = 0;
     const double nloops = applgrid.g->calculation() == appl::grid::AMCATNLO ? 4 : 2;
@@ -426,12 +412,14 @@ namespace PINE
         Q2max = max(Q2max, igrid->getQ2max());
       }
     return Q2max;
+    */
   }
 
 
   // Return the minimum x used in the subgrid
   double SubGrid::GetXmin() const
   {
+    /*
     const double nloops = applgrid.g->calculation() == appl::grid::AMCATNLO ? 4 : 2;
     double xmin = 1.0;
     for(int i=0; i<nloops; i++)
@@ -457,10 +445,12 @@ namespace PINE
           }
         }
     return xmin;
+    */
   };
 
   double SubGrid::GetComputeXmin() const
   {
+    /*
     const double nloops = applgrid.g->calculation() == appl::grid::AMCATNLO ? 4 : 2;
     double xmin = 1.0;
     for(int i=0; i<nloops; i++)
@@ -474,6 +464,7 @@ namespace PINE
       }
 
     return xmin;
+    */
   };
 
 
diff --git a/src/core/fk_utils.cc b/src/core/fk_utils.cc
index a990e6c..2c10354 100644
--- a/src/core/fk_utils.cc
+++ b/src/core/fk_utils.cc
@@ -39,6 +39,12 @@ std::string exec(const char* cmd) {
     return result;
 }
 
+std::string pineapplPath()
+{
+  std::string pineapplDir(STR(PINEAPPL_PATH));
+  return pineapplDir;
+};
+
 std::string applPath()
 {
   std::string applDir(STR(APPL_PATH));

From ac1d09f81a63e72ff04459730eeb828a31c6bea0 Mon Sep 17 00:00:00 2001
From: Stefano Carrazza <stefano.carrazza@cern.ch>
Date: Wed, 26 Aug 2020 16:21:44 +0200
Subject: [PATCH 03/14] adding prototype for q2max and xmin

---
 src/core/fk_pine.cc | 48 +++++++--------------------------------------
 1 file changed, 7 insertions(+), 41 deletions(-)

diff --git a/src/core/fk_pine.cc b/src/core/fk_pine.cc
index 9b56f28..0c1165f 100644
--- a/src/core/fk_pine.cc
+++ b/src/core/fk_pine.cc
@@ -395,57 +395,23 @@ namespace PINE
   void SubGrid::Splash(ostream& o) const
   {
     FKSubGrid::Splash(o);
-    o << "- APPLgrid: " << pineapplfile << endl;
+    o << "- PineAPPL: " << pineapplfile << endl;
   }
 
     // Get the maximum scale of an applgrid
   double SubGrid::GetQ2max() const
   {
-    /*
-    // Find maximum required scale
-    double Q2max = 0;
-    const double nloops = applgrid.g->calculation() == appl::grid::AMCATNLO ? 4 : 2;
-    for(int i=0; i<nloops; i++)  // pto
-      for (int j=0; j<applgrid.g->Nobs(); j++) // bin
-      {
-        appl::igrid const *igrid = applgrid.g->weightgrid(i, j);
-        Q2max = max(Q2max, igrid->getQ2max());
-      }
-    return Q2max;
-    */
+    vector<double> q2values(pineappl_grid_subgrid_q2_count(grid));
+    pineappl_grid_subgrid_q2(grid, q2values.data());
+    return *std::max_element(q2values.begin(), q2values.end());
   }
 
-
   // Return the minimum x used in the subgrid
   double SubGrid::GetXmin() const
   {
-    /*
-    const double nloops = applgrid.g->calculation() == appl::grid::AMCATNLO ? 4 : 2;
-    double xmin = 1.0;
-    for(int i=0; i<nloops; i++)
-      for (int j=0; j<applgrid.g->Nobs(); j++)
-      {
-        appl::igrid const *igrid = applgrid.g->weightgrid(i, j);
-        const size_t nsubproc = applgrid.g->subProcesses(i);
-
-          for (int t=0; t<igrid->Ntau(); t++) // Loop over Q^2 integral
-          {
-            const std::pair<int,int> l1 = get_igrid_limits_x1(igrid, nsubproc, t);
-            if (l1.first <= l1.second)
-            {
-              xmin = min(xmin, igrid->fx(igrid->gety1(l1.first)));
-              xmin = min(xmin, igrid->fx(igrid->gety1(l1.second)));
-            }
-            const std::pair<int,int> l2 = get_igrid_limits_x2(igrid, nsubproc, t);
-            if (l2.first <= l2.second)
-            {
-              xmin = min(xmin, igrid->fx(igrid->gety2(l2.first)));
-              xmin = min(xmin, igrid->fx(igrid->gety2(l2.second)));
-            }
-          }
-        }
-    return xmin;
-    */
+    vector<double> xvalues(pineappl_grid_subgrid_x_count(grid));
+    pineappl_grid_subgrid_x(grid, xvalues.data());
+    return *std::min_element(xvalues.begin(), xvalues.end());
   };
 
   double SubGrid::GetComputeXmin() const

From 9cf2d089538347834c66e96d7fbcbc11443f1594 Mon Sep 17 00:00:00 2001
From: Stefano Carrazza <stefano.carrazza@cern.ch>
Date: Wed, 26 Aug 2020 16:23:28 +0200
Subject: [PATCH 04/14] removing the annoying exit

---
 src/apfel_comb.cc | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/apfel_comb.cc b/src/apfel_comb.cc
index 345efad..6211c8f 100644
--- a/src/apfel_comb.cc
+++ b/src/apfel_comb.cc
@@ -136,6 +136,5 @@ int main(int argc, char* argv[]) {
 
   cout << "                      APPLComb Complete "<<endl;
   DisplayHR();
-  exit(0);
 }
 

From 87c2472bf261e716a8f083834c1cca94228db4f9 Mon Sep 17 00:00:00 2001
From: Stefano Carrazza <stefano.carrazza@cern.ch>
Date: Thu, 27 Aug 2020 17:00:15 +0200
Subject: [PATCH 05/14] plugin compute and combine

---
 src/apfelcomb/fk_qcd.h |   1 +
 src/core/fk_pine.cc    | 301 ++++++++++++++++++-----------------------
 src/core/fk_qcd.cc     |  10 ++
 3 files changed, 144 insertions(+), 168 deletions(-)

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/core/fk_pine.cc b/src/core/fk_pine.cc
index 0c1165f..8ddfb33 100644
--- a/src/core/fk_pine.cc
+++ b/src/core/fk_pine.cc
@@ -19,22 +19,17 @@ using NNPDF::FKHeader;
 
 namespace PINE
 {
-  #ifdef __APPL_PHOTON
-  const int applgrid_nfl = 14;
-  #else
-  const int applgrid_nfl = 13;
-  #endif
-
 
 // // ********************* Evolution factors ******************************
 
+
   class EvolutionFactors
   {
   public:
     EvolutionFactors(const int nxin, const int nxout):
-    b1(applgrid_nfl),
-    b2(applgrid_nfl*14),
-    b3(applgrid_nfl*14*nxin),
+    b1(14),
+    b2(14*14),
+    b3(14*14*nxin),
     data(new double[nxout*b3])
     {
       for (int i=0; i<nxout*b3; i++)
@@ -54,29 +49,16 @@ namespace PINE
 // ********************************* Basis rotation helpers *************************************
 
   // Rotates APFEL flavour basis into APPLgrid flavour basis (photon moves from 0 to 13)
-  void evolpdf_applgrid(const double& x, const double& Q, double* pdf)
+  double evolpdf_pineapplgrid(int32_t pdgid, double x, double q2, void*)
   {
-    for (int i=0; i<applgrid_nfl; i++) pdf[i]=0;
-    if (x<APFEL::xGrid(0))  // A nice trick of APPLgrid is to request PDF x-values smaller than are actually used
-      return;
-
-    double *APFEL_basis = new double[14]();
-    QCD::evolpdf(x, Q, APFEL_basis);
-
-    for (int i=0; i<13; i++)
-      pdf[i]=APFEL_basis[i+1];
-    if (applgrid_nfl == 14)
-      pdf[13] = APFEL_basis[0];
-
-    delete[] APFEL_basis;
+    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);
   }
 
-  // antiproton version of evolpdf_applgrid, swaps q<->qbar
-  void evolpdf_applgrid_pbar(const double& x, const double& Q, double* pdf)
+  double alphas(double q2, void*)
   {
-    evolpdf_applgrid(x,Q,pdf);
-    for (int i=1; i<7; i++)
-      std::swap(pdf[6+i], pdf[6-i]);
+    return QCD::alphas(sqrt(q2));
   }
 
   // ********************************* Kinematics Helpers *************************************
@@ -227,19 +209,15 @@ namespace PINE
 
   void SubGrid::Compute(qcd_param const& par, vector<double>& results) const
   {
-    /*
-    if (ppbar == true && par.xiF != 1)
-      std::cout << "WARNING: ppbar ROTATION NOT TESTED - APPLgrid does not support fac. scale variation with ppbar so I cannot cross-check" <<std::endl;
     if (par.evol_pto == 2)
       std::cout << "WARNING: APPLgrid does not currently support NNLO convolutions, fixing convolution to NLO" <<std::endl;
 
-    // Compute with applgrid interface
-    const int pto = (ptmin == 1) ? -1: min(par.evol_pto,(size_t)1);
-    vector<double> xsec;
-    if ( ppbar == true && par.xiF == 1)
-      xsec = applgrid.g->vconvolute( evolpdf_applgrid, evolpdf_applgrid_pbar, QCD::alphas, pto, par.xiR, par.xiF );
-    else
-      xsec = applgrid.g->vconvolute( evolpdf_applgrid, QCD::alphas, pto, par.xiR, par.xiF);
+    // TODO: fix pto
+    vector<double> 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;
 
@@ -247,7 +225,6 @@ namespace PINE
     for (size_t i=0; i<maskmap.size(); i++)
       for (int const& j : datamap[i])
         results[j] += xsec[maskmap[i]];
-        */
   }
 
 
@@ -255,138 +232,140 @@ namespace PINE
 
   void SubGrid::Combine(QCD::qcd_param const& par, NNPDF::FKGenerator* fk) const
   {
-    /*
     if (par.evol_pto == 2)
       std::cout << "WARNING: APPLgrid does not currently support NNLO convolutions, fixing convolution to NLO" <<std::endl;
 
-    if (applgrid_nfl == 14)
-      std::cout << "WARNING: Combining with photon channel in APPLgrid" <<std::endl;
-
-    // APPLgrid pointer
-    const appl::grid* g = applgrid.g;
-    const int ptmax = min(par.evol_pto,(size_t)1); // Maximum perturbative order limited to NLO
+    // TODO: fit pto
 
     // Progress monitoring
     int completedElements = 0;
-    const int nXelements = countElements(maskmap, ptmin, ptmax, g);
-    const time_point t1 = std::chrono::system_clock::now();
+    //const int nXelements = countElements(maskmap, ptmin, ptmax, g);
+    //const time_point t1 = std::chrono::system_clock::now();
     const vector<size_t> afl = QCD::active_flavours(par);
 
+    vector<uint32_t> orders(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<double> q2values(pineappl_subgrid_q2_grid_count(grid));
+    pineappl_subgrid_q2_grid(grid, q2values.data());
+
+    vector<double> xvalues(pineappl_subgrid_x_grid_count(grid));
+    pineappl_subgrid_x_grid(grid, xvalues.data());
+    const size_t nxpi = xvalues.size();
+
+    vector<double> bin_sizes(pineappl_grid_bin_count(grid));
+    pineappl_grid_bin_sizes(grid, bin_sizes.data());
+
+    vector<double> weight_matrix(nxpi * nxpi);
+
     for (size_t d=0; d<maskmap.size(); d++)
     {
       // Fetch associated applgrid info
       const size_t bin = maskmap[d];
-      for (size_t pto=0; pto<=(ptmax-ptmin); pto++) // Loop over perturbative order
+
+      for (size_t pto = 0; pto < orders.size()/4; pto++) // Loop over perturbative order
       {
-        // Determine grid index, allocate subprocess arrays
-        const int gidx = get_grid_idx(g, pto+ptmin);
-        appl::appl_pdf *genpdf = get_appl_pdf( g, gidx );
-        const size_t nsubproc = g->subProcesses(gidx);
-        double *W = new double[nsubproc];   // Weight array
-        double *H = new double[nsubproc];   // Evolution factor array
-        double *H1 = new double[nsubproc];  // Split evolution factor array 1
-        double *H2 = new double[nsubproc];  // Split evolution factor array 2
-
-        // Fetch grid pointer and loop over Q
-        appl::igrid const *igrid = g->weightgrid(gidx, bin);
-        for (int t=0; t<igrid->Ntau(); t++)
-        {
-          // Scales and strong coupling
-          const double Q2  = igrid->fQ2( igrid->gettau(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 fA1(nxin, igrid->Ny1());   // PDF 1 Evolution factors
-          EvolutionFactors fA2(nxin, igrid->Ny2());   // PDF 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
-          const std::pair<int,int> l1 = get_igrid_limits_x1(igrid, nsubproc, t);
-          const std::pair<int,int> l2 = get_igrid_limits_x2(igrid, nsubproc, t);
-
-          for (int ix = 0; ix < nxin; ix++)
-          for (size_t fl : afl)
-          {
-            for (int ox=l1.first; ox<=l1.second; ox++) // Loop over applgrid x1
-            {
-              QCD::EvolutionOperator(false,applgrid_nfl==14,ix,igrid->fx(igrid->gety1(ox)),fl,QF,fA1(ox,ix,fl));
-              if (vary_fac) QCD::DerivativeOperator(false,applgrid_nfl==14,ix,igrid->fx(igrid->gety1(ox)),fl,QF,fdA1(ox,ix,fl));
-            }
-            for (int ox=l2.first; ox<=l2.second; ox++) // Loop over applgrid x2
-            {
-              QCD::EvolutionOperator(ppbar,applgrid_nfl==14,ix,igrid->fx(igrid->gety2(ox)),fl,QF,fA2(ox,ix,fl));
-              if (vary_fac) QCD::DerivativeOperator(ppbar,applgrid_nfl==14,ix,igrid->fx(igrid->gety2(ox)),fl,QF,fdA2(ox,ix,fl));
-            }
-          }
+        // TODO: unSkip renormalization and factorization scale variation.
+        if (orders.at(4 * pto + 2) > 0 || orders.at(4 * pto + 3) > 0)
+          continue;
 
-          for (int a=l1.first; a<=l1.second; a++ )
+        for (size_t lumi = 0; lumi < lumi_count; lumi++)
+        {
+          // Fetch grid pointer and loop over Q
+          size_t slice_indices[2];
+          pineappl_subgrid_filled_q2_slices(grid, pto, bin, lumi, slice_indices);
+          for (size_t t = slice_indices[0]; t < slice_indices[1]; t++)
           {
-            const double x1 = igrid->fx(igrid->gety1(a));
-            const std::pair<int,int> limits = get_slice_limits(igrid, nsubproc, t, a);
-            for (int b=limits.first; b<=limits.second; b++) // Loop over applgrid x2
-            {
-              // fetch weight values
-              for (size_t ip=0; ip<nsubproc; ip++)
-                W[ip] = (*(const SparseMatrix3d*) const_cast<appl::igrid*>(igrid)->weightgrid(ip))(t,a,b);
-
-              // Calculate normalisation factors
-              const double x2 = igrid->fx(igrid->gety2(b));
-              const double pdfnrm = pdfwgt ? igrid->weightfun(x1)*igrid->weightfun(x2) : 1.0;
-              const double norm = pdfnrm*compute_wgt_norm(g, bin, pto+ptmin, as, x1, x2)*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
+            // 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[a * nxpi + b];
+
+                // Calculate normalisation factors
+                const double norm = pow(as, orders.at(4 * pto + 0))/(xvalues.at(a)*xvalues.at(b)*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
+                  {
+                    // Rotate to subprocess basis
+                    const size_t combinations = pineappl_lumi_combinations(grid_lumi, lumi);
+                    vector<int32_t> pdgids(2 * combinations);
+                    vector<double> factors(combinations);
+                    pineappl_lumi_entry(grid_lumi, lumi, pdgids.data(), factors.data());
+
+                    if (W != 0)
                     {
-                      // Rotate to subprocess basis
-                      genpdf->evaluate(fA1(a,i,k),fA2(b,j,l),H);
-
-                      if (vary_fac)
+                      for (size_t m = 0; m < combinations; m++)
                       {
-                        genpdf->evaluate(fdA1(a,i,k),fA2(b,j,l),H1);
-                        genpdf->evaluate(fA1(a,i,k),fdA2(b,j,l),H2);
-                      }
-
-                      for (size_t ip=0; ip<nsubproc; ip++)
-                        if (W[ip] != 0)
-                        {
-                          double fill = H[ip];                              // 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[ip]);
+                        const int32_t k = pdgids.at(2 * m + 0);
+                        const int32_t l = pdgids.at(2 * m + 1);
+
+                        const double H = factors.at(m) * (*fA(a, i, k)) * (*fA(b, j, l));
+
+                        // 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);
-            }
+                // Update progress
+                completedElements++;
+                //StatusUpdate(t1, (double)completedElements/(double)nXelements, cout);
+              }
           }
         }
-
-        // Free subprocess arrays
-        delete[] W;
-        delete[] H;
-        delete[] H1;
-        delete[] H2;
-
       } // /pto
     } // /data
+
+    pineappl_lumi_delete(grid_lumi);
+
     return;
-    */
   }
 
   //  *************************************** Metadata methods ********************************************
@@ -401,39 +380,25 @@ namespace PINE
     // Get the maximum scale of an applgrid
   double SubGrid::GetQ2max() const
   {
-    vector<double> q2values(pineappl_grid_subgrid_q2_count(grid));
-    pineappl_grid_subgrid_q2(grid, q2values.data());
-    return *std::max_element(q2values.begin(), q2values.end());
+    vector<double> 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
   {
-    vector<double> xvalues(pineappl_grid_subgrid_x_count(grid));
-    pineappl_grid_subgrid_x(grid, xvalues.data());
-    return *std::min_element(xvalues.begin(), xvalues.end());
+    // TODO: check this implementation.
+    return GetComputeXmin();
   };
 
   double SubGrid::GetComputeXmin() const
   {
-    /*
-    const double nloops = applgrid.g->calculation() == appl::grid::AMCATNLO ? 4 : 2;
-    double xmin = 1.0;
-    for(int i=0; i<nloops; i++)
-      for (int j=0; j<applgrid.g->Nobs(); j++)
-      {
-        appl::igrid const *igrid = applgrid.g->weightgrid(i, j);
-        xmin = min(xmin, igrid->fx(igrid->gety1(0)));
-        xmin = min(xmin, igrid->fx(igrid->gety1(igrid->Ny1()-1)));
-        xmin = min(xmin, igrid->fx(igrid->gety2(0)));
-        xmin = min(xmin, igrid->fx(igrid->gety2(igrid->Ny2()-1)));
-      }
-
-    return xmin;
-    */
+    vector<double> xvalues(pineappl_subgrid_x_grid_count(grid));
+    pineappl_subgrid_x_grid(grid, xvalues.data());
+    return std::min(xvalues.front(), xvalues.back());
   };
 
-
   vector<int> SubGrid::parse_maskmap(string const& mask)
   {
     const vector<string> masksplit = ssplit(mask);
diff --git a/src/core/fk_qcd.cc b/src/core/fk_qcd.cc
index f4fba02..8fd4123 100644
--- a/src/core/fk_qcd.cc
+++ b/src/core/fk_qcd.cc
@@ -285,6 +285,16 @@ namespace QCD
       pdf[i+7]= APFEL::xPDF(i,fx);
   }
 
+  // APFEL PDF in the evolution basis
+  double flvpdf(int32_t pdgid, double x, double q2)
+  {
+    updateEvol(sqrt(q2));
+    const double fx = fabs(x-APFEL::xGrid(0))/x < 1E-6 ? APFEL::xGrid(0):x;  // Fuzzy x
+    if (pdgid == 22)
+      return APFEL::xgamma(fx);
+    return APFEL::xPDF(pdgid, fx);
+  }
+
   // APFEL strong coupling
   double alphas(const double& Q)
   {

From 5922567ffa33750d8e4a6f61db001df9da32e13d Mon Sep 17 00:00:00 2001
From: Stefano Carrazza <stefano.carrazza@cern.ch>
Date: Wed, 2 Sep 2020 16:08:41 +0200
Subject: [PATCH 06/14] adding changes

---
 src/core/fk_pine.cc | 62 ++++++++++++++++++++++++++++++++-------------
 1 file changed, 44 insertions(+), 18 deletions(-)

diff --git a/src/core/fk_pine.cc b/src/core/fk_pine.cc
index 8ddfb33..7d99a41 100644
--- a/src/core/fk_pine.cc
+++ b/src/core/fk_pine.cc
@@ -13,6 +13,7 @@
 #include <chrono>
 #include <ctime>
 #include <limits>
+#include <stdexcept>
 
 using namespace std;
 using NNPDF::FKHeader;
@@ -30,20 +31,16 @@ namespace PINE
     b1(14),
     b2(14*14),
     b3(14*14*nxin),
-    data(new double[nxout*b3])
-    {
-      for (int i=0; i<nxout*b3; i++)
-        data[i] = 0;
-    };
-    ~EvolutionFactors() {delete[] data;};
+    data(nxout*b3)
+    {}
 
-    double* operator()(int const& ox, int const& ix, int const& fi ) {return data+b3*ox+b2*ix+b1*fi;};
-    const double* operator()(int const& ox, int const& ix, int const& fi ) const {return data+b3*ox+b2*ix+b1*fi;};
+    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 int b1;
-    const int b2;
-    const int b3;
-    double* data;
+    const size_t b1;
+    const size_t b2;
+    const size_t b3;
+    std::vector<double> data;
   };
 
 // ********************************* Basis rotation helpers *************************************
@@ -61,6 +58,36 @@ namespace PINE
     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 14;
+
+    case 21:
+      return 0;
+
+    case -6:
+    case -5:
+    case -4:
+    case -3:
+    case -2:
+    case -1:
+    case 6:
+    case 5:
+    case 4:
+    case 3:
+    case 2:
+    case 1:
+      return pdgid + 7;
+
+    default:
+      throw std::runtime_error("pdgid not available.");
+    }
+  }
+
   // ********************************* Kinematics Helpers *************************************
 /*
   // Returns the minimum and maximum x-grid points for a specified subgrid slice.
@@ -318,7 +345,7 @@ namespace PINE
                 const double W = weight_matrix[a * nxpi + b];
 
                 // Calculate normalisation factors
-                const double norm = pow(as, orders.at(4 * pto + 0))/(xvalues.at(a)*xvalues.at(b)*bin_sizes.at(bin))*nrmdat;
+                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
@@ -329,12 +356,11 @@ namespace PINE
                     vector<double> factors(combinations);
                     pineappl_lumi_entry(grid_lumi, lumi, pdgids.data(), factors.data());
 
-                    if (W != 0)
-                    {
-                      for (size_t m = 0; m < combinations; m++)
+                    for (size_t k : afl)         // loop over flavour 1
+                      for (size_t l : afl)       // loop over flavour 2
                       {
-                        const int32_t k = pdgids.at(2 * m + 0);
-                        const int32_t l = pdgids.at(2 * m + 1);
+                        //const int32_t k = pdgid_to_apfel_array_id(pdgids.at(2 * m + 0));
+                        //const int32_t l = pdgid_to_apfel_array_id(pdgids.at(2 * m + 1));
 
                         const double H = factors.at(m) * (*fA(a, i, k)) * (*fA(b, j, l));
 

From a8b0c71e111651be2d86cf667173c82db70b8691 Mon Sep 17 00:00:00 2001
From: Christopher Schwan <handgranaten-herbert@posteo.de>
Date: Thu, 3 Sep 2020 18:36:12 +0200
Subject: [PATCH 07/14] Fix two mistakes

---
 Makefile            |  2 +-
 src/core/fk_pine.cc | 12 +++++++-----
 2 files changed, 8 insertions(+), 6 deletions(-)

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/src/core/fk_pine.cc b/src/core/fk_pine.cc
index 7d99a41..33b78af 100644
--- a/src/core/fk_pine.cc
+++ b/src/core/fk_pine.cc
@@ -359,10 +359,13 @@ namespace PINE
                     for (size_t k : afl)         // loop over flavour 1
                       for (size_t l : afl)       // loop over flavour 2
                       {
-                        //const int32_t k = pdgid_to_apfel_array_id(pdgids.at(2 * m + 0));
-                        //const int32_t l = pdgid_to_apfel_array_id(pdgids.at(2 * m + 1));
-
-                        const double H = factors.at(m) * (*fA(a, i, k)) * (*fA(b, j, l));
+                        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)
                         // {
@@ -376,7 +379,6 @@ namespace PINE
                         if ( fill != 0.0 )
                           for (int const& td : datamap[d])
                             fk->Fill(td, i, j, k, l, norm*fill*W);
-                        }
                     }
                   }
 

From c03c5f03c2e8babfcd078dfffcbb0bf650198ff4 Mon Sep 17 00:00:00 2001
From: Stefano Carrazza <stefano.carrazza@cern.ch>
Date: Thu, 3 Sep 2020 20:26:40 +0200
Subject: [PATCH 08/14] adjustments

---
 src/core/fk_pine.cc | 48 ++++++++++++++++++++++++++++++---------------
 1 file changed, 32 insertions(+), 16 deletions(-)

diff --git a/src/core/fk_pine.cc b/src/core/fk_pine.cc
index 33b78af..f927698 100644
--- a/src/core/fk_pine.cc
+++ b/src/core/fk_pine.cc
@@ -64,10 +64,10 @@ namespace PINE
     switch (pdgid)
     {
     case 22:
-      return 14;
+      return 13;
 
     case 21:
-      return 0;
+      return 6;
 
     case -6:
     case -5:
@@ -81,7 +81,7 @@ namespace PINE
     case 3:
     case 2:
     case 1:
-      return pdgid + 7;
+      return pdgid + 6;
 
     default:
       throw std::runtime_error("pdgid not available.");
@@ -263,11 +263,7 @@ namespace PINE
       std::cout << "WARNING: APPLgrid does not currently support NNLO convolutions, fixing convolution to NLO" <<std::endl;
 
     // TODO: fit pto
-
     // Progress monitoring
-    int completedElements = 0;
-    //const int nXelements = countElements(maskmap, ptmin, ptmax, g);
-    //const time_point t1 = std::chrono::system_clock::now();
     const vector<size_t> afl = QCD::active_flavours(par);
 
     vector<uint32_t> orders(pineappl_grid_order_count(grid));
@@ -288,6 +284,28 @@ namespace PINE
 
     vector<double> weight_matrix(nxpi * nxpi);
 
+    int completedElements = 0;
+    int nXelements = 0;
+    for (size_t d=0; d<maskmap.size(); d++)
+    {
+      const size_t bin = maskmap[d];
+      for (size_t pto = 0; pto < orders.size()/4; pto++) // Loop over perturbative order
+      {
+        // TODO: unSkip renormalization and factorization scale variation.
+        if (orders.at(4 * pto + 2) > 0 || orders.at(4 * pto + 3) > 0)
+          continue;
+
+        for (size_t lumi = 0; lumi < lumi_count; lumi++)
+        {
+          size_t slice_indices[2];
+          pineappl_subgrid_filled_q2_slices(grid, pto, bin, lumi, slice_indices);
+          nXelements += (slice_indices[1] - slice_indices[0]) * nxpi * nxpi;
+        }
+      }
+    }
+
+    const time_point t1 = std::chrono::system_clock::now();
+
     for (size_t d=0; d<maskmap.size(); d++)
     {
       // Fetch associated applgrid info
@@ -301,6 +319,12 @@ namespace PINE
 
         for (size_t lumi = 0; lumi < lumi_count; lumi++)
         {
+          // prepare luminosity combinations
+          const size_t combinations = pineappl_lumi_combinations(grid_lumi, lumi);
+          vector<int32_t> pdgids(2 * combinations);
+          vector<double> factors(combinations);
+          pineappl_lumi_entry(grid_lumi, lumi, pdgids.data(), factors.data());
+
           // Fetch grid pointer and loop over Q
           size_t slice_indices[2];
           pineappl_subgrid_filled_q2_slices(grid, pto, bin, lumi, slice_indices);
@@ -349,13 +373,6 @@ namespace PINE
 
                 for (int i = 0; i < nxin; i++)    // Loop over input pdf x1
                   for (int j = 0; j < nxin; j++)  // Loop over input pdf x2
-                  {
-                    // Rotate to subprocess basis
-                    const size_t combinations = pineappl_lumi_combinations(grid_lumi, lumi);
-                    vector<int32_t> pdgids(2 * combinations);
-                    vector<double> factors(combinations);
-                    pineappl_lumi_entry(grid_lumi, lumi, pdgids.data(), factors.data());
-
                     for (size_t k : afl)         // loop over flavour 1
                       for (size_t l : afl)       // loop over flavour 2
                       {
@@ -380,11 +397,10 @@ namespace PINE
                           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);
+                StatusUpdate(t1, (double)completedElements/(double)nXelements, cout);
               }
           }
         }

From 08bcde24683fceada0fe76e0bfa0e310bf51f5c2 Mon Sep 17 00:00:00 2001
From: Stefano Carrazza <stefano.carrazza@cern.ch>
Date: Fri, 4 Sep 2020 22:24:45 +0200
Subject: [PATCH 09/14] updating nx points

---
 db/apfelcomb.dat | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/db/apfelcomb.dat b/db/apfelcomb.dat
index c73dafd..4a43cd4 100644
--- a/db/apfelcomb.dat
+++ b/db/apfelcomb.dat
@@ -614,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');
@@ -916,7 +916,7 @@ CREATE TABLE IF NOT EXISTS "fixmegrids" (
 	`positivity`	INTEGER NOT NULL,
 	`source`	TEXT NOT NULL
 );
-INSERT INTO fixmegrids VALUES(0,'ATLASZHIGHMASS49FB','ATLASZHIGHMASS49FB','ATLAS 49FB High mass Z',20,0,'PINE');
+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);

From aa7e2889115c5822f40f9760c43df4d1180653ba Mon Sep 17 00:00:00 2001
From: Christopher Schwan <handgranaten-herbert@posteo.de>
Date: Mon, 7 Sep 2020 10:05:57 +0200
Subject: [PATCH 10/14] Replace array access operator with `at`

---
 src/core/fk_pine.cc | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/src/core/fk_pine.cc b/src/core/fk_pine.cc
index f927698..bf2d7d7 100644
--- a/src/core/fk_pine.cc
+++ b/src/core/fk_pine.cc
@@ -10,6 +10,7 @@
 #include "NNPDF/nnpdfdb.h"
 #include "NNPDF/commondata.h"
 
+#include <array>
 #include <chrono>
 #include <ctime>
 #include <limits>
@@ -251,7 +252,7 @@ namespace PINE
     // Relate back to results
     for (size_t i=0; i<maskmap.size(); i++)
       for (int const& j : datamap[i])
-        results[j] += xsec[maskmap[i]];
+        results.at(j) += xsec.at(maskmap[i]);
   }
 
 
@@ -297,9 +298,9 @@ namespace PINE
 
         for (size_t lumi = 0; lumi < lumi_count; lumi++)
         {
-          size_t slice_indices[2];
-          pineappl_subgrid_filled_q2_slices(grid, pto, bin, lumi, slice_indices);
-          nXelements += (slice_indices[1] - slice_indices[0]) * nxpi * nxpi;
+          std::array<size_t, 2> 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;
         }
       }
     }
@@ -326,9 +327,9 @@ namespace PINE
           pineappl_lumi_entry(grid_lumi, lumi, pdgids.data(), factors.data());
 
           // Fetch grid pointer and loop over Q
-          size_t slice_indices[2];
-          pineappl_subgrid_filled_q2_slices(grid, pto, bin, lumi, slice_indices);
-          for (size_t t = slice_indices[0]; t < slice_indices[1]; t++)
+          std::array<size_t, 2> 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);
@@ -366,7 +367,7 @@ namespace PINE
               for (size_t b = 0; b < nxpi; b++) // Loop over applgrid x2
               {
                 // fetch weight values
-                const double W = weight_matrix[a * nxpi + b];
+                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;

From c0d3366d3120b885e92e5d197d592aefcc99fec5 Mon Sep 17 00:00:00 2001
From: Christopher Schwan <handgranaten-herbert@posteo.de>
Date: Mon, 7 Sep 2020 10:08:50 +0200
Subject: [PATCH 11/14] Reduce heap allocations

---
 src/core/fk_pine.cc | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/core/fk_pine.cc b/src/core/fk_pine.cc
index bf2d7d7..460f8f2 100644
--- a/src/core/fk_pine.cc
+++ b/src/core/fk_pine.cc
@@ -285,6 +285,9 @@ namespace PINE
 
     vector<double> weight_matrix(nxpi * nxpi);
 
+    vector<int32_t> pdgids;
+    vector<double> factors;
+
     int completedElements = 0;
     int nXelements = 0;
     for (size_t d=0; d<maskmap.size(); d++)
@@ -322,8 +325,8 @@ namespace PINE
         {
           // prepare luminosity combinations
           const size_t combinations = pineappl_lumi_combinations(grid_lumi, lumi);
-          vector<int32_t> pdgids(2 * combinations);
-          vector<double> factors(combinations);
+          pdgids.resize(2 * combinations);
+          factors.resize(combinations);
           pineappl_lumi_entry(grid_lumi, lumi, pdgids.data(), factors.data());
 
           // Fetch grid pointer and loop over Q

From 6100be5b2311b48db2384eff72a824e0e384eebf Mon Sep 17 00:00:00 2001
From: Christopher Schwan <handgranaten-herbert@posteo.de>
Date: Mon, 7 Sep 2020 10:47:37 +0200
Subject: [PATCH 12/14] Fix memory bug

---
 src/core/fk_pine.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/core/fk_pine.cc b/src/core/fk_pine.cc
index 460f8f2..33f0742 100644
--- a/src/core/fk_pine.cc
+++ b/src/core/fk_pine.cc
@@ -267,7 +267,7 @@ namespace PINE
     // Progress monitoring
     const vector<size_t> afl = QCD::active_flavours(par);
 
-    vector<uint32_t> orders(pineappl_grid_order_count(grid));
+    vector<uint32_t> orders(4 * pineappl_grid_order_count(grid));
     pineappl_grid_order_params(grid, orders.data());
 
     auto *grid_lumi = pineappl_grid_lumi(grid);

From 899d7c73bdae32cbc3de198b19c7905346c8dfa6 Mon Sep 17 00:00:00 2001
From: Christopher Schwan <handgranaten-herbert@posteo.de>
Date: Mon, 7 Sep 2020 16:11:09 +0200
Subject: [PATCH 13/14] Improve exception message

---
 src/core/fk_pine.cc | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/core/fk_pine.cc b/src/core/fk_pine.cc
index 33f0742..0931ebf 100644
--- a/src/core/fk_pine.cc
+++ b/src/core/fk_pine.cc
@@ -15,6 +15,7 @@
 #include <ctime>
 #include <limits>
 #include <stdexcept>
+#include <string>
 
 using namespace std;
 using NNPDF::FKHeader;
@@ -85,7 +86,7 @@ namespace PINE
       return pdgid + 6;
 
     default:
-      throw std::runtime_error("pdgid not available.");
+      throw std::runtime_error("pdgid `" + std::to_string(pdgid) + "` not available.");
     }
   }
 

From 7122e6c217848a9effdd3e91adcce50f7e537292 Mon Sep 17 00:00:00 2001
From: Christopher Schwan <handgranaten-herbert@posteo.de>
Date: Mon, 7 Sep 2020 16:12:14 +0200
Subject: [PATCH 14/14] Add gluon

---
 src/core/fk_pine.cc | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/core/fk_pine.cc b/src/core/fk_pine.cc
index 0931ebf..4913fca 100644
--- a/src/core/fk_pine.cc
+++ b/src/core/fk_pine.cc
@@ -77,6 +77,7 @@ namespace PINE
     case -3:
     case -2:
     case -1:
+    case 0:
     case 6:
     case 5:
     case 4: