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)