diff --git a/Reactions.txt b/Reactions.txt
index 2fd18edaa4ee7508bc2a81b82e18beb17c7868f0..a3de58155e229b94734b0dfc06fae86eeecbadb3 100644
--- a/Reactions.txt
+++ b/Reactions.txt
@@ -17,3 +17,4 @@ FONLL_DISCC libfonll_discc_xfitter.so
 BaseEvolution libbaseevolution_xfitter.so
 APFELxx libapfelxx_xfitter.so
 UvDvubardbars libuvdvubardbars_xfitter.so
+GRV_PionPdfDecomposition libgrv_pionpdfdecomposition_xfitter.so
diff --git a/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h b/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h
new file mode 100644
index 0000000000000000000000000000000000000000..465644ff1cc2725821c3508752d433b5865819c7
--- /dev/null
+++ b/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h
@@ -0,0 +1,39 @@
+#pragma once
+#include "BasePdfDecomposition.h"
+/**
+  @class GRV_PionPdfDecomposition 
+
+  @brief A class for GRV_Pion pdf decomposition
+
+  Used for pi-
+  Assumes that at starting scale:
+    ubar=d
+    dbar=u
+    s=sbar=0
+  Parametrised distributions are:
+    v   := dval-uval
+    qbar:=(ubar+dbar)/2
+    g
+  Therefore, transformations to physical basis:
+    u=dbar=qbar-v/4
+    d=ubar=qbar+v/4
+    g=g
+    others=0
+  And sum rules for pi- are:
+    \int_0^1 v dx=2
+    \int_0^1 x*(4*qbar+g) dx=1
+  @version 0.2
+  @date 2018-08-14
+  */
+namespace xfitter{
+class GRV_PionPdfDecomposition:public BasePdfDecomposition{
+	public:
+		GRV_PionPdfDecomposition();
+		GRV_PionPdfDecomposition(const std::string& inName);
+		~GRV_PionPdfDecomposition();
+		virtual void initAtStart(const std::string & pars) override final;
+		virtual std::function<std::map<int,double>(const double& x)>f0()const override final; 
+	private:
+		BasePdfParam*par_v,*par_qbar,*par_g;
+	};
+}
diff --git a/pdfdecompositions/GRV_PionPdfDecomposition/src/GRV_PionPdfDecomposition.cc b/pdfdecompositions/GRV_PionPdfDecomposition/src/GRV_PionPdfDecomposition.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8cad3faa5d44c1c559b32d840025275454314f80
--- /dev/null
+++ b/pdfdecompositions/GRV_PionPdfDecomposition/src/GRV_PionPdfDecomposition.cc
@@ -0,0 +1,98 @@
+/*
+   @file GRV_PionPdfDecomposition.cc
+   @date 2018-08-07
+   @author  AddPdfDecomposition.py
+   Created by  AddPdfDecomposition.py on GRV_Pion
+*/
+#include"GRV_PionPdfDecomposition.h"
+using uint=unsigned int;
+using namespace std;
+void initParameterisation(BasePdfParam*prm,const string&decompositionName){
+	const uint N=prm->getNPar();
+	prm->pars=new double*[N];
+	const string prefix=decompositionName+"_"+prm->getName()+"_p";
+	for(uint i=0;i<N;++i){
+		prm->pars[i]=XFITTER_PARS::gParameters[prefix+to_string(i)];
+	}
+}
+namespace xfitter{
+//For dynamic loading:
+extern "C" GRV_PionPdfDecomposition*create(){
+	return new GRV_PionPdfDecomposition();
+}
+GRV_PionPdfDecomposition::GRV_PionPdfDecomposition():BasePdfDecomposition("GRV_Pion"){
+	par_v   =nullptr;
+	par_qbar=nullptr;
+	par_g   =nullptr;
+}
+GRV_PionPdfDecomposition::GRV_PionPdfDecomposition(const std::string& inName):BasePdfDecomposition(inName){}
+GRV_PionPdfDecomposition::~GRV_PionPdfDecomposition(){
+	if(par_v   )delete par_v;
+	if(par_qbar)delete par_qbar;
+	if(par_g   )delete par_g;
+}
+// Init at start:
+void GRV_PionPdfDecomposition::initAtStart(const std::string & pars){
+	//HARDCODE copied from UvDvUbarDbarS, then modified
+	//TO BE REPLACED LATER
+	for(const auto&node:XFITTER_PARS::gParametersY.at("GRV_PionPdfDecomposition")["HERAPDF_pdfparam"]){
+		const string prmzName=node.first.as<string>();//Name of parameterisation
+		BasePdfParam*pParam=new HERAPDF_PdfParam(prmzName);
+		double*parValues=pParam->initFromYaml(node.second);//Parameters of the parameterisation
+		//Register parameters of this parameterisation in the global parameter map
+		for(int i=0;i<pParam->getNPar();i++){
+			double val     =parValues[i];
+			double step    =std::fabs(val)/100.;       /// if 0, parameter is fixed !!! 
+			double minv    =0;//0 means no bounds for minuit
+			double maxv    =0;
+			double priorVal=0;
+			double priorUnc=0;
+			int    add     =true;
+			std::ostringstream ss;
+			ss<<getName()<<'_'<<prmzName<<"_p"<<i;
+			const sring pnam=ss.str();
+			std::cout<<"INFO[GRV_PionPdfDecomposition]: Registering parameter "<<ss<<"="<<val<<std::endl;
+			/// Here it goes to minuit:
+			addexternalparam_(pnam.c_str(),val,step,minv,maxv,priorVal,priorUnc,add,&XFITTER_PARS::gParameters,pnam.size() );
+		}
+    addParameterisation(pdfName,pParam);
+  }
+//The following is not very nice: Decomposition should not create or initialize parameterisations
+	par_v   =initParameterisation(getPdfParam("v   "),getName());
+	par_qbar=initParameterisation(getPdfParam("qbar"),getName());
+	par_g   =initParameterisation(getPdfParam("g   "),getName());
+}
+void GRV_PionPdfDecomposition:initAtIteration() {
+	//Enforce sum rules
+	//Valence sum
+	par_v->setMoment(-1,2);
+	//Momentum sum
+	par_g->setMoment(0,1-4*par_qbar->moment(0));
+}
+// Returns a LHAPDF-style function, that returns PDFs in a physical basis for given x
+std::function<std::map<int,double>(const double& x)>GRV_PionPdfDecomposition::f0()const{
+	return [=](double const& x)->std::map<int, double>{
+		double v   =(*par_v)(x);
+		double qbar=(*par_qbar)(x);
+		double g   =(*par_g)(x);
+		double u=qbar-v/4;
+		double d=qbar+v/4;
+		std::map<int,double>res_={
+			{-6,0},
+			{-5,0},
+			{-4,0},
+			{-3,0},
+			{-2,d},//ubar
+			{-1,u},//dbar
+			{ 1,d},
+			{ 2,u},
+			{ 3,0},
+			{ 4,0},
+			{ 5,0},
+			{ 6,0},
+			{21,g}
+		};
+		return res_;
+	};
+}
+}
diff --git a/pdfdecompositions/GRV_PionPdfDecomposition/src/Makefile.am b/pdfdecompositions/GRV_PionPdfDecomposition/src/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..49ee9e52af54f000e966a81084ed8d4d623985c2
--- /dev/null
+++ b/pdfdecompositions/GRV_PionPdfDecomposition/src/Makefile.am
@@ -0,0 +1,12 @@
+
+# Created by AddPdfDecomposition.py on 2018-08-07
+
+AM_CXXFLAGS = -I$(srcdir)/../include  -I$(srcdir)/../../../include  -I$(srcdir)/../../../pdfparams/BasePdfParam/include/   -I$(srcdir)/../../BasePdfDecomposition/include -Wall -fPIC -Wno-deprecated 
+
+lib_LTLIBRARIES = libGRV_PionPdfDecomposition_xfitter.la
+libGRV_PionPdfDecomposition_xfitter_la_SOURCES = GRV_PionPdfDecomposition.cc
+
+datadir = ${prefix}/yaml/pdfdecompositions/GRV_Pion
+data_DATA = ../yaml/parameters.yaml
+
+dist_noinst_HEADERS = ../include ../yaml
diff --git a/pdfdecompositions/GRV_PionPdfDecomposition/yaml/parameters.yaml b/pdfdecompositions/GRV_PionPdfDecomposition/yaml/parameters.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/pdfdecompositions/UvDvUbarDbarS/include/UvDvUbarDbarS.h b/pdfdecompositions/UvDvUbarDbarS/include/UvDvUbarDbarS.h
index 89e55a10b07afc39672ea50b9d5aaa4867c1c362..74e06bcc22460d09d7ba99308483211d22b78f1a 100644
--- a/pdfdecompositions/UvDvUbarDbarS/include/UvDvUbarDbarS.h
+++ b/pdfdecompositions/UvDvUbarDbarS/include/UvDvUbarDbarS.h
@@ -40,7 +40,7 @@ namespace xfitter
     std::unique_ptr<double[]> getParValues(BasePdfParam const* param) const;
 
     /// Get valence
