From 0862932ccdff43baa0610787cbb40c0913bcc995 Mon Sep 17 00:00:00 2001 From: Ivan Novikov <novivanya@yandex.ru> Date: Mon, 13 Aug 2018 14:36:50 +0200 Subject: [PATCH 1/5] Add GRV_PionPdfDecomposition --- Reactions.txt | 1 + .../include/GRV_PionPdfDecomposition.h | 50 ++++++++ .../src/GRV_PionPdfDecomposition.cc | 119 ++++++++++++++++++ .../GRV_PionPdfDecomposition/src/Makefile.am | 12 ++ .../yaml/parameters.yaml | 0 5 files changed, 182 insertions(+) create mode 100644 pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h create mode 100644 pdfdecompositions/GRV_PionPdfDecomposition/src/GRV_PionPdfDecomposition.cc create mode 100644 pdfdecompositions/GRV_PionPdfDecomposition/src/Makefile.am create mode 100644 pdfdecompositions/GRV_PionPdfDecomposition/yaml/parameters.yaml diff --git a/Reactions.txt b/Reactions.txt index 2fd18edaa..a3de58155 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 000000000..40ba22a99 --- /dev/null +++ b/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h @@ -0,0 +1,50 @@ + +#pragma once + +#include "BasePdfDecomposition.h" + +/** + @class GRV_PionPdfDecomposition + + @brief A class for GRV_Pion pdf decomposition + + @version 0.1 + @date 2018-08-07 + */ + +/* +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 +*/ + +namespace xfitter { +class ParameterisationWrapper; +class GRV_PionPdfDecomposition:public BasePdfDecomposition{ + public: + GRV_PionPdfDecomposition(); + /// Default constructor. Name is the PDF name + GRV_PionPdfDecomposition(const std::string& inName); + ~GRV_PionPdfDecomposition(); + /// Optional initialization at the first call + virtual void initAtStart(const std::string & pars) override final; + /// Compute PDF in a physical base in LHAPDF format for given x and Q + virtual std::function<std::map<int,double>(const double& x)>f0()const override final; + private: + ParameterisationWrapper*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 000000000..e337eb6d2 --- /dev/null +++ b/pdfdecompositions/GRV_PionPdfDecomposition/src/GRV_PionPdfDecomposition.cc @@ -0,0 +1,119 @@ +/* + @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; +namespace xfitter{ +class ParameterisationWrapper{ +public: + ParameterisationWrapper(BasePdfParam*base,const string decompositionName):base{base}{ + const uint N=base->getNPar(); + pars=new double[N]; + const string prefix=decompositionName+"_"+base->getName()+"_p"; + for(uint i=0;i<N;++i){ + pars[i]=XFITTER_PARS::gParameters[prefix+to_string(i)]; + } + } + ~ParameterisationWrapper(){delete[]pars;} + double operator(double x){ + std::unique_ptr<double[]>p(new double[N]); + for(uint i=0;i<N;++i)p[i]=*pars[i]; + return base->compute(x,p); + } + double moment(int i){//return \int_0^1 x^i*f(x) dx + const uint N=base->getNPar(); + std::unique_ptr<double[]>p(new double[N]); + for(uint i=0;i<N;++i)p[i]=*pars[i]; + return base->moment(p,i); + } + void setMoment(int i,double S){ + *pars[0]=1; + *pars[0]=S/moment(i); + } +private: + BasePdfParam*base; + double**pars; +}; +//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); + } + par_v =new ParameterisationWrapper(getPdfParam("v "),getName()); + par_qbar=new ParameterisationWrapper(getPdfParam("qbar"),getName()); + par_g =new ParameterisationWrapper(getPdfParam("g "),getName()); +} +void GRV_PionPdfDecomposition:initAtIteration() { + //Enforce sum rules + //Valence sum + par_v->setMoment(0,2); + //Momentum sum + par_g->setMoment(1,1-4*par_qbar->moment(1)); +} +// 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 000000000..49ee9e52a --- /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 000000000..e69de29bb -- GitLab From 1553e13b64a273fba0d6413a7146d88ef365da27 Mon Sep 17 00:00:00 2001 From: Ivan Novikov <novivanya@yandex.ru> Date: Tue, 14 Aug 2018 16:30:11 +0200 Subject: [PATCH 2/5] Rewrite BasePdfParam to keep pointers to parameters, and add setMoment to BasePdfParam --- .../include/GRV_PionPdfDecomposition.h | 15 +- .../src/GRV_PionPdfDecomposition.cc | 49 ++----- .../UvDvUbarDbarS/include/UvDvUbarDbarS.h | 7 +- .../UvDvUbarDbarS/src/UvDvUbarDbarS.cc | 138 ++++-------------- pdfparams/BasePdfParam/include/BasePdfParam.h | 72 ++++----- pdfparams/BasePdfParam/src/BasePdfParam.cc | 74 +++++++--- pdfparams/BasePdfParam/src/Makefile.am | 2 +- .../include/HERAPDF_PdfParam.h | 16 +- .../HERAPDF_PdfParam/src/HERAPDF_PdfParam.cc | 7 +- tools/AddPdfParam.py | 55 +++---- 10 files changed, 166 insertions(+), 269 deletions(-) diff --git a/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h b/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h index 40ba22a99..a0f0b0a03 100644 --- a/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h +++ b/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h @@ -1,15 +1,12 @@ - #pragma once - #include "BasePdfDecomposition.h" - /** @class GRV_PionPdfDecomposition @brief A class for GRV_Pion pdf decomposition - @version 0.1 - @date 2018-08-07 + @version 0.2 + @date 2018-08-14 */ /* @@ -32,19 +29,15 @@ And sum rules for pi- are: \int_0^1 x*(4*qbar+g) dx=1 */ -namespace xfitter { -class ParameterisationWrapper; +namespace xfitter{ class GRV_PionPdfDecomposition:public BasePdfDecomposition{ public: GRV_PionPdfDecomposition(); - /// Default constructor. Name is the PDF name GRV_PionPdfDecomposition(const std::string& inName); ~GRV_PionPdfDecomposition(); - /// Optional initialization at the first call virtual void initAtStart(const std::string & pars) override final; - /// Compute PDF in a physical base in LHAPDF format for given x and Q virtual std::function<std::map<int,double>(const double& x)>f0()const override final; private: - ParameterisationWrapper*par_v,*par_qbar,*par_g; + 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 index e337eb6d2..8cad3faa5 100644 --- a/pdfdecompositions/GRV_PionPdfDecomposition/src/GRV_PionPdfDecomposition.cc +++ b/pdfdecompositions/GRV_PionPdfDecomposition/src/GRV_PionPdfDecomposition.cc @@ -7,37 +7,15 @@ #include"GRV_PionPdfDecomposition.h" using uint=unsigned int; using namespace std; -namespace xfitter{ -class ParameterisationWrapper{ -public: - ParameterisationWrapper(BasePdfParam*base,const string decompositionName):base{base}{ - const uint N=base->getNPar(); - pars=new double[N]; - const string prefix=decompositionName+"_"+base->getName()+"_p"; - for(uint i=0;i<N;++i){ - pars[i]=XFITTER_PARS::gParameters[prefix+to_string(i)]; - } - } - ~ParameterisationWrapper(){delete[]pars;} - double operator(double x){ - std::unique_ptr<double[]>p(new double[N]); - for(uint i=0;i<N;++i)p[i]=*pars[i]; - return base->compute(x,p); - } - double moment(int i){//return \int_0^1 x^i*f(x) dx - const uint N=base->getNPar(); - std::unique_ptr<double[]>p(new double[N]); - for(uint i=0;i<N;++i)p[i]=*pars[i]; - return base->moment(p,i); +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)]; } - void setMoment(int i,double S){ - *pars[0]=1; - *pars[0]=S/moment(i); - } -private: - BasePdfParam*base; - double**pars; -}; +} +namespace xfitter{ //For dynamic loading: extern "C" GRV_PionPdfDecomposition*create(){ return new GRV_PionPdfDecomposition(); @@ -79,16 +57,17 @@ void GRV_PionPdfDecomposition::initAtStart(const std::string & pars){ } addParameterisation(pdfName,pParam); } - par_v =new ParameterisationWrapper(getPdfParam("v "),getName()); - par_qbar=new ParameterisationWrapper(getPdfParam("qbar"),getName()); - par_g =new ParameterisationWrapper(getPdfParam("g "),getName()); +//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(0,2); + par_v->setMoment(-1,2); //Momentum sum - par_g->setMoment(1,1-4*par_qbar->moment(1)); + 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{ diff --git a/pdfdecompositions/UvDvUbarDbarS/include/UvDvUbarDbarS.h b/pdfdecompositions/UvDvUbarDbarS/include/UvDvUbarDbarS.h index 89e55a10b..74e06bcc2 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 d799d8591..73184464b 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,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++) { + for (unsigned 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]; } 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 503f81c7f..d258030e6 100644 --- a/pdfparams/BasePdfParam/include/BasePdfParam.h +++ b/pdfparams/BasePdfParam/include/BasePdfParam.h @@ -7,43 +7,49 @@ /** @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: + //Default constructor. Name is the unique name of instance + 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 parametrise 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 6f9e01822..bdaf74584 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; + (*pars[i])=val; + addexternalparam_(pnam.c_str(),val,step,minv,maxv,priorVal,priorUnc,add,&XFITTER_PARS::gParameters,pnam.size()); + cout<<pnam<<"="<<val<<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 be9d7b4b4..0a0c1a13f 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 = -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 883529a6e..66b37959c 100644 --- a/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h +++ b/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h @@ -10,18 +10,12 @@ Based on the ReactionTheory class. Reads options produces 3d cross section. - @version 0.1 - @date 2018-07-11 + @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 ec9eaae7f..dd9336dd1 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 01bc7e829..227e9684f 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) -- GitLab From 730a839cde493b6299d0c9f9ffb96efad6db6e58 Mon Sep 17 00:00:00 2001 From: Ivan Novikov <novivanya@yandex.ru> Date: Wed, 15 Aug 2018 13:42:27 +0200 Subject: [PATCH 3/5] Rewrite comments for doxygen --- .../include/GRV_PionPdfDecomposition.h | 38 +++++++++---------- pdfparams/BasePdfParam/include/BasePdfParam.h | 23 ++++++----- .../include/HERAPDF_PdfParam.h | 9 ++++- 3 files changed, 36 insertions(+), 34 deletions(-) diff --git a/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h b/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h index a0f0b0a03..465644ff1 100644 --- a/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h +++ b/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h @@ -5,30 +5,26 @@ @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 */ - -/* -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 -*/ - namespace xfitter{ class GRV_PionPdfDecomposition:public BasePdfDecomposition{ public: diff --git a/pdfparams/BasePdfParam/include/BasePdfParam.h b/pdfparams/BasePdfParam/include/BasePdfParam.h index d258030e6..fa279bfdc 100644 --- a/pdfparams/BasePdfParam/include/BasePdfParam.h +++ b/pdfparams/BasePdfParam/include/BasePdfParam.h @@ -26,30 +26,29 @@ class BasePdfParam{ public: - //Default constructor. Name is the unique name of instance 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 + //!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 parametrise xf(x), not f(x) - //Therefore, moment(-1) should be used for valence sums, and moment(0) for momentum sum + //!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] + //!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 + //!Get name of the instance const std::string getName()const{return _name;} - //Initialize from a yaml node. Uses node[getName] as the basis + //!Initialize from a yaml node. Uses node[getName] as the basis virtual void initFromYaml(YAML::Node value); protected: - //Unique name of instance + //!Unique name of instance std::string _name; - //Array of pointers to some global locations where minimization parameters are stored + //!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 + //!Number of parameters, which is also the size of the array **parameters defined above unsigned int Npars; }; diff --git a/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h b/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h index 66b37959c..f16ce7f9e 100644 --- a/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h +++ b/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h @@ -8,7 +8,14 @@ @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 + + 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 -- GitLab From 6c08a888280fb4b3a68f203c942ecd982b8280da Mon Sep 17 00:00:00 2001 From: Ivan Novikov <novivanya@yandex.ru> Date: Wed, 15 Aug 2018 15:14:03 +0200 Subject: [PATCH 4/5] Fix segfaults --- pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc | 4 ++-- pdfparams/BasePdfParam/src/BasePdfParam.cc | 4 ++-- pdfparams/BasePdfParam/src/Makefile.am | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc b/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc index 73184464b..16aab69fc 100644 --- a/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc +++ b/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc @@ -83,8 +83,8 @@ namespace xfitter std::unique_ptr<double[]> UvDvUbarDbarS::getParValues(BasePdfParam const* param) const { std::unique_ptr<double[]> pars( new double[param->getNPar()] ); for (unsigned 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]; + const std::string pnam=param->getName()+ "_p" + std::to_string(i) ; + pars[i] = *XFITTER_PARS::gParameters.at(pnam); } return pars; } diff --git a/pdfparams/BasePdfParam/src/BasePdfParam.cc b/pdfparams/BasePdfParam/src/BasePdfParam.cc index bdaf74584..1f321516d 100644 --- a/pdfparams/BasePdfParam/src/BasePdfParam.cc +++ b/pdfparams/BasePdfParam/src/BasePdfParam.cc @@ -74,9 +74,9 @@ void BasePdfParam::initFromYaml(YAML::Node value) { double priorVal=0; double priorUnc=0; int add = true; - (*pars[i])=val; addexternalparam_(pnam.c_str(),val,step,minv,maxv,priorVal,priorUnc,add,&XFITTER_PARS::gParameters,pnam.size()); - cout<<pnam<<"="<<val<<endl; + pars[i]=XFITTER_PARS::gParameters.at(pnam); + cout<<pnam<<"="<<(*pars[i])<<endl; } }else{ cout<<"ERROR["<<_name<<"]: initFromYaml: parameter is not a sequence!"<<endl; diff --git a/pdfparams/BasePdfParam/src/Makefile.am b/pdfparams/BasePdfParam/src/Makefile.am index 0a0c1a13f..7cf756040 100644 --- a/pdfparams/BasePdfParam/src/Makefile.am +++ b/pdfparams/BasePdfParam/src/Makefile.am @@ -1,4 +1,4 @@ -AM_CXXFLAGS = -I$(srcdir)/../include -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 -- GitLab From f1e0fc4abbc651a5a97a880712b0ae902ab080bd Mon Sep 17 00:00:00 2001 From: Ivan Novikov <novivanya@yandex.ru> Date: Thu, 16 Aug 2018 10:20:17 +0200 Subject: [PATCH 5/5] Fix momentum sum rule in UvDvUbarDbarS --- pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc b/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc index 16aab69fc..6ca48015e 100644 --- a/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc +++ b/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc @@ -53,7 +53,7 @@ namespace xfitter xsumq+=2*getPdfParam("xdbar")->moment(0); xsumq+=2*getPdfParam("xs" )->moment(0); // gluon part - getPdfParam("xg")->setMoment(0,xsumq); + getPdfParam("xg")->setMoment(0,1-xsumq); printParams(); } -- GitLab