-    double valence(double x, std::string const& name, double sum) const;
+    double valence(double x,const std::string&name)const;
 
     /// Get sea
     double sea(double x, std::string const& name) const;
@@ -62,10 +62,5 @@ namespace xfitter
 
     /// Get g
     double g(double x) const;
-
-  private:
-    /// sum-rule fixed parameters
-    double _uSum,  _dSum, _gSum;
-    
   };
 }
diff --git a/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc b/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc
index d799d85911fb69af45075ef4777c7a801c7ac045..6ca48015e4d6d9414e4609e9454fe1c7592cef04 100644
--- a/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc
+++ b/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc
@@ -1,7 +1,7 @@
  
 /*
    @file UvDvUbarDbarS_PdfDecomposition.cc
-   @date 2018-07-11
+   @date 2018-08-14
    @author  AddPdfDecomposition.py
    Created by  AddPdfDecomposition.py on 2018-07-11
 */
@@ -13,19 +13,6 @@
 #include <iomanip>
 #include <cmath>
 
-/// TEMPORARY XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
-extern "C" {
-  /// Interface to minuit parameters
-  void addexternalparam_(const char name[],  const double &val, 
-                         const double  &step,
-                         const double &min, const double &max, 
-                         const double &prior, const double &priorUnc,
-                         const int &add, 
-                         map<std::string,double*> *map,
-                         int len);
-}
-
-
 namespace xfitter
 {
   /// the class factories
@@ -38,83 +25,35 @@ namespace xfitter
   UvDvUbarDbarS::UvDvUbarDbarS(): BasePdfDecomposition{"UvDvUbarDbarS"} { }
 
   //_________________________________________________________________________________
-  void UvDvUbarDbarS::initAtStart(const std::string & pars)
-  {
+  void UvDvUbarDbarS::initAtStart(const std::string & pars){
     // HARDWIRE A BIT FOR NOW XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
+    // Currently we create and initialise parameterisations
+    // This is not great, it is not Decomposition's job, they should be created globally
     const YAML::Node pdfs    = XFITTER_PARS::gParametersY.at("UvDvubardbars")["HERAPDF_pdfparam"];
 
     for (const auto& pdf : pdfs) {
       const string pdfName = pdf.first.as<string>();
       BasePdfParam* pParam = new HERAPDF_PdfParam(pdfName);
-      // read its parameters:
-      double *parValues = pParam->initFromYaml(pdf.second);
-
-
-      /// HARDWIRE old-way for now:  XXXXXXXXXXXXXXXXXXXXXXXXXXXXX
-      for (int i = 0; i< pParam->getNPar(); i++) {
-	// get a good name
-	const std::string pnam = getName() + "_"+ pdfName + "_p" +  std::to_string(i) ;
-	std::cout << pnam ;
-	std::cout << " Pars: " << parValues[i] <<std::endl;
-
-	double val = parValues[i];
-	double step = std::fabs(val)/100.;       /// if 0, parameter is fixed !!! 
-
-	
-	double minv  = 0;
-	double maxv  = 0;
-	double priorVal = 0;
-	double priorUnc = 0;
-	int add = true;
-
-	
-	/// Here it goes to minuit:
-	addexternalparam_(pnam.c_str(), val, step, minv, maxv, priorVal, priorUnc, add, &XFITTER_PARS::gParameters,  pnam.size() );
-      }
-      
-      addParameterisation(pdfName,pParam);
+      pParam->initFromYaml(pdf.second);
+			addParameterisation(pdfName,pParam);
     }
   }
 
   void UvDvUbarDbarS::initAtIteration() {
+    //Enforce sum rules
     // counting sum-rules for uv and dv
-    const BasePdfParam* pdfParam = getPdfParam("xuv");
-    std::unique_ptr<double[]> pars = getParValues(pdfParam);  
-    pars[0] = 1.;
-    double sum0 = pdfParam->moment(pars.get(),-1);
-    _uSum = 2.0 / sum0;
-
-    pdfParam = getPdfParam("xdv");
-    pars = getParValues(pdfParam);  
-    pars[0] = 1.;
-    sum0 = pdfParam->moment(pars.get(),-1);
-    _dSum = 1.0 / sum0;
-
+    getPdfParam("xuv")->setMoment(-1,2.0);
+    getPdfParam("xdv")->setMoment(-1,1.0);
     // momentum sum-rule
-
     // quark part
-    double xsumq = 0;
-    const std::vector<std::string> nameS{"xubar","xdbar","xs","xuv","xdv"};
-    for ( auto const& name : nameS ) {
-      const BasePdfParam* pdfParam = getPdfParam(name);
-      std::unique_ptr<double[]> pars = getParValues(pdfParam);
-
-      if ( name == "xuv") {
-	pars[0] = _uSum;
-      }
-      if (name == "xdv" ) {
-	pars[0] = _dSum;
-      }
-      xsumq += pdfParam->moment(pars.get(),0);
-    }
-    
+    double xsumq=0;
+    xsumq+=  getPdfParam("xuv"  )->moment(0);
+    xsumq+=  getPdfParam("xdv"  )->moment(0);
+    xsumq+=2*getPdfParam("xubar")->moment(0);
+    xsumq+=2*getPdfParam("xdbar")->moment(0);
+    xsumq+=2*getPdfParam("xs"   )->moment(0);
     // gluon part
-    pdfParam = getPdfParam("xg");
-    pars = getParValues(pdfParam);
-    pars[0] = 1.;
-    double xsumg = pdfParam->moment(pars.get(),0);
-    
-    _gSum = (1. - xsumq)/xsumg;
+    getPdfParam("xg")->setMoment(0,1-xsumq);
 
     printParams();
   }
@@ -131,19 +70,10 @@ namespace xfitter
       auto pars = getParValues(pdfParam);
       int  npar = pdfParam->getNPar();
 
-      if (name == "xg") {
-	pars[0] = _gSum;
-      }
-      if (name == "xuv") {
-	pars[0] = _uSum;
-      }
-      if (name == "xdv") {
-	pars[0] = _dSum;
-      }
       //      std::cout << "name :" << name << " npar: " << npar << "\n";
       std::cout << std::left<< std::setw(8) << name << std::right;
       for ( int i =0 ; i<npar; i++) {
-	std::cout << std::setw(12) << pars[i];
+        std::cout << std::setw(12) << pars[i];
       }
       std::cout << "\n";
     }
@@ -152,39 +82,27 @@ namespace xfitter
   
   std::unique_ptr<double[]> UvDvUbarDbarS::getParValues(BasePdfParam const* param) const {
     std::unique_ptr<double[]> pars( new double[param->getNPar()] );
-    for (int i=0; i<param->getNPar(); i++) {
-      const std::string pnam = getName() + "_" + param->getName() + "_p" +  std::to_string(i) ;
-      pars[i] = *XFITTER_PARS::gParameters[pnam];
+    for (unsigned int i=0; i<param->getNPar(); i++) {
+      const std::string pnam=param->getName()+ "_p" +  std::to_string(i) ;
+      pars[i] = *XFITTER_PARS::gParameters.at(pnam);
     }
     return pars;
   }
   
-  
-  double UvDvUbarDbarS::valence(double x, std::string const& name, double sum) const {
-    const BasePdfParam* pdfParam = getPdfParam(name);
-    // get parameters
-    std::unique_ptr<double[]> pars = getParValues(pdfParam);  
-
-    pars[0] = sum;
-    
-    return pdfParam->compute(x,pars.get());
+  double UvDvUbarDbarS::valence(double x, const std::string&name) const {
+    return(*getPdfParam(name))(x);
   }
 
-  double UvDvUbarDbarS::sea(double x, std::string const& name) const {
-    const BasePdfParam* pdfParam = getPdfParam(name);
-    // get parameters
-    std::unique_ptr<double[]> pars = getParValues(pdfParam);  
-    
-    return pdfParam->compute(x,pars.get());
+  double UvDvUbarDbarS::sea(double x,const std::string&name)const{
+    return(*getPdfParam(name))(x);
   }
 
-  
   double UvDvUbarDbarS::uv(double x) const {
-    return valence(x,"xuv",_uSum); 
+    return valence(x,"xuv"); 
   }
 
   double UvDvUbarDbarS::dv(double x) const {
-    return valence(x,"xdv",_dSum); 
+    return valence(x,"xdv"); 
   }
 
   double UvDvUbarDbarS::ubar(double x) const {
@@ -199,12 +117,8 @@ namespace xfitter
     return sea(x,"xs");
   }
 
-  double UvDvUbarDbarS::g(double x) const {
-    const BasePdfParam* pdfParam = getPdfParam("xg");
-    // get parameters
-    std::unique_ptr<double[]> pars = getParValues(pdfParam);  
-    pars[0] = _gSum;
-    return pdfParam->compute(x,pars.get());
+  double UvDvUbarDbarS::g(double x)const{
+    return (*getPdfParam("xg"))(x);
   }
 
   
diff --git a/pdfparams/BasePdfParam/include/BasePdfParam.h b/pdfparams/BasePdfParam/include/BasePdfParam.h
index 503f81c7f2ca556f4702cadaa84f1c1e18dcb915..fa279bfdcf1cda4e8808c9ed33a283eda6ddb77c 100644
--- a/pdfparams/BasePdfParam/include/BasePdfParam.h
+++ b/pdfparams/BasePdfParam/include/BasePdfParam.h
@@ -7,43 +7,48 @@
 /**
   @class BasePdfParam
 
-  @brief Base class to develop PDF parameterisations
+  @brief Base class for PDF parameterisations
 
-  Contains methods to compute PDF value at given x as well as integrals
+  Represents a function xf(x) that depends on x and some external parameters --- xf(x,A,B,C,...)
+  Has methods to
+  *Evaluate the function at given x and current (globally stored) parameters
+  *Calculate n-th moment, e.g. \int_0^1 x^n*xf(x) dx
+  *Rescale some parameters so that the n-th moment has a given value
 
-  @version 0.1
-  @date 2018-07-11
-  */
-
-class BasePdfParam {
- public:
-
-  /// Default constructor. Name is the PDF name
-  BasePdfParam(const std::string& inName): _name(inName),_npar(0) {}
+  Parameterisations keeps pointers to its parameters.
 
-  /// Set number of parameters
-  void SetNPar(int npar) { _npar = npar;}
-  
-  /// Compute xf(x,pars). Pure virtual method
-  virtual double compute( double const x, double const* pars) const = 0;
+  Parameterisations are initialised by parsing instance-specific YAML nodes
+  Concrete implementations of BasePdfParam must be able to get their parameters and any/all additional options
 
-  /// Compute moments of xf(x) ( for i=-1  of f(x) ), needed for sum-rules
-  virtual double moment( double const* pars, int const iMoment = 0) const ;
-
-  /// Return number of parameters
-  const int getNPar() const {return _npar;} 
-
-  /// Get name of the PDF object
-  const std::string getName() const {return _name;} 
+  @version 0.2
+  @date 2018-08-13
+  */
 
-  /// Get initial values from a yaml node. Uses node[getName] as the basis
-  double* initFromYaml(YAML::Node value) ;
-  
- private:
-  /// Name of the PDF object (e.g. uv, dv ...)
+class BasePdfParam{
+public:
+  BasePdfParam(const std::string&instance_name):_name(instance_name),pars{nullptr},Npars(0){}
+  ~BasePdfParam(){if(pars)delete[]pars;}
+  void               setNPar(unsigned int N){Npars=N;}
+  const unsigned int getNPar()const{return Npars;}
+  //!Evaluate xf(x) at given x with current parameters, pure virtual
+  virtual double operator()(double x)const=0;
+  //!Calculate n-th moment, e.g. \int_0^1 x^n*xf(x) dx
+  //!Note that we parameterise xf(x), not f(x)
+  //!Therefore, moment(-1) should be used for valence sums, and moment(0) for momentum sum
+  virtual double  moment(int nMoment=-1)const;
+  //!Rescale some parameters so that the n-th moment has a given value
+  //!A typical implementation will probably do this by setting *pars[0]
+  virtual void setMoment(int nMoment,double value);
+  //!Get name of the instance
+  const std::string getName()const{return _name;} 
+  //!Initialize from a yaml node. Uses node[getName] as the basis
+  virtual void initFromYaml(YAML::Node value);
+protected:
+  //!Unique name of instance
   std::string _name;
-  /// Number of parameters
-  int _npar;
-  
-  /// Vector of parameters. Parameters are pointers to doubles
+  //!Array of pointers to some global locations where minimization parameters are stored
+  //TODO: Implement proper pars initialization from YAML
+  double**pars;
+  //!Number of parameters, which is also the size of the array **parameters defined above
+  unsigned int Npars;
 };
diff --git a/pdfparams/BasePdfParam/src/BasePdfParam.cc b/pdfparams/BasePdfParam/src/BasePdfParam.cc
index 6f9e018225f1611d72f979833573c3d17bd50af4..1f321516d051c7563c94443508e5fd25d538b1da 100644
--- a/pdfparams/BasePdfParam/src/BasePdfParam.cc
+++ b/pdfparams/BasePdfParam/src/BasePdfParam.cc
@@ -1,10 +1,22 @@
-#include "BasePdfParam.h"
-#include <cmath>
-#include <memory>
-#include <iostream>
+#include"BasePdfParam.h"
+#include"xfitter_pars.h"
+#include<cmath>
+#include<memory>
+#include<iostream>
 
+/// TEMPORARY XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
+extern "C" {
+  /// Interface to minuit parameters
+  void addexternalparam_(const char name[],  const double &val, 
+                         const double  &step,
+                         const double &min, const double &max, 
+                         const double &prior, const double &priorUnc,
+                         const int &add, 
+                         map<std::string,double*> *map,
+                         int len);
+}
 /// Implement numeric integration
-double BasePdfParam::moment( double const* pars, int const iMoment) const {
+double BasePdfParam::moment(int iMoment)const{
   /// Simple rule, split log/lin spacing at xsplit=0.1
 
   const double xsplit = 0.1;
@@ -23,38 +35,52 @@ double BasePdfParam::moment( double const* pars, int const iMoment) const {
   // log x part:
   for (int i=0; i<nlog; i++) {
     double dx = pow(10,xminlog+(i+1)*xsteplog) - pow(10,xminlog+i*xsteplog);
-    double val = compute(x+dx/2.,pars)*pow(x+dx/2.,iMoment);
+    double val = (*this)(x+dx/2.)*pow(x+dx/2.,iMoment);
     x += dx;
     sum += dx*val;
   }
   // lin x part:
   for (int i=0; i<nlin; i++) {
     double dx = xsteplin;
-    double val = compute(x+dx/2.,pars)*pow(x+dx/2.,iMoment);
+    double val = (*this)(x+dx/2.)*pow(x+dx/2.,iMoment);
     x += dx;
     sum += dx*val;
   }
   
   return sum;
 }
-
-
-double*  BasePdfParam::initFromYaml(YAML::Node value) {
-  std::cout << " here here " << value << std::endl;
-  if ( value.IsSequence() ) {
-    size_t len = value.size();
-
-    std::cout << len << std::endl;
-    SetNPar(len);
-    double *pars = new double[len];
-
-    for (size_t i=0; i<len; i++) {
-      pars[i] = value[i].as<double>();
+void BasePdfParam::setMoment(int nMoment,double value){
+  *pars[0]=1;
+  *pars[0]=value/moment(nMoment);
+}
+void BasePdfParam::initFromYaml(YAML::Node value) {
+  // HARDWIRE A BIT FOR NOW XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
+  //TODO rewrite this for a different, new YAML format
+  using namespace std;
+  using uint=unsigned int;
+  cout<<"DEBUG["<<_name<<"]: initFromYaml: value="<<value<<endl;
+  if(value.IsSequence()){
+    Npars=value.size();
+    cout<<Npars<<endl;
+    pars=new double*[Npars];
+    // HARDWIRE old-way for now:  XXXXXXXXXXXXXXXXXXXXXXXXXXXXX
+    for(uint i=0;i<Npars;++i){
+      // get a good name
+      const std::string pnam=_name+"_p"+std::to_string(i);
+      double val     =value[i].as<double>();
+      double step    =fabs(val)/100.;       /// if 0, parameter is fixed !!! 
+      double minv    =0;
+      double maxv    =0;
+      double priorVal=0;
+      double priorUnc=0;
+      int add = true;
+      addexternalparam_(pnam.c_str(),val,step,minv,maxv,priorVal,priorUnc,add,&XFITTER_PARS::gParameters,pnam.size());
+			pars[i]=XFITTER_PARS::gParameters.at(pnam);
+      cout<<pnam<<"="<<(*pars[i])<<endl;
     }
-    return pars;
-  }
-  else {
-    return nullptr;
+  }else{
+    cout<<"ERROR["<<_name<<"]: initFromYaml: parameter is not a sequence!"<<endl;
   }
+  //TODO: Handle possible errors
 }
 
diff --git a/pdfparams/BasePdfParam/src/Makefile.am b/pdfparams/BasePdfParam/src/Makefile.am
index be9d7b4b42f75725f0a4572015685af8c7420638..7cf75604084672d8cc0c5e9944c579d5ffccb499 100644
--- a/pdfparams/BasePdfParam/src/Makefile.am
+++ b/pdfparams/BasePdfParam/src/Makefile.am
@@ -1,4 +1,4 @@
-AM_CXXFLAGS = -I$(srcdir)/../include  -Wall -fPIC -Wno-deprecated 
+AM_CXXFLAGS =-O0 -g -I$(srcdir)/../include -I$(srcdir)/../../../include -Wall -fPIC -Wno-deprecated 
 
 lib_LTLIBRARIES = libBasePdfParam_xfitter.la
 
diff --git a/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h b/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h
index 883529a6e24ee2f87374cecfe9048957cdeefd5b..f16ce7f9e0d9a093f5743ff675f60722e52ed5e9 100644
--- a/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h
+++ b/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h
@@ -8,20 +8,21 @@
 
   @brief A class for HERAPDF_ pdf parameterisation
 
-  Based on the ReactionTheory class. Reads options produces 3d cross section.
+  HERAPDF - style parameterisation:
+  xf(x)=A*x^B*(1-x)^C*(1+P(x))
+  where P(x) is a polynomial with other parameters as coefficients
 
-  @version 0.1
-  @date 2018-07-11
+  Number of parameters may vary, but at least 3, corresponding to A,B,C
+
+  In terms of indexed parameters:
+  xf(x)=par_0*x^par_1*(1-x)^par_2*(1+sum_{i=3}^{N}{par_i*x^(i-2)})
+
+  @version 0.2
+  @date 2018-08-14
   */
 
-class HERAPDF_PdfParam : public BasePdfParam
-{
+class HERAPDF_PdfParam:public BasePdfParam{
   public:
-     /// Default constructor. Name is the PDF name
-     HERAPDF_PdfParam (const std::string& inName) : BasePdfParam(inName) {}
-     /// Compute xf(x,pars). Pure virtual method
-     virtual double compute( double const x, double const* pars) const override final;
-    
-     /// (optionally) compute moments:
-     // virtual double moment( double const* pars, int const iMoment = 1) override final;
+    HERAPDF_PdfParam(const std::string&inName):BasePdfParam(inName){}
+    virtual double operator()(double x)const override final;
 };
diff --git a/pdfparams/HERAPDF_PdfParam/src/HERAPDF_PdfParam.cc b/pdfparams/HERAPDF_PdfParam/src/HERAPDF_PdfParam.cc
index ec9eaae7f8c736ea96a426f02b21ec48aea09a66..dd9336dd1d4fd7cddba0ecc76f071373c48a11d7 100644
--- a/pdfparams/HERAPDF_PdfParam/src/HERAPDF_PdfParam.cc
+++ b/pdfparams/HERAPDF_PdfParam/src/HERAPDF_PdfParam.cc
@@ -9,18 +9,17 @@
 #include <cmath>
 
 // Main function to compute PDF
-double HERAPDF_PdfParam::compute(double const x, double const* pars) const
-{
+double HERAPDF_PdfParam::operator()(double x)const{
   const int npar = getNPar();
   if (npar<3) {
     return NAN;
   }
-  double power = pars[0]*pow(x,pars[1])*pow((1-x),pars[2]);
+  double power=(*pars[0])*pow(x,(*pars[1]))*pow((1-x),(*pars[2]));
   double poly = 1;
   double xx = 1;
   for (int i = 3; i<npar; i++) {
     xx *= x;
-    poly += pars[i]*xx;
+    poly+=(*pars[i])*xx;
   }
   return power*poly;
 }
diff --git a/tools/AddPdfParam.py b/tools/AddPdfParam.py
index 01bc7e829ce249ae0f527dda09904a7c29af3fe3..227e9684fcf7e544ee93289f39fc17933f79b55a 100644
--- a/tools/AddPdfParam.py
+++ b/tools/AddPdfParam.py
@@ -35,26 +35,28 @@ with open(hFile,"w+") as f:
 #include "BasePdfParam.h"
 
 /**
-  @class {:s}PdfParam 
+  @class {name:s}PdfParam 
 
-  @brief A class for {:s} pdf parameterisation
+  @brief A class for {name:s} pdf parameterisation
 
   @version 0.1
-  @date {:s}
+  @date {date:s}
   */
 
-class {:s}PdfParam : public BasePdfParam
-{{
+class {name:s}PdfParam:public BasePdfParam{{
   public:
-     /// Default constructor. Name is the PDF name
-    {:s}PdfParam (const std::string& inName) : BasePdfParam(inName) {{}}
-     /// Compute xf(x,pars). Pure virtual method
-     virtual double compute( double const x, double const* pars) override final;
-    
-     /// (optionally) compute moments:
-     // virtual double moment( double const* pars, int const iMoment = 1) override final;
+     {name:s}PdfParam(const std::string&inName):BasePdfParam(inName){{}}
+     //Evaluate xf(x) at given x with current parameters
+     virtual double operator()(double x)const override final;
+     // (Optional) compute moments:
+     // virtual double  moment(int nMoment=-1)const override final;
+     // (Optional) set moments:
+     // virtual void setMoment(int nMoment,double value)override final;
+     // (Optional)
+     //Initialize from a yaml node. Uses node[getName] as the basis
+     // virtual void initFromYaml(YAML::Node value)override final;
 }};
-'''.format( name, name, datetime.date.today().isoformat(),name,name)
+'''.format(name=name,date=datetime.date.today().isoformat())
 )
 
 
@@ -65,29 +67,18 @@ print "Creating source file "+sFile
 with open(sFile,"w+") as f:
     f.write(''' 
 /*
-   @file {:s}PdfParam.cc
-   @date {:s}
-   @author  AddPdfParam.py
-   Created by  AddPdfParam.py on {:s}
+   @file {name:s}PdfParam.cc
+   @date {date:s}
+   @author AddPdfParam.py
+   Created by AddPdfParam.py on {date:s}
 */
 
-#include "{:s}PdfParam.h"
+#include "{name:s}PdfParam.h"
 
-
-// Main function to compute PDF
-double {:s}PdfParam::compute( double const x, double const* pars)
-{{
-  return 0;
+double {name:s}PdfParam::operator()(double x){{
+  //Your code here
 }}
-
-// Optional function to compute integrals:
-
-// double {:s}PdfParam::moment( double const* pars, int const iMoment)
-// {{
-//   return 0
-// }}
-
-'''.format(name,datetime.date.today().isoformat(),datetime.date.today().isoformat(),name,name,name)
+'''.format(name=name,date=datetime.date.today().isoformat())
 )
 
 aFile = "pdfparams/{:s}PdfParam/src/Makefile.am".format(name)