diff --git a/Makefile.am b/Makefile.am index 81480791c5b88051201979dccf2d572b1c3fa5bd..480f0e4628577d4fc673ce4385d6f6123d775242 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,32 +1,35 @@ ACLOCAL_AMFLAGS=-I m4 AUTOMAKE_OPTIONS = foreign SUBDIRS = minuit/src interfaces/src DY/src DIPOLE/src RT/src EW/src common common/linalg \ - common/num_utils pdf2yaml tools/process evolutions/LHAPDF/src evolutions/QCDNUM/src \ - tools/draw tools/MakeLHAPDF FastNLO/src DiffDIS/src ACOT/src SACOT/src ABM/src FONLL/src Cascade/src \ - genetic/mixmax_r004 genetic/src QEDevol/src \ - include interfaces/include FastNLO/include FastNLO/include/fastnlotk DiffDIS/include \ - DY/include tools/draw/include \ - reactions/KFactor/src \ + common/num_utils pdf2yaml tools/process\ + evolutions/LHAPDF/src evolutions/QCDNUM/src\ + tools/draw tools/MakeLHAPDF FastNLO/src DiffDIS/src ACOT/src SACOT/src ABM/src FONLL/src Cascade/src \ + genetic/mixmax_r004 genetic/src QEDevol/src \ + include interfaces/include FastNLO/include FastNLO/include/fastnlotk DiffDIS/include \ + DY/include tools/draw/include \ + reactions/KFactor/src \ reactions/Fractal_DISNC/src \ reactions/BaseDISCC/src \ reactions/Hathor/src \ reactions/BaseDISNC/src \ - reactions/RT_DISNC/src \ + reactions/RT_DISNC/src \ reactions/FFABM_DISNC/src \ reactions/FFABM_DISCC/src \ reactions/APPLgrid/src \ reactions/BaseHVQMNR/src \ - reactions/HVQMNR_LHCb_7TeV_beauty/src \ + reactions/HVQMNR_LHCb_7TeV_beauty/src \ reactions/HVQMNR_LHCb_7TeV_charm/src \ reactions/testZMVFNS/src \ reactions/fastNLO/src/ \ - reactions/FONLL_DISCC/src \ + reactions/FONLL_DISCC/src \ reactions/FONLL_DISNC/src \ pdfparams/BasePdfParam/src \ - pdfparams/HERAPDF_PdfParam/src \ + pdfparams/HERAPDF_PdfParam/src \ + pdfparams/PolySqrtPdfParam/src \ pdfdecompositions/BasePdfDecomposition/src \ pdfdecompositions/LHAPDFDecomposition/src \ pdfdecompositions/UvDvUbarDbarS/src \ + pdfdecompositions/SU3_PionPdfDecomposition/src \ evolutions/BaseEvolution/src \ evolutions/APFELxx/src \ minimizers/BaseMinimizer/src minimizers/CERESMinimizer/src minimizers/MINUITMinimizer/src diff --git a/Reactions.txt b/Reactions.txt index f7f8292f80287907a86f30d31c6f1d5d7a875b28..a87355f2e9da9cf6b78bfbbfed9a2f086467954d 100644 --- a/Reactions.txt +++ b/Reactions.txt @@ -18,7 +18,9 @@ BaseEvolution libbaseevolution_xfitter.so APFELxx libapfelxx_xfitter.so QCDNUM libqcdnum_xfitter.so UvDvubardbars libuvdvubardbars_xfitter.so -GRV_PionPdfDecomposition libgrv_pionpdfdecomposition_xfitter.so +SU3_Pion libSU3_PionPdfDecomposition_xfitter.so MINUIT libMINUITMinimizer_xfitter.so CERES libCERESMinimizer_xfitter.so LHAPDF liblhapdf_xfitter.so +HERAPDF libHERAPDF_PdfParam_xfitter.so +PolySqrt libPolySqrtPdfParam_xfitter.so diff --git a/configure.ac b/configure.ac index f4f4b5313c7b59e4e28345a4cfb2fee072a3e4cf..00f7f37ef560137b24c1603736b4b7bcb6e8fde5 100644 --- a/configure.ac +++ b/configure.ac @@ -605,8 +605,10 @@ AC_CONFIG_FILES([include/Makefile evolutions/QCDNUM/src/Makefile pdfparams/BasePdfParam/src/Makefile pdfparams/HERAPDF_PdfParam/src/Makefile + pdfparams/PolySqrtPdfParam/src/Makefile pdfdecompositions/LHAPDFDecomposition/src/Makefile pdfdecompositions/UvDvUbarDbarS/src/Makefile + pdfdecompositions/SU3_PionPdfDecomposition/src/Makefile pdfdecompositions/BasePdfDecomposition/src/Makefile reactions/KFactor/src/Makefile reactions/BaseDISCC/src/Makefile diff --git a/doxygen.cfg b/doxygen.cfg index 3b8eca87b41cd1d8410e17c7fa039a5a6d5ba949..829df2fc246c5c0ba350e1493028a72eadcaf976 100644 --- a/doxygen.cfg +++ b/doxygen.cfg @@ -582,10 +582,16 @@ INPUT = include src \ reactions/BaseHVQMNR/src reactions/HVQMNR_LHCb_7TeV_beauty/include \ reactions/HVQMNR_LHCb_7TeV_beauty/src \ - pdfparams/BasePdfParam/include pdfparams/HERAPDF_PdfParam/include \ - pdfdecompositions/BasePdfDecomposition/include pdfdecompositions/LHAPDFDecomposition/include \ + pdfparams/BasePdfParam/include \ + pdfparams/HERAPDF_PdfParam/include \ + pdfparams/PolySqrtPdfParam/include \ + pdfdecompositions/BasePdfDecomposition/include \ + pdfdecompositions/LHAPDFDecomposition/include \ pdfdecompositions/UvDvUbarDbarS/include \ - minimizers/BaseMinimizer/include minimizers/CERESMinimizer/include minimizers/MINUITMinimizer/include \ + pdfdecompositions/GRV_PionPdfDecomposition/include \ + minimizers/BaseMinimizer/include \ + minimizers/CERESMinimizer/include \ + minimizers/MINUITMinimizer/include \ # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is # also the default input encoding. Doxygen uses libiconv (or the iconv built diff --git a/evolutions/APFELxx/include/EvolutionAPFELxx.h b/evolutions/APFELxx/include/EvolutionAPFELxx.h index f5e3114215f60344ba669e0b5e20524aa0e44452..563866dce415c6c62d54df234eff040a9719abea 100644 --- a/evolutions/APFELxx/include/EvolutionAPFELxx.h +++ b/evolutions/APFELxx/include/EvolutionAPFELxx.h @@ -18,28 +18,20 @@ namespace xfitter @brief Derived class of BaseEvolution for using APFEL++ as an evolution code. - @version 0.1 - @date 2018-07-11 + @version 0.2 + @date 2018-09-29 */ class EvolutionAPFELxx: public BaseEvolution { public: - /// Empty constructor (needed for the dynamic loading) - EvolutionAPFELxx(): BaseEvolution{"APFELxx",nullptr} {} - - /// Constructor wit PDF decomposition - EvolutionAPFELxx(std::function<std::map<int,double>(double const& x)> const& inPDFs): BaseEvolution{"APFELxx", inPDFs} {} + EvolutionAPFELxx(const char*name):BaseEvolution{name}{} + virtual const char*getClassName()const final override; /** * @brief Function that initialises the evolution in APFEL++. */ - void initAtStart(); - - /** - * @brief Function that updates the relevant parameters of the - * evolution at each iteration of the fitting procedure. - */ - void initAtIteration(); + virtual void atStart()override final; + virtual void atConfigurationChange()override final; /** * @name Getters @@ -76,6 +68,7 @@ namespace xfitter ///@} private: + std::function<std::map<int,double>(const double& x)> _inPDFs; std::vector<double> _Masses; std::vector<double> _Thresholds; std::unique_ptr<const apfel::Grid> _Grid; diff --git a/evolutions/APFELxx/src/EvolutionAPFELxx.cc b/evolutions/APFELxx/src/EvolutionAPFELxx.cc index 6acaf7d1fcfdfbef1af37d44b5adbf13af4149f9..bdc18aea17a2b55d993ab5a4987173eb79065bac 100644 --- a/evolutions/APFELxx/src/EvolutionAPFELxx.cc +++ b/evolutions/APFELxx/src/EvolutionAPFELxx.cc @@ -10,22 +10,25 @@ namespace xfitter { // the class factories - extern "C" EvolutionAPFELxx* create() { - return new EvolutionAPFELxx(); + extern "C" EvolutionAPFELxx*create(const char*name){ + return new EvolutionAPFELxx(name); } + const char*EvolutionAPFELxx::getClassName()const{return "APFELxx";} //_________________________________________________________________________________ - void EvolutionAPFELxx::initAtStart() + void EvolutionAPFELxx::atStart() { // APFEL++ banner apfel::Banner(); + const YAML::Node yamlNode=XFITTER_PARS::getEvolutionNode(_name); + _inPDFs=XFITTER_PARS::getInputFunctionFromYaml(yamlNode); // Retrieve parameters needed to initialize APFEL++. const double* MCharm = XFITTER_PARS::gParameters.at("mch"); const double* MBottom = XFITTER_PARS::gParameters.at("mbt"); const double* MTop = XFITTER_PARS::gParameters.at("mtp"); - const YAML::Node xGrid = XFITTER_PARS::gParametersY.at("APFELxx")["xGrid"]; + const YAML::Node xGrid = yamlNode["xGrid"]; vector<apfel::SubGrid> sgv; for(auto const& sg : xGrid) @@ -41,17 +44,19 @@ namespace xfitter // Initialize QCD evolution objects _DglapObj = apfel::InitializeDglapObjectsQCD(*_Grid, _Masses, _Thresholds); + atConfigurationChange(); } //_________________________________________________________________________________ - void EvolutionAPFELxx::initAtIteration() + void EvolutionAPFELxx::atConfigurationChange() { + const YAML::Node yamlNode=XFITTER_PARS::getEvolutionNode(_name); // Retrieve the relevant parameters needed to compute the evolutions const int PtOrder = OrderMap(XFITTER_PARS::gParametersS.at("Order")) - 1; const double* Q0 = XFITTER_PARS::gParameters.at("Q0"); const double* Q_ref = XFITTER_PARS::gParameters.at("Mz"); const double* Alphas_ref = XFITTER_PARS::gParameters.at("alphas"); - const YAML::Node QGrid = XFITTER_PARS::gParametersY.at("APFELxx")["QGrid"]; + const YAML::Node QGrid = yamlNode["QGrid"]; // Reinitialise and tabulate the running coupling at every // iteration. This is fast enough and allows for the reference @@ -62,17 +67,17 @@ namespace xfitter // Construct the DGLAP objects const auto Dglap = BuildDglap(_DglapObj, - [=] (double const& x, double const&)->std::map<int,double>{ return apfel::PhysToQCDEv(this->_inPDFs(x)); }, - *Q0, PtOrder, _AlphaQCD); + [=] (double const& x, double const&)->std::map<int,double>{ return apfel::PhysToQCDEv(this->_inPDFs(x)); }, + *Q0, PtOrder, _AlphaQCD); // Tabulate PDFs (ideally the parameters of the tabulation should // be read from parameters.yaml). _TabulatedPDFs = std::unique_ptr<apfel::TabulateObject<apfel::Set<apfel::Distribution>>> (new apfel::TabulateObject<apfel::Set<apfel::Distribution>>{*Dglap, - QGrid[0].as<int>(), - QGrid[1].as<double>(), - QGrid[2].as<double>(), - QGrid[3].as<int>()}); + QGrid[0].as<int>(), + QGrid[1].as<double>(), + QGrid[2].as<double>(), + QGrid[3].as<int>()}); } //_________________________________________________________________________________ diff --git a/evolutions/BaseEvolution/include/BaseEvolution.h b/evolutions/BaseEvolution/include/BaseEvolution.h index a5c7c5ebf39c2389d942bb1893be4aec00bd1487..a1feedfe30075ddd85bc4320129b12c38fd60b02 100644 --- a/evolutions/BaseEvolution/include/BaseEvolution.h +++ b/evolutions/BaseEvolution/include/BaseEvolution.h @@ -3,79 +3,70 @@ #include <string> #include <map> #include <functional> +#include <yaml-cpp/yaml.h> namespace xfitter { /** @class BaseEvolution - @brief Base class for the evolving quantities + @brief Base class for the evolved PDFs and other evolved quantities Contains methods to compute the evolution of PDFs, alpha_s, and other possible evolving quantites. - @version 0.1 - @date 2018-07-11 + @version 0.2 + @date 2018-09-29 */ - - - - class BaseEvolution { public: + /// Unique name of instance + const std::string _name; /** * @brief The BaseEvolution default constructor. - * @param name: the name assignet to the instance + * @param name: the unique name used to identify the instance */ - BaseEvolution(const std::string& name, std::function<std::map<int,double>(double const& x)> const& inPDFs): _name(name), _inPDFs(inPDFs) { } - - /// Explicitly set PDF decomposition - void SetPdfDecomposition( std::function<std::map<int,double>(double const& x)> const& inPDFs) { _inPDFs = inPDFs;} + BaseEvolution(const char*name):_name(name){} /** * @brief Function to be called at the begining to initialise the - * evolution code. + * evolution code based on its YAML node + * + * This function is called only once */ - virtual void initAtStart() = 0; + virtual void atStart(){}; /** - * @brief Function to be call at each iteration to update the - * relevant evolution parameters. + * @brief Function to be called at each iteration */ - virtual void initAtIteration() = 0; + virtual void atIteration(){}; /** - * @name Setters + * @brief This function should be called when at least one parameter in the YAML node of given evolution changes */ - ///@{ - /** - * @brief Function to set the PdfDecomposition object to be used - * as initial scale distributions. - */ - void SetInitialPDFs(std::function<std::map<int,double>(double const& x)> const& inPDFs) { _inPDFs = inPDFs; } - ///@} + virtual void atConfigurationChange(){}; /** * @name Getters */ ///@{ /** - * @brief Function that returns a std::function that in turns + * @brief Function that returns a std::function that in turn * returns a map<int, double> as a function of x and Q. * @return map<int, double>-valued function of x and Q. */ virtual std::function<std::map<int,double>(double const& x, double const& Q)> xfxQMap() = 0; /** - * @brief Function that returns a std::function that in turns + * @brief Function that returns a std::function that in turn * returns a double as a function of the pdf index i, x and Q. * @return double-valued function of i, x and Q. */ virtual std::function<double(int const& i, double const& x, double const& Q)> xfxQDouble() = 0; /** - * @brief Function that returns a std::function that in turns + * @brief Function that returns a std::function that in turn * returns a void as a function of the pdf index x, Q, and pdfs, * where pdfs is the array of PDFs. * @return void-valued function of x, Q and pdfs. @@ -89,9 +80,6 @@ namespace xfitter */ virtual std::function<double(double const& Q)> AlphaQCD() = 0; - /// Get name - const std::string getName() const { return _name; } - /// Get generic property of the evolution virtual std::string getPropertyS(std::string const& propertyName ) const { return "" ; } @@ -101,15 +89,16 @@ namespace xfitter /// Get generic property of the evolution virtual double getPropertyD(std::string const& propertyName ) const { return 0.; } + /// Get class name, can be used to verify that the correct concrete class is being used + virtual const char*getClassName()const=0; ///@} protected: - const std::string _name; std::function<std::map<int,double>(double const& x)> _inPDFs; }; /// For dynamic loader - typedef BaseEvolution* create_evolution(); + typedef BaseEvolution*create_evolution(const char*name); } diff --git a/evolutions/BaseEvolution/src/BaseEvolution.cc b/evolutions/BaseEvolution/src/BaseEvolution.cc index f3443f6f07841787d61f6a4439684f548aff86d5..0081c72b2b96d13e502e23e183517d98c294aa35 100644 --- a/evolutions/BaseEvolution/src/BaseEvolution.cc +++ b/evolutions/BaseEvolution/src/BaseEvolution.cc @@ -1 +1,2 @@ #include "BaseEvolution.h" +//Why is this here??? diff --git a/evolutions/LHAPDF/include/EvolutionLHAPDF.h b/evolutions/LHAPDF/include/EvolutionLHAPDF.h index 59bdaa59791ee25d64630a5e80e0c24295505814..89c0935872c6e8e05da9f1b2acfe28adedf0b190 100644 --- a/evolutions/LHAPDF/include/EvolutionLHAPDF.h +++ b/evolutions/LHAPDF/include/EvolutionLHAPDF.h @@ -11,7 +11,12 @@ */ #include "BaseEvolution.h" +//Try to suppress unused-local-typedef warning from boost 1.53.0 for gcc +//Apparently these warnings have been fixed in later versions of boost +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-local-typedefs" #include "LHAPDF/LHAPDF.h" +#pragma GCC diagnostic pop namespace xfitter { @@ -20,15 +25,14 @@ class EvolutionLHAPDF : BaseEvolution { public: /// Empty constructor (needed for the dynamic loading) - EvolutionLHAPDF(): BaseEvolution("LHAPDF",nullptr) {}; + EvolutionLHAPDF(const char*name); + virtual const char*getClassName()const override final; public: - /// Constructor setting the name - virtual std::string getEvolutionName() const { return "LHAPDF" ;}; /// Global initialization - virtual void initAtStart() override final; - /// Init at each iteration - virtual void initAtIteration() override final; + virtual void atStart() override final; + /// Init at each change of at least one parameter + virtual void atConfigurationChange() override final; /// Return PDFs as a map <int,double> where int is PDF ID (-6, ... 6, 21) virtual std::function<std::map<int,double>(double const& x, double const& Q)> xfxQMap() override final; diff --git a/evolutions/LHAPDF/src/EvolutionLHAPDF.cc b/evolutions/LHAPDF/src/EvolutionLHAPDF.cc index bf8c57a4d988414e04028921b9cceeb10a2325c2..e9c613c04a18becad424648342bfc763d944012a 100644 --- a/evolutions/LHAPDF/src/EvolutionLHAPDF.cc +++ b/evolutions/LHAPDF/src/EvolutionLHAPDF.cc @@ -15,34 +15,41 @@ namespace xfitter { // the class factories -extern "C" EvolutionLHAPDF* create() { - return new EvolutionLHAPDF(); +extern "C" EvolutionLHAPDF*create(const char*name){ + return new EvolutionLHAPDF(name); } +EvolutionLHAPDF::EvolutionLHAPDF(const char*name):BaseEvolution(name){ + _pdf=nullptr; +} +const char*EvolutionLHAPDF::getClassName()const{return "LHAPDF";} /// Global initialization - void EvolutionLHAPDF::initAtStart() { - - // Initialize LHAPDF XXXXXXXXXXXXXXXXXXXXXXxxxxxxxxxxxx - auto pars = XFITTER_PARS::gParametersY["LHAPDF"]; - - _set_name = pars["set"].as<std::string>(); + void EvolutionLHAPDF::atStart() { + YAML::Node pars=XFITTER_PARS::getEvolutionNode(_name); + try{ + _set_name = pars["set"].as<std::string>(); + }catch(YAML::TypedBadConversion<std::string>&ex){ + hf_errlog(18090310,"F: In EvolutionLHAPDF::atStart: failed to convert YAML node \"set\" to string; printing node to stderr"); + std::cerr<<pars<<std::endl; + } // check if exists first CheckForPDF(_set_name.c_str()); - return ; + _member = pars["member"].as<int>(); + if(_pdf)delete _pdf; + _pdf = LHAPDF::mkPDF(_set_name,_member); }; /// Init at each iteration - void EvolutionLHAPDF::initAtIteration() { - - // Initialize LHAPDF XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX - auto pars = XFITTER_PARS::gParametersY["LHAPDF"]; - + void EvolutionLHAPDF::atConfigurationChange() { + YAML::Node pars=XFITTER_PARS::getEvolutionNode(_name); + //TODO: check for errors while parsing YAML _set_name = pars["set"].as<std::string>(); _member = pars["member"].as<int>(); - _pdf = LHAPDF::mkPDF(_set_name,_member); + if(_pdf)delete _pdf; + _pdf=LHAPDF::mkPDF(_set_name,_member); return ; }; diff --git a/evolutions/QCDNUM/include/EvolutionQCDNUM.h b/evolutions/QCDNUM/include/EvolutionQCDNUM.h index 44d9cfefe9eacd9e2b829948883e60f64e7e2dcc..646d09af475ce59f004ed81b4c07609baf0f1d4f 100644 --- a/evolutions/QCDNUM/include/EvolutionQCDNUM.h +++ b/evolutions/QCDNUM/include/EvolutionQCDNUM.h @@ -13,8 +13,8 @@ @brief A wrapper class for QCDNUM evolution - @version 0.1 - @date 2018-08-14 + @version 0.2 + @date 2018-09-29 */ namespace xfitter @@ -23,19 +23,18 @@ namespace xfitter class EvolutionQCDNUM: public BaseEvolution { public: - /// Empty constructor (needed for the dynamic loading) - EvolutionQCDNUM(): BaseEvolution{"QCDNUM",nullptr}{}; - - public: - virtual std::string getEvolutionName() const { return "QCDNUM" ;}; - virtual void initAtStart() override final; - virtual void initAtIteration() override final; + EvolutionQCDNUM(const char*name):BaseEvolution{name}{}; + + virtual const char*getClassName()const override final; + virtual void atStart()override final; + virtual void atIteration()override final; + virtual void atConfigurationChange()override final; virtual std::function<std::map<int,double>(double const& x, double const& Q)> xfxQMap() override final; virtual std::function<void(double const& x, double const& Q, double* pdfs)> xfxQArray() override final; virtual std::function<double(int const& i, double const& x, double const& Q)> xfxQDouble() override final; virtual std::function<double(double const& Q)> AlphaQCD() override final; - protected: - virtual int parseOptions(){ return 0;}; +// protected: +// virtual int parseOptions(){ return 0;}; private: /// PDFs called outside boundaries check: int _icheck{0}; diff --git a/evolutions/QCDNUM/src/EvolutionQCDNUM.cc b/evolutions/QCDNUM/src/EvolutionQCDNUM.cc index 055ce7b49f4c0cf3248a73775bfe2674072870ca..3855a334dcc7a4a96b52f550b389b27fb9d95573 100644 --- a/evolutions/QCDNUM/src/EvolutionQCDNUM.cc +++ b/evolutions/QCDNUM/src/EvolutionQCDNUM.cc @@ -23,8 +23,11 @@ double funcPDF(int *ipdf, double *x) { // helper to parse yaml sequences of uniform type template <class T> -vector<T> getSeq(std::string name) { - const YAML::Node node = XFITTER_PARS::gParametersY.at("QCDNUM")[name]; +vector<T> getSeq(const YAML::Node node) { + if(!node.IsSequence()){ + std::cerr<<"[ERROR]getSeq: node=\n"<<node<<std::endl; + hf_errlog(180829150,"F: In QCDNUM in function getSeq: wrong node type, expected sequence"); + } size_t len = node.size(); vector<T> v(len); for (size_t i=0; i<len; i++) { @@ -55,14 +58,13 @@ namespace xfitter { // the class factories - extern "C" EvolutionQCDNUM* create() { - return new EvolutionQCDNUM(); + extern "C" EvolutionQCDNUM* create(const char*name) { + return new EvolutionQCDNUM(name); } - - + const char*EvolutionQCDNUM::getClassName()const{return "QCDNUM";} // Initialize at the start of the computation - void EvolutionQCDNUM::initAtStart() + void EvolutionQCDNUM::atStart() { QCDNUM::qcinit(6," "); @@ -81,9 +83,9 @@ namespace xfitter QCDNUM::getint(memT[i].c_str(),itest); if (itest < siz[i]) { - std::string mess = "F: QCDNUM memory allocation insufficient. Recompile QCDNUM with at least "+memT[i] + " = " + std::to_string(siz[i]) + " in qcdnum.inc."; - hf_errlog(231020152,mess.c_str()); - } + std::string mess = "F: QCDNUM memory allocation insufficient. Recompile QCDNUM with at least "+memT[i] + " = " + std::to_string(siz[i]) + " in qcdnum.inc."; + hf_errlog(231020152,mess.c_str()); + } } // Evolution order: @@ -100,17 +102,16 @@ namespace xfitter // const double* mtp = XFITTER_PARS::gParameters.at("mtp"); // no top PDF treatment yet XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX - const YAML::Node yQCDNUM = XFITTER_PARS::gParametersY.at("QCDNUM"); - + YAML::Node yQCDNUM=XFITTER_PARS::getEvolutionNode(_name); _icheck = yQCDNUM["ICheck"].as<int>(); _splineOrder = yQCDNUM["SplineOrder"].as<int>(); _readTables = yQCDNUM["Read_QCDNUM_Tables"].as<int>(); // Get grids - vector<double> xGrid = getSeq<double>("xGrid"); - vector<int> xGridW = getSeq<int>("xGridW"); - vector<double> Q2Grid = getSeq<double>("Q2Grid"); - vector<double> Q2GridW = getSeq<double>("Q2GridW"); + vector<double> xGrid = getSeq<double>(yQCDNUM["xGrid"]); + vector<int> xGridW = getSeq<int> (yQCDNUM["xGridW"]); + vector<double> Q2Grid = getSeq<double>(yQCDNUM["Q2Grid"]); + vector<double> Q2GridW = getSeq<double>(yQCDNUM["Q2GridW"]); int nxgrid = yQCDNUM["NXbins"].as<int>(); int nxout = 0; @@ -177,25 +178,25 @@ namespace xfitter QCDNUM::dmpwgt(1,22,"unpolarised.wgt"); } - return ; + //Evolution gets its decomposition from YAML + gPdfDecomp=XFITTER_PARS::getInputFunctionFromYaml(yQCDNUM); + atConfigurationChange(); } - // Initialize at - void EvolutionQCDNUM::initAtIteration() + void EvolutionQCDNUM::atConfigurationChange() { - gPdfDecomp = _inPDFs; // XXXXXXXXXXXXXX - const double* q0 = XFITTER_PARS::gParameters.at("Q0"); - int iq0 = QCDNUM::iqfrmq( (*q0) * (*q0) ); const double* Mz = XFITTER_PARS::gParameters.at("Mz"); const double* alphas = XFITTER_PARS::gParameters.at("alphas"); - double epsi = 0; - QCDNUM::setalf(*alphas,(*Mz)*(*Mz)); + } + void EvolutionQCDNUM::atIteration(){ + const double* q0 = XFITTER_PARS::gParameters.at("Q0"); + int iq0 = QCDNUM::iqfrmq( (*q0) * (*q0) ); + double epsi = 0; QCDNUM::evolfg(_itype,funcPDF,qcdnumDef,iq0,epsi); - return ; } std::function<std::map<int,double>(double const& x, double const& Q)> EvolutionQCDNUM::xfxQMap() { @@ -203,8 +204,8 @@ namespace xfitter const auto _f0 = [=] (double const& x, double const& Q) -> std::map<int, double> { std::map<int, double> res; for (int ipdf =-6; ipdf<7; ipdf++) { - int ii = ( ipdf == 0 ) ? 21 : ipdf ; - res[ii] = QCDNUM::fvalxq(_itype,ipdf,x,Q*Q,_icheck); + int ii = ( ipdf == 0 ) ? 21 : ipdf ; + res[ii] = QCDNUM::fvalxq(_itype,ipdf,x,Q*Q,_icheck); } return res; }; diff --git a/include/CheckForPDF.h b/include/CheckForPDF.h index b741706a3e91f9c31ce533585690d20174772cb2..d7ce9c7407071dea6bff69d1be04b858bc91d288 100644 --- a/include/CheckForPDF.h +++ b/include/CheckForPDF.h @@ -1,7 +1,7 @@ #ifndef __CHECKFORPDF_H #define __CHECKFORPDF_H -/// for Fortran calling, check presense of PDF file. +/// for Fortran calling, check presence of PDF file. void CheckForPDF(char const*pdfname); diff --git a/include/ReactionTheory.h b/include/ReactionTheory.h index 59fcd895597b83c836f98711e9d3da05e762d7d2..01facec836b60bafec707daeb3c8c1f37d7305ae 100644 --- a/include/ReactionTheory.h +++ b/include/ReactionTheory.h @@ -24,7 +24,7 @@ typedef double (*pTwoParFunc)(const double&, const double& ); typedef void (*pThreeParSub)(const double& , const double&, const double&); // Function to emulate LHAPDF xfx behavior: -typedef void (*pXFXlike)(const double&, const double&, double*); +typedef void (*pXFXlike)(const double&x,const double&Q,double*results); //using pZeroParFunc = std::function< double() >; //using pOneParFunc = std::function< double(const double&) >; @@ -60,6 +60,7 @@ class ReactionTheory using super = ReactionTheory; virtual string getReactionName() const =0; ///< Should return expected reaction name. Normally generated automatically by AddReaction.py + //A better name would be atStart virtual int initAtStart(const string &) =0; ///< Initialization first time ReactionTheory implementation is called virtual void setxFitterParameters(map<string,double*> &xfitter_pars) {_xfitter_pars = xfitter_pars; }; ///< Set environment map for doubles @@ -81,6 +82,7 @@ class ReactionTheory /// Perform optional re-initialization for a given iteration. Interface for old-style pdf functions + //A better name would be atIteration virtual void initAtIteration(); //! Perform optional action when minuit fcn 3 is called (normally after fit) @@ -128,7 +130,7 @@ class ReactionTheory /// Set evolution name void setEvolution(std::string& evolution) { _evolution = evolution; } - /// Retrieve evolition + /// Retrieve evolution name //A better name would be "getEvolutionName" -- Ivan const std::string getEvolution() const { return _evolution; } protected: diff --git a/include/pdfparam.inc b/include/pdfparam.inc index d83015dc0907afe2eada72848d149ba699d949db..f4b077271d84cb307df746960a75cd2cf83c4e0b 100644 --- a/include/pdfparam.inc +++ b/include/pdfparam.inc @@ -1,6 +1,10 @@ !> PDF parameterisations +C !> Common parameters for sumrules + double precision::uvalSum + double precision::dvalSum + common/sumrule_sums_common/uvalSum,dvalSum C !> Standard parameterisation double precision paruval(1:10), pardval(1:10) diff --git a/include/xfitter_pars.h b/include/xfitter_pars.h index 73f9a56e289a47bc135538d6a8a96d9a05ca9c9a..baba52b5744691fb5839f90e5276184d6c451bbe 100644 --- a/include/xfitter_pars.h +++ b/include/xfitter_pars.h @@ -22,17 +22,23 @@ using std::vector; namespace xfitter{ // to be defined in evolutions/ - class BaseEvolution; + class BaseEvolution; // to be defined in pdfdecompositions/ - class BasePdfDecomposition; + class BasePdfDecomposition; + // to be defined in pdfparams/ + class BasePdfParam; // to be defined in minimizers/ class BaseMinimizer; + using InitialPDFfunction=std::function<std::map<int,double>(const double&x)>; + using EvolvedPDFfunction=std::function<void(double const&x,double const&Q,double*pdfs)>; } namespace XFITTER_PARS { /// Global pointer to the mimimizer extern xfitter::BaseMinimizer* gMinimizer; + /// Globally available YAML node pointing to root of YAML parameters tree, read from parameters.yaml. Might be modified during runtime + extern YAML::Node rootNode; /// Global map of double parameters. They can be used by the minimizer. Initialized based on parameters.yaml extern map<string,double*> gParameters; @@ -49,16 +55,33 @@ namespace XFITTER_PARS { extern map<string,YAML::Node> gParametersY; /// Global map of PDF functions produced by evolutions. - extern map<string,std::function<void(double const& x, double const& Q, double* pdfs)> > gXfxQArrays; + extern map<string,xfitter::EvolvedPDFfunction> gXfxQArrays; /// Global map to store evolutions extern map<string,xfitter::BaseEvolution*> gEvolutions; - - /// Global map to store evolutions + /// Global map to store decompositions extern map<string,xfitter::BasePdfDecomposition*> gPdfDecompositions; + /// Global map to store parameterisations + extern map<string,xfitter::BasePdfParam*>gParameterisations; + /// Helper function to get input function from a yaml node + /// + /// It finds a "decomposition" subnode in given node, extracts a decomposition name from it, finds this decomposition and returns its output function + /// This function will either return a valid function, or issue a fatal error and never return + /// This is used in evolution initialization + xfitter::InitialPDFfunction getInputFunctionFromYaml(const YAML::Node&); + /// Helper function to get a yaml node corresponding to an evolution, by this evolutions's instance name + YAML::Node getEvolutionNode(const std::string&name=""); + /// Helper function to get a yaml node corresponding to a decomposition, by this decomposition's instance name + YAML::Node getDecompositionNode(const std::string&name=""); + /// Helper function to get a yaml node corresponding to a parameterisation, by this parameterisation's instance name + YAML::Node getParameterisationNode(const std::string&name=""); /// Helper function to get string parameters std::string getParameterS(std::string name); + + /// Helper functions + string getDefaultEvolutionName(); + string getDefaultDecompositionName(); /// Parse yaml file @param name void parse_file(const std::string& name); diff --git a/include/xfitter_steer.h b/include/xfitter_steer.h index 28b7d49b7b05ca7eb08cb0c836e9505ba2cea52a..387ea8eb3476b5284ff03f9a38c188d01bd8167c 100644 --- a/include/xfitter_steer.h +++ b/include/xfitter_steer.h @@ -14,12 +14,14 @@ namespace xfitter { class BaseEvolution; class BasePdfDecomposition; + class BasePdfParam; class BaseMinimizer; /// Load named evolution code. BaseEvolution* get_evolution(std::string name=""); /// Load named pdfDecomposition code. BasePdfDecomposition* get_pdfDecomposition(std::string name=""); + BasePdfParam*getParameterisation(const std::string&name=""); /// Load the minimizer BaseMinimizer* get_minimizer(); } diff --git a/minimizers/BaseMinimizer/include/BaseMinimizer.h b/minimizers/BaseMinimizer/include/BaseMinimizer.h index 5002a47f0394f7f51ae3cbb4b9e517042e1ae470..f76069a31256ded0bf2be023dd4e2e90d1a5de70 100644 --- a/minimizers/BaseMinimizer/include/BaseMinimizer.h +++ b/minimizers/BaseMinimizer/include/BaseMinimizer.h @@ -25,7 +25,7 @@ namespace xfitter {} /// Initialization - virtual void initAtStart() = 0; + virtual void atStart() = 0; /// Provide some information virtual void printInfo(){}; @@ -41,7 +41,7 @@ namespace xfitter virtual void addParameter(double par, std::string const &name, double step = 0.01, double const* bounds = nullptr, double const* priors = nullptr ); /// Action at each iteration - virtual void initAtIteration(){}; + virtual void atIteration(){}; /// Miniminzation loop virtual void doMimimization() = 0; diff --git a/minimizers/CERESMinimizer/include/CERESMinimizer.h b/minimizers/CERESMinimizer/include/CERESMinimizer.h index bb16c1695163445fa3042a2808a54d34f0cf9599..fe95f48cfd31925b3fc2e17a019d93468826bf14 100644 --- a/minimizers/CERESMinimizer/include/CERESMinimizer.h +++ b/minimizers/CERESMinimizer/include/CERESMinimizer.h @@ -25,7 +25,7 @@ namespace xfitter { CERESMinimizer (const std::string& inName); /// Optional initialization at the first call - virtual void initAtStart() override final; + virtual void atStart() override final; /// Miniminzation loop virtual void doMimimization() override final; diff --git a/minimizers/CERESMinimizer/src/CERESMinimizer.cc b/minimizers/CERESMinimizer/src/CERESMinimizer.cc index d4b23521c6d958c5c6b6a58b9bc5341d5da6b283..7ae82228bf4b99fccc00ae063792bab0c956ac2a 100644 --- a/minimizers/CERESMinimizer/src/CERESMinimizer.cc +++ b/minimizers/CERESMinimizer/src/CERESMinimizer.cc @@ -81,8 +81,10 @@ CERESMinimizer::CERESMinimizer(const std::string& inName) : BaseMinimizer(inName } // Init at start: -void CERESMinimizer::initAtStart() { - google::InitGoogleLogging(""); +void CERESMinimizer::atStart() { + //There's a suppressed warning here + //warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings] + google::InitGoogleLogging((char*)""); return; } diff --git a/minimizers/MINUITMinimizer/include/MINUITMinimizer.h b/minimizers/MINUITMinimizer/include/MINUITMinimizer.h index 32057c9ec15ff3aedbaf658041db7be1e7d67c30..d34c06f82ab976f7dd4a1cb2708a09d6b8b01b5c 100644 --- a/minimizers/MINUITMinimizer/include/MINUITMinimizer.h +++ b/minimizers/MINUITMinimizer/include/MINUITMinimizer.h @@ -24,7 +24,7 @@ class MINUITMinimizer : public BaseMinimizer MINUITMinimizer (const std::string& inName); /// Optional initialization at the first call - virtual void initAtStart() override final; + virtual void atStart() override final; /// Miniminzation loop virtual void doMimimization() override final; diff --git a/minimizers/MINUITMinimizer/src/MINUITMinimizer.cc b/minimizers/MINUITMinimizer/src/MINUITMinimizer.cc index 53c8c42f9c0a8eed4976207a7352d405a225554a..60bd9928b78e1bfcbdc94218e645ffe4c85dcfb1 100644 --- a/minimizers/MINUITMinimizer/src/MINUITMinimizer.cc +++ b/minimizers/MINUITMinimizer/src/MINUITMinimizer.cc @@ -91,7 +91,7 @@ MINUITMinimizer::MINUITMinimizer(const std::string& inName) : BaseMinimizer(inNa } // Init at start: -void MINUITMinimizer::initAtStart() { +void MINUITMinimizer::atStart() { // Call fortran interface generate_io_filenames_(); minuit_ini_(); diff --git a/parameters.yaml b/parameters.yaml index 70f954334f48cacdfe9c1f3665f1013af93d6f58..28adb9bf56a3882c9cb8f08cd5528bfb0c70ef3a 100644 --- a/parameters.yaml +++ b/parameters.yaml @@ -1,41 +1,88 @@ -# File to store parameters for theory modules. Overwrites -# fortran-based ewparam.txt. - -# -# default minimizer: -# -Minimizer : MINUIT # CERES - -MINUIT : - Commands : | - fix 6 7 9 11 12 13 14 15 16 17 18 19 20 22 +Minimizer: MINUIT # CERES +MINUIT: + Commands: | + fix 6 7 9 11 12 13 14 15 16 17 18 19 20 22 call fcn 3 # doErrors : Hesse # None - - -Profiler: - alphas : [ 0.118, 0.119, 0.117 ] - evolutions : - proton : - sets : [ NNPDF31_nnlo_as_0118, CT10 ] - members : [ [ 0, 1, end ], [0, 1, end] ] - -# -# default evolution: -# -Evolution : LHAPDF - -LHAPDF : - set : "CT10" - member : 0 - - -# -# default PDF decomposition -# -PDFDecomposition : UvDvubardbars -include : "yaml/pdfdecompositions/UvDvUbarDbarS/parameters.yaml" - + #value step min max priors +Parameters: + Ag : DEPENDENT + Bg : [-0.061953] + Cg : [ 5.562367, 0.55] + Auv : DEPENDENT + Buv : [ 0.810476] + Cuv : [ 4.823512] + Duv : [ 0 ] + Euv : [ 9.921366] + Adv : DEPENDENT + Bdv : [ 1.029995] + Cdv : [ 4.846279] + Aubar: [ 0.27 ] + Bubar: [-0.13 ] + Cubar: [ 7.059694] + Dubar: [ 1.548098] + Adbar: [ 0.27 ] + Bdbar: [-0.13 ] + Cdbar: + value: 9.586246 + step: 1.2345 + #min + #max + #pr_mean + #pr_sigma + As : [ 0.13 ] + Bs : [ 0 ] + Cs : [10.0 ] + +Parameterisations: + par_uv: + class: HERAPDF + parameters: [Auv,Buv,Cuv,Duv,Euv] + par_dv: + class: HERAPDF + parameters: [Adv,Bdv,Cdv] + par_ubar: + class: HERAPDF + parameters: [Aubar,Bubar,Cubar,Dubar] + par_dbar: + class: HERAPDF + parameters: [Adbar,Bdbar,Cdbar] + par_s: + class: HERAPDF + parameters: [As,Bs,Cs] + par_g: + class: HERAPDF + parameters: [Ag,Bg,Cg] + +DefaultDecomposition: proton +Decompositions: + proton: + class: UvDvubardbars + xuv: par_uv + xdv: par_dv + xubar: par_ubar + xdbar: par_dbar + xs: par_s + xg: par_g + +DefaultEvolution: proton-QCDNUM +Evolutions: + proton-QCDNUM: + class: QCDNUM + decomposition: proton #this could be omitted, as the default decomposition is set + xGrid : [9.9e-7, 0.01, 0.1, 0.4, 0.7] + xGridW : [1, 2, 4, 8, 16] + Q2Grid : [1., 2.05e8 ] + Q2GridW : [1., 1.] + NQ2bins : 120 + NXbins : 200 + Read_QCDNUM_Tables : 1 + SplineOrder : 2 + ICheck : 0 + proton-LHAPDF: + class: LHAPDF + set: "CT10" + member: 0 # # Initial scale @@ -119,13 +166,13 @@ fcharm: 0. # RT DIS scheme settings: #--------------------------- include : yaml/reactions/RT_DISNC/parameters.yaml - # # FONLL scheme settings: #--------------------------- include : yaml/reactions/FONLL_DISNC/parameters.yaml include : yaml/reactions/FONLL_DISCC/parameters.yaml +# # FF ABM scheme settings: #--------------------------- include : yaml/reactions/FFABM_DISNC/parameters.yaml @@ -136,7 +183,6 @@ include : yaml/reactions/FFABM_DISCC/parameters.yaml #--------------------------- include : yaml/reactions/APPLgrid/parameters.yaml - # # (optional) Fractal module settings: #--------------------------- @@ -146,8 +192,8 @@ include : yaml/reactions/APPLgrid/parameters.yaml # Speficify HF scheme used for DIS NC processes: #--------------------------- hf_scheme_DISNC : - defaultValue : 'BaseDISNC' # global specification -# defaultValue : 'RT_DISNC' # global specification +# defaultValue : 'BaseDISNC' # global specification + defaultValue : 'RT_DISNC' # global specification # defaultValue : 'FONLL_DISNC' # global specification # defaultValue : 'FFABM_DISNC' # 'HERA1+2 NCep 920' : 'BaseDISNC' # datafile specific (based on name) diff --git a/pdfdecompositions/BasePdfDecomposition/include/BasePdfDecomposition.h b/pdfdecompositions/BasePdfDecomposition/include/BasePdfDecomposition.h index 94f3141a7ae1325c238f93b26ebf22631c4894fa..7b99701bc748e3284624018a0b048a2e94053c6f 100644 --- a/pdfdecompositions/BasePdfDecomposition/include/BasePdfDecomposition.h +++ b/pdfdecompositions/BasePdfDecomposition/include/BasePdfDecomposition.h @@ -21,45 +21,35 @@ namespace xfitter { class BasePdfDecomposition { public: + const std::string _name;//Unique name used to identify this decomposition instance - /// Default constructor. Name is the PDF name - BasePdfDecomposition(const std::string& inName): _name(inName) { }; - virtual ~BasePdfDecomposition() {}; + /// Default constructor. + BasePdfDecomposition(const char*name):_name(name){} + virtual ~BasePdfDecomposition(){} /// Initialization at the first call - virtual void initAtStart(const std::string& pars) = 0; - + virtual void atStart(){} /// Optional initialization at each iteration. Can be used to compute sum-rules - virtual void initAtIteration() {} + virtual void atIteration(){} + /// This function should be called when at least one parameter in the YAML node of given decomposition changes + virtual void atConfigurationChange(){} /// Print pdf parameters - virtual void printParams() {} - + //This shouldn't be here, printing parameters should be just a global function --Ivan + virtual void printParams(){} + /// Returns a LHAPDF-style function, that returns PDFs in a physical basis for given x virtual std::function<std::map<int,double>(const double& x)> f0() const = 0; - - void addParameterisation(const std::string& pname, BasePdfParam* pdfParam) { - _pdfParams[pname] = pdfParam; - } - - BasePdfParam* getPdfParam(std::string const& name) const { - return _pdfParams.at(name); - } - - const std::string& getName() const { return _name; } + /// Get class name, can be used to verify that the correct concrete class is being used + virtual const char*getClassName()const=0; protected: /// PDF parameterisations + //Not really needed in this form --Ivan std::map<std::string,BasePdfParam*> _pdfParams; - - - - private: - /// Name of PDF decomposition - std::string _name; - }; /// For dynamic loader + //Wait, is this even used anywhere? --Ivan typedef BasePdfDecomposition* create_pdfDecomposition(); } diff --git a/pdfdecompositions/BasePdfDecomposition/src/BasePdfDecomposition.cc b/pdfdecompositions/BasePdfDecomposition/src/BasePdfDecomposition.cc index 4f9c53666e90961d8b6710106c9c7ddb5c88b7d8..b835543f822f968e1af1ce927b9fb35636599fa8 100644 --- a/pdfdecompositions/BasePdfDecomposition/src/BasePdfDecomposition.cc +++ b/pdfdecompositions/BasePdfDecomposition/src/BasePdfDecomposition.cc @@ -1,3 +1 @@ #include "BasePdfDecomposition.h" - - diff --git a/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h b/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h deleted file mode 100644 index 465644ff1cc2725821c3508752d433b5865819c7..0000000000000000000000000000000000000000 --- a/pdfdecompositions/GRV_PionPdfDecomposition/include/GRV_PionPdfDecomposition.h +++ /dev/null @@ -1,39 +0,0 @@ -#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 deleted file mode 100644 index 8cad3faa5d44c1c559b32d840025275454314f80..0000000000000000000000000000000000000000 --- a/pdfdecompositions/GRV_PionPdfDecomposition/src/GRV_PionPdfDecomposition.cc +++ /dev/null @@ -1,98 +0,0 @@ -/* - @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 deleted file mode 100644 index 49ee9e52af54f000e966a81084ed8d4d623985c2..0000000000000000000000000000000000000000 --- a/pdfdecompositions/GRV_PionPdfDecomposition/src/Makefile.am +++ /dev/null @@ -1,12 +0,0 @@ - -# 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/LHAPDFDecomposition/include/LHAPDFDecomposition.h b/pdfdecompositions/LHAPDFDecomposition/include/LHAPDFDecomposition.h index 83c2d2bec32ea6e3fced4f6135ff69fbc07d3adc..0700d5d96660bb03d0f05fda8c074ac072247acf 100644 --- a/pdfdecompositions/LHAPDFDecomposition/include/LHAPDFDecomposition.h +++ b/pdfdecompositions/LHAPDFDecomposition/include/LHAPDFDecomposition.h @@ -2,8 +2,13 @@ #pragma once #include "BasePdfDecomposition.h" - +//Try to suppress unused-local-typedef warning from boost 1.53.0 for gcc +//Apparently these warnings have been fixed in later versions of boost +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-local-typedefs" #include "LHAPDF/LHAPDF.h" +#pragma GCC diagnostic pop + /** @class LHAPDFDecomposition @@ -19,17 +24,17 @@ namespace xfitter class LHAPDFDecomposition : public BasePdfDecomposition { public: - /// Default constructor. Name is the PDF name - LHAPDFDecomposition(const std::string& PDFset, const int& mem = 0); + LHAPDFDecomposition(const char*name); + ~LHAPDFDecomposition(); + virtual const char*getClassName()const override final; /// Optional initialization at the first call - virtual void initAtStart(const std::string& pars) override final; + virtual void atStart()override final; - /// Compute PDF in a physical base in LHAPDF format at the initial scale + /// Compute PDF in a physical basis in LHAPDF format at the initial scale virtual std::function<std::map<int,double>(const double& x)> f0() const override final; private: - const int _mem; - std::vector<LHAPDF::PDF*> _dist; + LHAPDF::PDF*_pdf{nullptr}; }; } diff --git a/pdfdecompositions/LHAPDFDecomposition/src/LHAPDFDecomposition.cc b/pdfdecompositions/LHAPDFDecomposition/src/LHAPDFDecomposition.cc index 00ceaac8519030d087b01233f6a80055253a97ed..946385d1a6657f3d21706b9df05bb7f3c671d744 100644 --- a/pdfdecompositions/LHAPDFDecomposition/src/LHAPDFDecomposition.cc +++ b/pdfdecompositions/LHAPDFDecomposition/src/LHAPDFDecomposition.cc @@ -6,23 +6,36 @@ Created by AddPdfDecomposition.py on 2018-07-12 */ -#include "LHAPDFDecomposition.h" -#include "xfitter_pars.h" +#include"LHAPDFDecomposition.h" +#include"xfitter_pars.h" +#include"xfitter_cpp_base.h" +#include"CheckForPDF.h" -namespace xfitter -{ +namespace xfitter { + //For dynamic loading: + extern "C" LHAPDFDecomposition*create(const char*name){ + return new LHAPDFDecomposition(name); + } //_________________________________________________________________________________ - LHAPDFDecomposition::LHAPDFDecomposition(const std::string& PDFset, const int& mem): - BasePdfDecomposition{"LHAPDF"}, - _mem(mem) - { - // Upload PDF set from LHAPDF - _dist = LHAPDF::mkPDFs(PDFset); + LHAPDFDecomposition::LHAPDFDecomposition(const char*name):BasePdfDecomposition{name}{} + LHAPDFDecomposition::~LHAPDFDecomposition(){if(_pdf)delete _pdf;} + const char*LHAPDFDecomposition::getClassName()const{return"LHAPDF";} + void LHAPDFDecomposition::atStart(){ + YAML::Node pars=XFITTER_PARS::getDecompositionNode(_name); + string setName; + int member; + try{ + setName=pars["set"].as<std::string>(); + }catch(YAML::TypedBadConversion<std::string>&ex){ + hf_errlog(18090310,"F: In LHAPDFDecomposition::atStart: failed to convert YAML node \"set\" to string; printing node to stderr"); + std::cerr<<pars<<std::endl; } - //_________________________________________________________________________________ - void LHAPDFDecomposition::initAtStart(const std::string & pars) - { + // check if exists first + CheckForPDF(setName.c_str()); + member=pars["member"].as<int>(); + if(_pdf)delete _pdf; + _pdf=LHAPDF::mkPDF(setName,member); } //_________________________________________________________________________________ @@ -30,7 +43,7 @@ namespace xfitter { return [=] (const double& x)->std::map<int,double> { - return _dist[_mem]->xfxQ(x, *(XFITTER_PARS::gParameters.at("Q0"))); + return _pdf->xfxQ(x, *(XFITTER_PARS::gParameters.at("Q0"))); }; } } diff --git a/pdfdecompositions/SU3_PionPdfDecomposition/include/SU3_PionPdfDecomposition.h b/pdfdecompositions/SU3_PionPdfDecomposition/include/SU3_PionPdfDecomposition.h new file mode 100644 index 0000000000000000000000000000000000000000..010ba1c4e8cd5f515700266d3b0a04ef47b2fe7a --- /dev/null +++ b/pdfdecompositions/SU3_PionPdfDecomposition/include/SU3_PionPdfDecomposition.h @@ -0,0 +1,38 @@ +#pragma once +#include "BasePdfDecomposition.h" +/** + @class SU3_PionPdfDecomposition + + @brief A class for pdf decomposition for the pion with SU3-symmetric sea + + Used for pi- + Assumes that at starting scale: + ubar=d + dbar=u=s=sbar + Parametrised distributions are: + v:=(dval-uval)/2=d-u + S:=(u +dbar)/2=u + g + Therefore, transformations to physical basis: + d=ubar=v+S + u=dbar=s=sbar=S + g=g + others=0 + And sum rules for pi- are: + \int_0^1 v dx=1 + \int_0^1 x*(6S+2v+g) dx=1 + @version 0.2 + @date 2018-08-14 + */ +namespace xfitter{ +class SU3_PionPdfDecomposition:public BasePdfDecomposition{ + public: + SU3_PionPdfDecomposition(const char*name); + virtual const char*getClassName()const override final; + virtual void atStart()override final; + virtual void atIteration()override final; + virtual std::function<std::map<int,double>(const double& x)>f0()const override final; + private: + BasePdfParam*par_v{nullptr},*par_S{nullptr},*par_g{nullptr}; + }; +} diff --git a/pdfdecompositions/SU3_PionPdfDecomposition/src/Makefile.am b/pdfdecompositions/SU3_PionPdfDecomposition/src/Makefile.am new file mode 100644 index 0000000000000000000000000000000000000000..fb7bcd83cbc4215307e0e3736d002fc4d23e15e2 --- /dev/null +++ b/pdfdecompositions/SU3_PionPdfDecomposition/src/Makefile.am @@ -0,0 +1,16 @@ + +# 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 \ + -I$(srcdir)/../../../pdfparams/PolySqrtPdfParam/include/ #Hacks - remove later + +lib_LTLIBRARIES=libSU3_PionPdfDecomposition_xfitter.la +libSU3_PionPdfDecomposition_xfitter_la_SOURCES = SU3_PionPdfDecomposition.cc + +datadir = ${prefix}/yaml/pdfdecompositions/SU3_Pion +data_DATA = ../yaml/parameters.yaml + +dist_noinst_HEADERS = ../include ../yaml + +# XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXx These hacks are temporary +libSU3_PionPdfDecomposition_xfitter_la_LDFLAGS=-lPolySqrtPdfParam_xfitter -lBasePdfParam_xfitter -L$(libdir) diff --git a/pdfdecompositions/SU3_PionPdfDecomposition/src/SU3_PionPdfDecomposition.cc b/pdfdecompositions/SU3_PionPdfDecomposition/src/SU3_PionPdfDecomposition.cc new file mode 100644 index 0000000000000000000000000000000000000000..389d0de2513e57f3f56447f8f8445a8590833c36 --- /dev/null +++ b/pdfdecompositions/SU3_PionPdfDecomposition/src/SU3_PionPdfDecomposition.cc @@ -0,0 +1,75 @@ +/* + @file SU3_PionPdfDecomposition.cc + @date 2018-08-07 + @author AddPdfDecomposition.py + Created by AddPdfDecomposition.py on SU3_Pion +*/ +#include"xfitter_cpp_base.h" +#include<iostream> +#include"SU3_PionPdfDecomposition.h" +#include"xfitter_pars.h" +#include"xfitter_steer.h" +using uint=unsigned int; +using namespace std; +//For dynamic loading: +namespace xfitter{ +extern "C" SU3_PionPdfDecomposition*create(const char*name){ + return new SU3_PionPdfDecomposition(name); +} +SU3_PionPdfDecomposition::SU3_PionPdfDecomposition(const char*name):BasePdfDecomposition{name}{} +const char*SU3_PionPdfDecomposition::getClassName()const{return"SU3_Pion";} +BasePdfParam*getParam(const BasePdfDecomposition*self,const YAML::Node&node,const char*s){ + try{ + return getParameterisation(node[s].as<string>()); + }catch(const YAML::InvalidNode&ex){ + if(node[s].IsNull()){ + cerr<<"[ERROR] No \""<<s<<"\" parameterisation given for decomposition \""<<self->_name<<"\""<<endl; + hf_errlog(18092410,"F: Error in decomposition parameters, details written to stderr"); + }else throw ex; + }catch(const YAML::BadConversion&ex){ + cerr<<"[ERROR] Bad parameter \""<<s<<"\" given for decomposition \""<<self->_name<<"\""<<endl; + hf_errlog(18092410,"F: Error in decomposition parameters, details written to stderr"); + } + return nullptr;//unreachable, suppress warning +} +// Init at start: +void SU3_PionPdfDecomposition::atStart(){ + const YAML::Node node=XFITTER_PARS::getDecompositionNode(_name); + //get parameterisation usually doesn't throw + par_v=getParam(this,node,"valence"); + par_S=getParam(this,node,"sea"); + par_g=getParam(this,node,"gluon"); +} +void SU3_PionPdfDecomposition::atIteration() { + //Enforce sum rules + //Valence sum + par_v->setMoment(-1,1); + //Momentum sum + par_g->setMoment(0,1-6*par_S->moment(0)-2*par_v->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)>SU3_PionPdfDecomposition::f0()const{ + return [=](double const& x)->std::map<int, double>{ + double v=(*par_v)(x); + double S=(*par_S)(x); + double g=(*par_g)(x); + double d=S+v; + std::map<int,double>res_={ + {-6,0}, + {-5,0}, + {-4,0}, + {-3,S},//sbar + {-2,d},//ubar + {-1,S},//dbar + { 1,d},//d + { 2,S},//u + { 3,S},//s + { 4,0}, + { 5,0}, + { 6,0}, + {21,g} + }; + return res_; + }; +} +} diff --git a/pdfdecompositions/GRV_PionPdfDecomposition/yaml/parameters.yaml b/pdfdecompositions/SU3_PionPdfDecomposition/yaml/parameters.yaml similarity index 100% rename from pdfdecompositions/GRV_PionPdfDecomposition/yaml/parameters.yaml rename to pdfdecompositions/SU3_PionPdfDecomposition/yaml/parameters.yaml diff --git a/pdfdecompositions/UvDvUbarDbarS/include/UvDvUbarDbarS.h b/pdfdecompositions/UvDvUbarDbarS/include/UvDvUbarDbarS.h index 74e06bcc22460d09d7ba99308483211d22b78f1a..ae93e7c953665fcbfd9e18f71218f0073a75c514 100644 --- a/pdfdecompositions/UvDvUbarDbarS/include/UvDvUbarDbarS.h +++ b/pdfdecompositions/UvDvUbarDbarS/include/UvDvUbarDbarS.h @@ -20,13 +20,14 @@ namespace xfitter { public: /// Default constructor. Name is the PDF name - UvDvUbarDbarS(); + UvDvUbarDbarS(const char*name); + virtual const char*getClassName()const override final; /// Optional initialization at the first call - virtual void initAtStart(const std::string& pars) override final; + virtual void atStart()override final; /// Compute sum-rules - virtual void initAtIteration() override final; + virtual void atIteration() override final; /// print parameters virtual void printParams() override final; @@ -36,31 +37,11 @@ namespace xfitter private: - /// Get parameter values from the minimizer - std::unique_ptr<double[]> getParValues(BasePdfParam const* param) const; - - /// Get valence - double valence(double x,const std::string&name)const; - - /// Get sea - double sea(double x, std::string const& name) const; - - /// Get uv, apply sum-rule - double uv(double x) const; - - /// Get dv, apply sum-rule - double dv(double x) const; - - /// Get dbar - double dbar(double x) const; - - /// Get ubar - double ubar(double x) const; - - /// Get s - double s(double x) const; - - /// Get g - double g(double x) const; + BasePdfParam*par_xuv{nullptr}, + *par_xdv{nullptr}, + *par_xubar{nullptr}, + *par_xdbar{nullptr}, + *par_xs{nullptr}, + *par_xg{nullptr}; }; } diff --git a/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc b/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc index 6ca48015e4d6d9414e4609e9454fe1c7592cef04..630f3aed2ae931431f5c6b6c5c27c5d93d65c442 100644 --- a/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc +++ b/pdfdecompositions/UvDvUbarDbarS/src/UvDvUbarDbarS.cc @@ -7,8 +7,8 @@ */ #include "UvDvUbarDbarS.h" -#include "HERAPDF_PdfParam.h" #include "xfitter_pars.h" +#include "xfitter_steer.h" #include <iostream> #include <iomanip> #include <cmath> @@ -16,135 +16,98 @@ namespace xfitter { /// the class factories - extern "C" UvDvUbarDbarS* create() { - return new UvDvUbarDbarS(); + extern "C" UvDvUbarDbarS*create(const char*name){ + return new UvDvUbarDbarS(name); } //_________________________________________________________________________________ - UvDvUbarDbarS::UvDvUbarDbarS(): BasePdfDecomposition{"UvDvUbarDbarS"} { } + UvDvUbarDbarS::UvDvUbarDbarS(const char*name):BasePdfDecomposition{name}{} + const char*UvDvUbarDbarS::getClassName()const{return"UvDvUbarDbarS";} //_________________________________________________________________________________ - 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); - pParam->initFromYaml(pdf.second); - addParameterisation(pdfName,pParam); - } + void UvDvUbarDbarS::atStart(){ + const YAML::Node node=XFITTER_PARS::getDecompositionNode(_name); + //TODO: handle errors + par_xuv =getParameterisation(node["xuv"].as<string>()); + par_xdv =getParameterisation(node["xdv"].as<string>()); + par_xubar=getParameterisation(node["xubar"].as<string>()); + par_xdbar=getParameterisation(node["xdbar"].as<string>()); + par_xs =getParameterisation(node["xs"].as<string>()); + par_xg =getParameterisation(node["xg"].as<string>()); } - void UvDvUbarDbarS::initAtIteration() { + void UvDvUbarDbarS::atIteration() { //Enforce sum rules // counting sum-rules for uv and dv - getPdfParam("xuv")->setMoment(-1,2.0); - getPdfParam("xdv")->setMoment(-1,1.0); + par_xuv->setMoment(-1,2.0); + par_xdv->setMoment(-1,1.0); // momentum sum-rule // quark part 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); + xsumq+= par_xuv ->moment(0); + xsumq+= par_xdv ->moment(0); + xsumq+=2*par_xubar->moment(0); + xsumq+=2*par_xdbar->moment(0); + xsumq+=2*par_xs ->moment(0); // gluon part - getPdfParam("xg")->setMoment(0,1-xsumq); + par_xg->setMoment(0,1-xsumq); printParams(); } void UvDvUbarDbarS::printParams() { - std::cout << "\n" << std::left<< std::setw(8) << " Name " << std::right; - for ( int i =0 ; i<6; i++) { - std::cout << std::setw(12) << " par"+std::to_string(i) ; - } - std::cout << "\n"; - for ( auto p : _pdfParams) { - auto pdfParam = p.second; - auto name = p.first; - auto pars = getParValues(pdfParam); - int npar = pdfParam->getNPar(); - - // 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 << "\n"; - } - std::cout << "\n"; - } - - 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=param->getName()+ "_p" + std::to_string(i) ; - pars[i] = *XFITTER_PARS::gParameters.at(pnam); - } - return pars; - } - - double UvDvUbarDbarS::valence(double x, const std::string&name) const { - return(*getPdfParam(name))(x); - } - - double UvDvUbarDbarS::sea(double x,const std::string&name)const{ - return(*getPdfParam(name))(x); - } - - double UvDvUbarDbarS::uv(double x) const { - return valence(x,"xuv"); - } - - double UvDvUbarDbarS::dv(double x) const { - return valence(x,"xdv"); - } - - double UvDvUbarDbarS::ubar(double x) const { - return sea(x,"xubar"); + //Sorry, I broke this during a rewrite + //This should be a global function anyway + //--Ivan + //std::cout << "\n" << std::left<< std::setw(8) << " Name " << std::right; + //for ( int i =0 ; i<6; i++) { + //std::cout << std::setw(12) << " par"+std::to_string(i) ; + //} + //std::cout << "\n"; + //for ( auto p : _pdfParams) { + //auto pdfParam = p.second; + //auto name = p.first; + //auto pars = getParValues(pdfParam); + //int npar = pdfParam->getNPar(); + + //// 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 << "\n"; + //} + //std::cout << "\n"; } - - double UvDvUbarDbarS::dbar(double x) const { - return sea(x,"xdbar"); // XXXXXXXXXXXXX HARDWIRE - } - - double UvDvUbarDbarS::s(double x) const { - return sea(x,"xs"); - } - - double UvDvUbarDbarS::g(double x)const{ - return (*getPdfParam("xg"))(x); - } - - //_________________________________________________________________________________ std::function<std::map<int,double>(const double& x)> UvDvUbarDbarS::f0() const { // lambda function - const auto _f0 = [=] (double const& x)->std::map<int, double> { + return [=] (double const& x)->std::map<int, double> { + double ubar=(*par_xubar)(x); + double dbar=(*par_xdbar)(x); + double u=(*par_xuv)(x)+ubar; + double d=(*par_xdv)(x)+dbar; + double s=(*par_xs)(x); + double g=(*par_xg)(x); std::map<int, double> res = { {-6,0}, {-5,0}, {-4,0}, - {-3,s(x)}, - {-2,ubar(x)}, - {-1,dbar(x)}, - { 1,dbar(x)+dv(x)}, - { 2,ubar(x)+uv(x)}, - { 3,s(x)}, + {-3,s}, + {-2,ubar}, + {-1,dbar}, + { 1,d}, + { 2,u}, + { 3,s}, { 4,0}, { 5,0}, { 6,0}, - {21,g(x)} + {21,g} }; return res; }; - return _f0; } } diff --git a/pdfparams/BasePdfParam/include/BasePdfParam.h b/pdfparams/BasePdfParam/include/BasePdfParam.h index fa279bfdcf1cda4e8808c9ed33a283eda6ddb77c..2d4a8217073f60f921fa064d891e2aebd749d012 100644 --- a/pdfparams/BasePdfParam/include/BasePdfParam.h +++ b/pdfparams/BasePdfParam/include/BasePdfParam.h @@ -24,10 +24,11 @@ @date 2018-08-13 */ +namespace xfitter{ class BasePdfParam{ public: - BasePdfParam(const std::string&instance_name):_name(instance_name),pars{nullptr},Npars(0){} - ~BasePdfParam(){if(pars)delete[]pars;} + BasePdfParam(const std::string&instance_name):_name(instance_name){} + virtual~BasePdfParam(); 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 @@ -41,14 +42,13 @@ public: 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); + virtual void atStart(); protected: //!Unique name of instance - std::string _name; + const std::string _name; //!Array of pointers to some global locations where minimization parameters are stored - //TODO: Implement proper pars initialization from YAML - double**pars; + double**pars{nullptr}; //!Number of parameters, which is also the size of the array **parameters defined above - unsigned int Npars; + unsigned int Npars{0}; }; +} diff --git a/pdfparams/BasePdfParam/src/BasePdfParam.cc b/pdfparams/BasePdfParam/src/BasePdfParam.cc index ec680c177f7f8c916198ab2f2540ad08112a9bbf..468cf3883d8eceb5d27306d363f5fb8c50a4e34c 100644 --- a/pdfparams/BasePdfParam/src/BasePdfParam.cc +++ b/pdfparams/BasePdfParam/src/BasePdfParam.cc @@ -2,6 +2,7 @@ #include"BaseMinimizer.h" #include"xfitter_pars.h" #include"xfitter_steer.h" +#include"xfitter_cpp_base.h" #include<cmath> #include<memory> #include<iostream> @@ -17,7 +18,9 @@ extern "C" { map<std::string,double*> *map, int len); } +namespace xfitter{ /// Implement numeric integration +BasePdfParam::~BasePdfParam(){if(pars)delete[]pars;} double BasePdfParam::moment(int iMoment)const{ /// Simple rule, split log/lin spacing at xsplit=0.1 @@ -55,15 +58,37 @@ 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 +void BasePdfParam::atStart(){ using namespace std; + YAML::Node node=XFITTER_PARS::getParameterisationNode(_name); + YAML::Node parsNode=node["parameters"]; + if(!parsNode.IsSequence()){ + cerr<<"[ERROR] Bad \"parameters\" for parameterisation \""<<_name<<"\", expected sequence"<<endl; + hf_errlog(18092420,"F: Bad parameters in parameterisation, see stderr"); + } + Npars=parsNode.size(); + pars=new double*[Npars]; + //TODO: destructor + for(unsigned int i=0;i<Npars;++i){ + try{ + pars[i]=XFITTER_PARS::gParameters.at(parsNode[i].as<string>()); + }catch(const YAML::BadConversion&ex){ + cerr<<"[ERROR] Bad name of parameter "<<i<<" for parameterisation \""<<_name<<"\": "<<ex.what()<<endl; + hf_errlog(18092421,"F: Bad parameter name in parameterisation definition, see stderr"); + }catch(const out_of_range&ex){ + string parname=parsNode[i].as<string>(); + if(XFITTER_PARS::gParameters.count(parname)==0){ + cerr<<"[ERROR] Unknown parameter \""<<parname<<"\" in parameterisation \""<<_name<<"\""<<endl; + hf_errlog(18092422,"F: Unknown parameter in parameterisation definition, see stderr"); + } + } + } + /* using uint=unsigned int; - cout<<"DEBUG["<<_name<<"]: initFromYaml: value="<<value<<endl; + //cout<<"DEBUG["<<_name<<"]: initFromYaml: value="<<value<<endl; if(value.IsSequence()){ Npars=value.size(); - cout<<Npars<<endl; + //cout<<Npars<<endl; pars=new double*[Npars]; // HARDWIRE old-way for now: XXXXXXXXXXXXXXXXXXXXXXXXXXXXX for(uint i=0;i<Npars;++i){ @@ -78,15 +103,10 @@ void BasePdfParam::initFromYaml(YAML::Node value) { //int add = true; // addexternalparam_(pnam.c_str(),val,step,minv,maxv,priorVal,priorUnc,add,&XFITTER_PARS::gParameters,pnam.size()); - xfitter::BaseMinimizer* minimizer = xfitter::get_minimizer(); - minimizer->addParameter(val,pnam,step,nullptr,nullptr); - - pars[i]=XFITTER_PARS::gParameters.at(pnam); - cout<<pnam<<"="<<(*pars[i])<<endl; } }else{ cout<<"ERROR["<<_name<<"]: initFromYaml: parameter is not a sequence!"<<endl; } - //TODO: Handle possible errors + */ +} } - diff --git a/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h b/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h index f16ce7f9e0d9a093f5743ff675f60722e52ed5e9..b406f9bddb5bb4acd26808b436d30ca6317843de 100644 --- a/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h +++ b/pdfparams/HERAPDF_PdfParam/include/HERAPDF_PdfParam.h @@ -21,8 +21,11 @@ @date 2018-08-14 */ +namespace xfitter{ class HERAPDF_PdfParam:public BasePdfParam{ public: HERAPDF_PdfParam(const std::string&inName):BasePdfParam(inName){} virtual double operator()(double x)const override final; + virtual double moment(int nMoment=-1)const override final; }; +} diff --git a/pdfparams/HERAPDF_PdfParam/src/HERAPDF_PdfParam.cc b/pdfparams/HERAPDF_PdfParam/src/HERAPDF_PdfParam.cc index dd9336dd1d4fd7cddba0ecc76f071373c48a11d7..44ea567e620e1b7a6f539c84ca5799dbfe65bb33 100644 --- a/pdfparams/HERAPDF_PdfParam/src/HERAPDF_PdfParam.cc +++ b/pdfparams/HERAPDF_PdfParam/src/HERAPDF_PdfParam.cc @@ -8,6 +8,11 @@ #include "HERAPDF_PdfParam.h" #include <cmath> +namespace xfitter{ +//for dynamic loading +extern"C" HERAPDF_PdfParam*create(const char*name){ + return new HERAPDF_PdfParam(name); +} // Main function to compute PDF double HERAPDF_PdfParam::operator()(double x)const{ const int npar = getNPar(); @@ -23,11 +28,32 @@ double HERAPDF_PdfParam::operator()(double x)const{ } return power*poly; } - -// Optional function to compute integrals: - -// double HERAPDF_PdfParam::moment( double const* pars, int const iMoment) -// { -// return 0 -// } - +double HERAPDF_PdfParam::moment(int n)const{ + //Integral of HERAPDF-style function is expressed in terms of Euler beta function: + //beta(x,y)=int_0^1 t^(x-1)*(1-t)^(y-1) dx=gamma(x)*gamma(y)/gamma(x+y) + //moment(n)=int_0^1 P[0]*x^(P[1]+n)*(1-x)^P[2]*(1+sum_{i=3}^N{P[i]*x^(i-2)}) + //Let A:=P[0], B:=P[1]+n+1, C:=P[2]+1, then + //moment=int_0^1 A*x^(B-1)*(1-x)^(C-1)*(1+sum_{i=3}^N{P[i]*x^(i-2)}) + //=A*(beta(B,C)+sum_{i=3}^N{P[i]*beta(B+i-2,C)}) + //beta(B+1,C)=B/(B+C) + //=> beta(B+n,C)=beta(B,C)*product_{k=0}^{n-1}{(B+k)/(B+C+k)} + //=> beta(B+i-2,C)=beta(B,C)*product_{k=0}^{i-3}{(B+k)/(B+C+k)} + //=> moment=A*beta(B,C)*(1+sum_{i=3}^N{P[i]*product_{k=0}^{i-3}{(B+k)/(B+C+k)}})= + //=> moment=A*beta(B,C)*(1+P[3]*B/(B+C)+P[4]*B/(B+C)*(B+1)/(B+C+1)+...) + //beta(B,C)=exp(lgamma(B)+lgamma(C)-lgamma(B+C)) + using uint=unsigned int; + const double B=(*pars[1])+(n+1),C=(*pars[2])+1; + const uint N=getNPar(); + double sum=1; + double prod=1; + double a=B; + double b=B+C; + for(uint i=3;i<N;++i){ + prod=prod*a/b; + sum+=(*pars[i])*prod; + a++; + b++; + } + return (*pars[0])*exp(lgamma(B)+lgamma(C)-lgamma(B+C))*sum; +} +} diff --git a/pdfparams/HERAPDF_PdfParam/src/Makefile.am b/pdfparams/HERAPDF_PdfParam/src/Makefile.am index 7d0f214503f97facffaa045d2251977c7a9217c5..2913dd7129ce3c150d59c6abbaa87066ffd83560 100644 --- a/pdfparams/HERAPDF_PdfParam/src/Makefile.am +++ b/pdfparams/HERAPDF_PdfParam/src/Makefile.am @@ -10,3 +10,4 @@ datadir = ${prefix}/yaml/pdfparams/HERAPDF_PdfParam data_DATA = ../yaml/parameters.yaml dist_noinst_HEADERS = ../include ../yaml +libHERAPDF_PdfParam_xfitter_la_LDFLAGS=-lBasePdfParam_xfitter -L$(libdir) diff --git a/pdfparams/PolySqrtPdfParam/include/PolySqrtPdfParam.h b/pdfparams/PolySqrtPdfParam/include/PolySqrtPdfParam.h new file mode 100644 index 0000000000000000000000000000000000000000..c8073b1ab3487ff0a1ed85eb1244bd910f39dac5 --- /dev/null +++ b/pdfparams/PolySqrtPdfParam/include/PolySqrtPdfParam.h @@ -0,0 +1,36 @@ + +#pragma once + +#include "BasePdfParam.h" + +/** + @class PolySqrtPdfParam + + @brief A class for PolySqrt pdf parameterisation + + xf(x)=A*x^B*(1-x)^C*(1+P(sqrt(x))) + where P(sqrt(x)) is a polynomial of sqrt(x) with other parameters as coefficients: + P(sqrt(x))=D*sqrt(x)+E*x+F*x*sqrt(x)+... + + 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-1)}) + + @version 0.1 + @date 2018-08-16 + */ + +namespace xfitter{ +class PolySqrtPdfParam:public BasePdfParam{ + public: + PolySqrtPdfParam(const std::string&inName):BasePdfParam(inName){} + virtual double operator()(double x)const override final; + 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; +}; +} diff --git a/pdfparams/PolySqrtPdfParam/src/Makefile.am b/pdfparams/PolySqrtPdfParam/src/Makefile.am new file mode 100644 index 0000000000000000000000000000000000000000..6a51b542c004cb852176601c1d7680abb6b1716a --- /dev/null +++ b/pdfparams/PolySqrtPdfParam/src/Makefile.am @@ -0,0 +1,13 @@ + +# Created by AddPdfParam.py on 2018-08-16 + +AM_CXXFLAGS = -I$(srcdir)/../include -I$(srcdir)/../../../include -I$(srcdir)/../../BasePdfParam/include -Wall -fPIC -Wno-deprecated + +lib_LTLIBRARIES = libPolySqrtPdfParam_xfitter.la +libPolySqrtPdfParam_xfitter_la_SOURCES = PolySqrtPdfParam.cc + +datadir = ${prefix}/yaml/pdfparams/PolySqrt +data_DATA = ../yaml/parameters.yaml + +dist_noinst_HEADERS = ../include ../yaml +libPolySqrtPdfParam_xfitter_la_LDFLAGS=-lBasePdfParam_xfitter -L$(libdir) diff --git a/pdfparams/PolySqrtPdfParam/src/PolySqrtPdfParam.cc b/pdfparams/PolySqrtPdfParam/src/PolySqrtPdfParam.cc new file mode 100644 index 0000000000000000000000000000000000000000..e1258bd7cd1eef20594d8b06974d0d3210d02050 --- /dev/null +++ b/pdfparams/PolySqrtPdfParam/src/PolySqrtPdfParam.cc @@ -0,0 +1,72 @@ + +/* + @file PolySqrtPdfParam.cc + @date 2018-08-16 + @author AddPdfParam.py + Created by AddPdfParam.py on 2018-08-16 +*/ + +#include"PolySqrtPdfParam.h" +#include<cmath> +using namespace std; +using uint=unsigned int; +namespace xfitter{ +//for dynamic loading +extern"C" PolySqrtPdfParam*create(const char*name){ + return new PolySqrtPdfParam(name); +} +double PolySqrtPdfParam::operator()(double x)const{ + const uint N=getNPar(); + double pol=1; + double mulx=1; + double sqrtx=sqrt(x); + for(uint i=3;i<N;++i){ + mulx*=sqrtx; + pol+=(*pars[i])*mulx; + } + return(*pars[0])*pow(x,(*pars[1]))*pow((1-x),(*pars[2]))*pol; +} +double PolySqrtPdfParam::moment(int n)const{ + //Integral of PolySqrtPdfParam-style function is expressed in terms of Euler beta function: + //beta(x,y)=int_0^1 t^(x-1)*(1-t)^(y-1) dx=gamma(x)*gamma(y)/gamma(x+y) + //beta(B,C)=exp(lgamma(B)+lgamma(C)-lgamma(B+C)) + //moment(n)=int_0^1 P[0]*x^(P[1]+n)*(1-x)^P[2]*(1+sum_{i=3}^N{P[i]*x^(i/2-1)}) + //Let A:=P[0], B:=P[1]+n+1, C:=P[2]+1, then + //moment=int_0^1 A*x^(B-1)*(1-x)^(C-1)*(1+sum_{i=3}^N{P[i]*x^(i/2-1)}) + //=A*(beta(B,C)+sum_{i=3}^N{P[i]*beta(B+i/2-1,C)})= + //=A*(beta(B,C)+sum(i=3;i<=N;i+=2){P[i]*beta(B+i/2-1,C)}+sum(i=4;i<=N;i+=2){P[i]*beta(B+i/2-1,C)}) + //=A*(beta(B,C)+sum(k=0;3+2k<=N;k++){P[3+2k]*beta(B+k+1/2,C)}+sum(k=1;2+2k<=N;k++){P[2+2k]*beta(B+k,C)}) + //beta(B+1,C)=B/(B+C) + //=> beta(B+k,C)=beta(B,C)*product_{i=0}^{k-1}{(B+i)/(B+C+i)} + //=> beta(B+k+1/2,C)=beta(B+1/2,C)*product_{i=0}^{k-1}{(B+1/2+i)/(B+C+1/2+i)} + //=> moment=A*(beta(B,C)*(1+sum(k=1;2+2k<=N;k++){P[2+2k]*product_{i=0}^{k-1}{(B+i)/(B+C+i)}})+beta(B+1/2,C)*sum(k=0;3+2k<=N;k++){P[3+2k]*product_{i=0}^{k-1}{(B+1/2+i)/(B+C+1/2+i)}}) + using uint=unsigned int; + const double B=(*pars[1])+(n+1),C=(*pars[2])+1; + const uint N=getNPar(); + double sum=1; + double prod=1; + double a=B; + double b=B+C; + for(uint i=4;i<N;i+=2){ + prod=prod*a/b; + sum+=(*pars[i])*prod; + a++; + b++; + } + double lgammaC=lgamma(C); + double ret=exp(lgamma(B)+lgammaC-lgamma(B+C))*sum; + sum=0; + prod=1; + a=B+0.5; + b=a+C; + for(uint i=3;i<N;i+=2){ + sum+=(*pars[i])*prod; + prod=prod*a/b; + a++; + b++; + } + ret+=exp(lgamma(B+0.5)+lgammaC-lgamma(B+C+0.5))*sum; + ret*=(*pars[0]); + return ret; +} +} diff --git a/pdfparams/PolySqrtPdfParam/yaml/parameters.yaml b/pdfparams/PolySqrtPdfParam/yaml/parameters.yaml new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/reactions/APPLgrid/include/ReactionAPPLgrid.h b/reactions/APPLgrid/include/ReactionAPPLgrid.h index 1d53507c0cd965a80f59961c15235b4c30c61184..6645d3aeb6982839733b6af65801b15a0cf9b7fb 100644 --- a/reactions/APPLgrid/include/ReactionAPPLgrid.h +++ b/reactions/APPLgrid/include/ReactionAPPLgrid.h @@ -4,6 +4,7 @@ #include "ReactionTheory.h" #include "appl_grid/appl_grid.h" #include <memory> +#include "BaseEvolution.h" /** @class' ReactionAPPLgrid @@ -15,33 +16,28 @@ @version 0.1 @date 2017-03-28 */ - +struct DatasetData{ + std::vector<std::shared_ptr<appl::grid> >grids; + int order; + double muR,muF; // !> renormalisation and factorisation scales + bool flagNorm; // !> if true, multiply by bin width + bool flagUseReference; // !> if true, prediction will be calculated from reference histogram (for tests and grids validation) + std::vector<TH1D*> references; + std::vector<double> eScale; // !> CMS energy + xfitter::BaseEvolution*evolutions[2]; +}; class ReactionAPPLgrid : public ReactionTheory { public: - ReactionAPPLgrid(){}; - -// ~ReactionAPPLgrid(){}; -// ~ReactionAPPLgrid(const ReactionAPPLgrid &){}; -// ReactionAPPLgrid & operator =(const ReactionAAPPLgrid &r){return *(new ReactionAPPLgrid(r));}; - - public: + ReactionAPPLgrid(); + ~ReactionAPPLgrid(); virtual string getReactionName() const { return "APPLgrid" ;}; int initAtStart(const string &); virtual void setDatasetParameters( int dataSetID, map<string,string> pars, map<string,double> parsDataset) override ; virtual int compute(int dataSetID, valarray<double> &val, map<string, valarray<double> > &err); - protected: - virtual int parseOptions(){ return 0;}; - - private: - enum class collision { pp, ppbar, pn}; - map<int, collision> _collType; - map<int, std::vector<std::shared_ptr<appl::grid> > > _grids; - map<int, int> _order; - map<int, double> _muR, _muF; // !> renormalisation and factorisation scales - map<int, bool> _flagNorm; // !> if true, multiply by bin width - map<int, bool> _flagUseReference; // !> if true, prediction will be calculated from reference histogram (for tests and grids validation) - map<int, std::vector<TH1D*> > _references; - map<int, std::vector<double> > _eScale; // !> CMS energy + protected: + virtual int parseOptions(){ return 0;}; + private: + map<int,DatasetData>dataset_data; }; diff --git a/reactions/APPLgrid/src/Makefile.am b/reactions/APPLgrid/src/Makefile.am index a9cedc2ec35a63f4dd2f0746eeac2af39a9aed33..7b776413699927f42ebc5f560121209e75a322bd 100644 --- a/reactions/APPLgrid/src/Makefile.am +++ b/reactions/APPLgrid/src/Makefile.am @@ -1,14 +1,10 @@ -# Created by AddReaction.py on 2017-03-28 - if ENABLE_APPLGRID - AM_CXXFLAGS = -I$(srcdir)/../include -I$(srcdir)/../../../include -I$(srcdir)/../../../interfaces/include -Wall -fPIC -Wno-deprecated + AM_CXXFLAGS = -I$(srcdir)/../include -I$(srcdir)/../../../include -I$(srcdir)/../../../interfaces/include -I$(srcdir)/../../../evolutions/BaseEvolution/include -Wall -fPIC -Wno-deprecated lib_LTLIBRARIES = libapplgrid_xfitter.la libapplgrid_xfitter_la_SOURCES = ReactionAPPLgrid.cc - libapplgrid_xfitter_la_LDFLAGS = $(APPLGRID_LDFLAGS) $(GFRTLIB) libapplgrid_xfitter_la_CPPFLAGS = $(APPLGRID_CPPFLAGS) $(ROOT_CFLAGS) - endif datadir = ${prefix}/yaml/reactions/APPLgrid diff --git a/reactions/APPLgrid/src/ReactionAPPLgrid.cc b/reactions/APPLgrid/src/ReactionAPPLgrid.cc index fc929facae5dc2d3f04a038e18c38ac816a37d6a..55ed4cefdfb4cb2b39de72ae79b94ed20e2bbf77 100644 --- a/reactions/APPLgrid/src/ReactionAPPLgrid.cc +++ b/reactions/APPLgrid/src/ReactionAPPLgrid.cc @@ -5,97 +5,100 @@ Created by AddReaction.py on 2017-03-28 */ -#include "ReactionAPPLgrid.h" -#include "TFile.h" - +#include"ReactionAPPLgrid.h" +#include"TFile.h" +#include"xfitter_pars.h" +#include"xfitter_steer.h" +using namespace std; +using namespace xfitter; // the class factories extern "C" ReactionAPPLgrid* create() { return new ReactionAPPLgrid(); } - +function<void(double const&x,double const&Q,double*pdfs)>active_xfxQ_functions[2]; +void xfxWrapper0(const double&x,const double&Q,double*results){active_xfxQ_functions[0](x,Q,results);} +void xfxWrapper1(const double&x,const double&Q,double*results){active_xfxQ_functions[1](x,Q,results);} +ReactionAPPLgrid::ReactionAPPLgrid(){} +ReactionAPPLgrid::~ReactionAPPLgrid(){} // Initialize at the start of the computation -int ReactionAPPLgrid::initAtStart(const string &s ) -{ - return 0; -} +int ReactionAPPLgrid::initAtStart(const string &s){return 0;} - // Initialisze for a given dataset: + // Initialize for a given dataset: void ReactionAPPLgrid::setDatasetParameters(int dataSetID, map<string,string> pars, map<string, double> parsDataset) { -// Get grid name: - if ( pars.find("GridName") != pars.end() ) { - try { - std::istringstream ss(pars["GridName"]); - //std::cout << pars["GridName"] << '\n'; - std::string token; - while(std::getline(ss, token, ',')) - { - //std::cout << token << '\n'; - std::shared_ptr<appl::grid> g(new appl::grid(token)); - g->trim(); - _grids[dataSetID].push_back(g); - TFile file(token.c_str()); - _references[dataSetID].push_back((TH1D*)file.Get("grid/reference")); - if(!_references[dataSetID].back()) - hf_errlog(17033000, "W: no reference histogram grid/reference in " + token); - else - _references[dataSetID].back()->SetDirectory(0); - } - } - catch ( const std::exception& e ) { - std::string text = "F:Failed to read APPLgrid file(s) "+pars["GridName"]; - hf_errlog_(17032802,text.c_str(),text.size()); - } - } - else { - std::string text = "F:GridName must be specified for the Reaction APPLgrid"; - hf_errlog_(17032801,text.c_str(),text.size()); - } - -// Determine order - int order = OrderMap( GetParamS("Order")); // Global order - if (pars.find("Order") != pars.end() ) { // Local order - int localOrder = OrderMap( pars["Order"] ); - order = localOrder>order ? order : localOrder; - } - - _order[dataSetID] = order; -// Determine MuR and MuF. Use default - _muR[dataSetID] = pars.find("muR") == pars.end() ? GetParam("muR") : stod(pars["muR"]); - _muF[dataSetID] = pars.find("muF") == pars.end() ? GetParam("muF") : stod(pars["muF"]); - - if (_muR[dataSetID] == 0) _muR[dataSetID] = 1.0; - if (_muF[dataSetID] == 0) _muF[dataSetID] = 1.0; - -// Determine if pp or ppbar - _collType[dataSetID] = collision::pp; - if (parsDataset.find("ppbar") != parsDataset.end() ) { - _collType[dataSetID] = collision::ppbar; + /*DEBUG + cerr<<"[DEBUG]ReactionAPPLgrid::setDatasetParameters(dataSetID="<<dataSetID<<",pars,parsDataset):\n" + <<"pars={\n"; + for(map<string,string>::const_iterator it=pars.begin();it!=pars.end();++it){ + cerr<<" \""<<it->first<<"\":\""<<it->second<<"\"\n"; } - if (parsDataset.find("pn") != parsDataset.end() ) { - _collType[dataSetID] = collision::pn; + cerr<<"}"<<endl; + cerr<<"parsDataset={\n"; + for(map<string,double>::const_iterator it=parsDataset.begin();it!=parsDataset.end();++it){ + cerr<<" \""<<it->first<<"\":"<<it->second<<'\n'; } - // check if collision settings are provided in the new format key=value - map<string,string>::iterator it = pars.find("collision"); - if (it != pars.end() ) - { - if(it->second == "pp") - _collType[dataSetID] = collision::pp; - else if(it->second == "ppbar") - _collType[dataSetID] = collision::ppbar; - else if(it->second == "pn") - _collType[dataSetID] = collision::pn; - else - hf_errlog(17102101, "F: unrecognised collision type = " + it->second); + cerr<<"}"<<endl; + DEBUG*/ + DatasetData&data=dataset_data[dataSetID]; + vector<TH1D*>&references=data.references; + BaseEvolution*(&evolutions)[2]=data.evolutions; + map<string,string>::iterator it; + // Get grid name: + it=pars.find("GridName"); + if(it!=pars.end()){ + try{ + std::istringstream ss(it->second); + std::string token; + while(std::getline(ss, token, ',')){ + std::shared_ptr<appl::grid> g(new appl::grid(token)); + g->trim(); + data.grids.push_back(g); + TFile file(token.c_str()); + references.push_back((TH1D*)file.Get("grid/reference")); + if(!references.back()) + hf_errlog(17033000, "W: no reference histogram grid/reference in " + token); + else + references.back()->SetDirectory(0); + } + } + catch ( const std::exception& e ) { + std::string text = "F:Failed to read APPLgrid file(s) "+pars["GridName"]; + hf_errlog_(17032802,text.c_str(),text.size()); + } + } + else { + std::string text = "F:GridName must be specified for the Reaction APPLgrid"; + hf_errlog_(17032801,text.c_str(),text.size()); + } +// Get pointers to evolutions + for(int i=0;i<2;++i){ + string evolutionName=""; + it=pars.find("evolution"+to_string(i+1)); + if(it!=pars.end())evolutionName=it->second;//else default evolution name will be used + //cerr<<"[DEBUG]APPLgrid: use evolution \""<<evolutionName<<'\"'<<endl; + evolutions[i]=get_evolution(evolutionName); + } +// Determine order + int order = OrderMap( GetParamS("Order")); // Global order + if (pars.find("Order") != pars.end() ) { // Local order + int localOrder = OrderMap( pars["Order"] ); + order = localOrder>order ? order : localOrder; } + data.order=order; +// Determine MuR and MuF. Use default + data.muR=pars.find("muR") == pars.end() ? GetParam("muR") : stod(pars["muR"]); + data.muF=pars.find("muF") == pars.end() ? GetParam("muF") : stod(pars["muF"]); + + if(data.muR==0)data.muR=1.0; + if(data.muF==0)data.muF=1.0; // bin width normalisation (by default no rescaling) - _flagNorm[dataSetID] = false; + data.flagNorm=false; it = pars.find("norm"); if (it != pars.end() ) { if(stoi(it->second) == 0) - _flagNorm[dataSetID] = false; + data.flagNorm=false; else if(stoi(it->second) == 1) - _flagNorm[dataSetID] = true; + data.flagNorm=true; else hf_errlog(17102102, "F: unrecognised norm = " + it->second); } @@ -104,13 +107,14 @@ void ReactionAPPLgrid::setDatasetParameters(int dataSetID, map<string,string> pa if (it != pars.end() ) { if(stoi(it->second) == 0) - _flagUseReference[dataSetID] = false; + data.flagUseReference=false; else if(stoi(it->second) == 1) { - _flagUseReference[dataSetID] = true; + data.flagUseReference=true; // check that reference histogram is available - for(std::size_t i=0; i<_references[dataSetID].size(); i++) - if(!_references[dataSetID][i]) + size_t endi=data.references.size(); + for(size_t i=0;i<endi;++i) + if(!data.references[i]) hf_errlog(17033000, "W: no reference histogram is available"); } else @@ -118,61 +122,60 @@ void ReactionAPPLgrid::setDatasetParameters(int dataSetID, map<string,string> pa } // CMS energy (by default the one used to create the grid is used) it = pars.find("energy"); - for(unsigned int i = 0; i < _grids[dataSetID].size(); i++) + unsigned int endi=data.grids.size(); + for(unsigned int i=0;i<endi;++i) { double eScale = 1.0; if (it != pars.end()) { - if(_flagUseReference[dataSetID]) + if(data.flagUseReference) hf_errlog(17110300, "W: can not apply energy scaling when using predictions from reference histogram"); else { - double eStored = _grids[dataSetID][i]->getCMSScale(); + double eStored=data.grids[i]->getCMSScale(); if(eStored < 1e-3) hf_errlog(17110301, "W: can not apply energy scaling because stored getCMSScale = 0"); else - eScale = _grids[dataSetID][i]->getCMSScale() / stof(it->second); + eScale = eStored/ stof(it->second); } } - _eScale[dataSetID].push_back(eScale); + data.eScale.push_back(eScale); } } // Main function to compute results at an iteration int ReactionAPPLgrid::compute(int dataSetID, valarray<double> &val, map<string, valarray<double> > &err) { - // - // iterate over grids - int pos = 0; - for(unsigned int g = 0; g < _grids[dataSetID].size(); g++) + DatasetData&data=dataset_data.at(dataSetID); + const int order=data.order; + const double muR=data.muR; + const double muF=data.muF; + BaseEvolution*(&evolutions)[2]=data.evolutions; + unsigned int pos = 0; + for(unsigned int g=0;g<data.grids.size();g++) { - auto grid = _grids[dataSetID][g]; + auto grid=data.grids[g]; + double eScale=data.eScale[g]; std::vector<double> gridVals(grid->Nobs()); - if(!_flagUseReference[dataSetID]) - { - // Convolute the grid: - switch (_collType[dataSetID]) - { - case collision::pp : - gridVals = grid->vconvolute( getXFX(), getAlphaS(), _order[dataSetID]-1, _muR[dataSetID], _muF[dataSetID], _eScale[dataSetID][g] ); - break; - case collision::ppbar : - gridVals = grid->vconvolute( getXFX(), getXFX("pbar"), getAlphaS(), _order[dataSetID]-1, _muR[dataSetID], _muF[dataSetID], _eScale[dataSetID][g] ); - break; - case collision::pn : - gridVals = grid->vconvolute( getXFX(), getXFX("n"), getAlphaS(), _order[dataSetID]-1, _muR[dataSetID], _muF[dataSetID], _eScale[dataSetID][g] ); - break; + if(!data.flagUseReference){ + //For some reason we do not take alphaS from evolutions? + active_xfxQ_functions[0]=evolutions[0]->xfxQArray(); + if(evolutions[0]==evolutions[1]){ + gridVals=grid->vconvolute(xfxWrapper0,getAlphaS(),order-1,muR,muF,eScale); + }else{ + active_xfxQ_functions[1]=evolutions[1]->xfxQArray(); + gridVals=grid->vconvolute(xfxWrapper0,xfxWrapper1,getAlphaS(),order-1,muR,muF,eScale); } } else { // use reference histogram for(std::size_t i=0; i<gridVals.size(); i++) - gridVals[i] = _references[dataSetID][g]->GetBinContent(i + 1); + gridVals[i]=data.references[g]->GetBinContent(i + 1); } // scale by bin width if requested - if(_flagNorm[dataSetID]) + if(data.flagNorm) for (std::size_t i=0; i<gridVals.size(); i++) gridVals[i] *= grid->deltaobs(i); @@ -181,10 +184,10 @@ int ReactionAPPLgrid::compute(int dataSetID, valarray<double> &val, map<string, std::copy_n(gridVals.begin(), gridVals.size(), &val[pos]); pos += grid->Nobs(); } - - if ( val.size() != pos ) { - hf_errlog(18072311,"F: Sum of grid sizes does not correspond to the size of the datafile "); + if(val.size()!=pos){//TODO: number of data points actually doesn't have to match grid size in some cases, so this check should be replaced by something else + std::ostringstream s; + s<<"F: ReactionAPPLgrid: Number of data points ("<<val.size()<<") in dataset (ID="<<dataSetID<<") does not match total grid size ("<<pos<<")"; + hf_errlog(18072311,s.str().c_str()); } - return 0; } diff --git a/src/Profiler.cc b/src/Profiler.cc index 42619453f5f0e02d19cfd3cf8fe3b0bfb4388159..7b86b54bad398050fbba52422489a09ed8304482 100644 --- a/src/Profiler.cc +++ b/src/Profiler.cc @@ -75,20 +75,20 @@ namespace xfitter if ( node ) { for ( auto const& term : node ) { - std::string name = term.first.as<string>(); - // check if present on the global list - if ( XFITTER_PARS::gParameters.find(name) != XFITTER_PARS::gParameters.end() ) { - profileParameter(name, term.second); - } - - else if ( name == "evolutions") { - for ( auto const& evol : term.second ) { - std::string const name = term.first.as<string>(); - // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX - profilePDF("",evol.second ); - } - } - + std::string name = term.first.as<string>(); + // check if present on the global list + if ( XFITTER_PARS::gParameters.find(name) != XFITTER_PARS::gParameters.end() ) { + profileParameter(name, term.second); + } + + else if ( name == "evolutions") { + for ( auto const& evol : term.second ) { + std::string const name = term.first.as<string>(); + // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX + profilePDF("",evol.second ); + } + } + }; } return; @@ -99,7 +99,7 @@ namespace xfitter if ( node.IsSequence() ) { size_t len = node.size(); if ( (len != 2) && (len != 3) ) { - hf_errlog( 2018082301,"S: Expected 2 or 3 parameters for a profiled variable. Check variable "+name); + hf_errlog( 2018082301,"S: Expected 2 or 3 parameters for a profiled variable. Check variable "+name); } // Evaluate varied predictions: @@ -108,18 +108,18 @@ namespace xfitter double save = *ppar; for ( size_t i=0; i<len; i++) { - double val = node[i].as<double>(); - *ppar = val; - preds.push_back( evaluatePredictions() ); + double val = node[i].as<double>(); + *ppar = val; + preds.push_back( evaluatePredictions() ); } *ppar = save; // Add to systematics list: if ( len == 2) { - addSystematics(name+":T",(preds[1]-preds[0])/preds[0]); + addSystematics(name+":T",(preds[1]-preds[0])/preds[0]); } else { - addSystematics(name+":T",(preds[1]-preds[0])/preds[0], (preds[2]-preds[0])/preds[0]); + addSystematics(name+":T",(preds[1]-preds[0])/preds[0], (preds[2]-preds[0])/preds[0]); } } } @@ -127,150 +127,146 @@ namespace xfitter void Profiler::profilePDF( std::string const& evolName, YAML::Node const& node) { // get evolution - + auto *evol = get_evolution(evolName); // get corresponding yaml node: - std::string eName = evol->getName(); - YAML::Node *gNode = & XFITTER_PARS::gParametersY[eName]; + YAML::Node gNode=XFITTER_PARS::getEvolutionNode(evolName); - //std::cout << "DEBUG: " << eName << "\n" << (*gNode) << std::endl; + //std::cout << "DEBUG: " << evolName << "\n" << (*gNode) << std::endl; //std::cout << "DEBUG3: " << node << std::endl; YAML::Node const sets = node["sets"]; YAML::Node const members = node["members"]; + //Sanity checks if ( !sets || !members ) { - hf_errlog(2018082401,"S: Profiler: missing set or member parameters for evolution "+eName); // XXXXXXXXXXXXXXXX + hf_errlog(2018082401,"S: Profiler: missing set or member parameters for evolution "+evolName); // XXXXXXXXXXXXXXXX + } + if(!sets.IsSequence()||!members.IsSequence()){ + hf_errlog(2018082402,"S: Profiler: sets and members must be sequence"); // XXXXXXXXXXXXXXXX + } + if(sets.size()!=members.size()){ + hf_errlog(2018082403,"S: Profiler: sets and members must be the same length"); // XXXXXXXXXXXXXXXX } - - if ( sets.IsSequence() && members.IsSequence() ) { - if ( sets.size() == members.size() ) { - size_t endi = sets.size(); - for (size_t i=0; i< endi; i++) { - std::string pName = sets[i].as<string>(); - - if ( members[i].IsSequence() ) { - int msize = members[i].size(); - - if ( msize != 3) { - hf_errlog(2018082404,"S: Profiler: sets must be sequence of length 3"); // XXXXXXXXXXXXXXXX - } - - int central = members[i][0].as<int>(); - int first = members[i][1].as<int>(); - int last = 0; - - try { - last = members[i][2].as<int>(); - } - catch (...) { - last = 0; /// auto-decodez - } - - // save original - auto oSet = (*gNode)["set"] ; - auto oMember = (*gNode)["member"] ; - - if ( ! oSet || ! oMember ) { - hf_errlog(2018082410,"S: No set or member variables for evolution : "+eName); - } - - // all predictions - std::vector< std::valarray<double> > preds; - - // Set central PDF and init - (*gNode)["set"] = pName; - (*gNode)["member"] = central; - preds.push_back(evaluatePredictions() ); - - - // now we can get set properties: is it hessia asymmetric, MC or symmetric hessian - std::string errorType = evol->getPropertyS("ErrorType"); - std::cout << errorType << std::endl; - - if ( last == 0) { - // auto determine XXXXXXXXXXXXXXXX - last = evol->getPropertyI("NumMembers")-1; - } - - if ( last > evol->getPropertyI("NumMembers") ) { - hf_errlog(2018082431,"W: Profller: too many members requested, use max instead"); - }; - - if (errorType == "hessian") { - if ( (last - first + 1) % 2 != 0 ) { - hf_errlog(2018082431,"S: Profller: hessian sets need even number of members. Check your inputs"); - } - if ( first % 2 == 0 ) { - hf_errlog(2018082432,"S: Profller: hessian error members should start from odd number. Check your inputs"); - } - } - - // loop over all - for (int imember = first; imember<=last; imember++) { - (*gNode)["member"] = imember; - preds.push_back( evaluatePredictions() ); - // for ( double th : preds[imember] ) { - //std::cout << th << std::endl; - //} - //std::cout << imember << std::endl; - } - - // Restore original - - (*gNode)["set"] = oSet ; - (*gNode)["member"] = oMember; - - - // Depending on error type, do nuisance parameters addition - if ( errorType == "symmhessian" ) { - for (int imember = first; imember<=last; imember++) { - addSystematics("PDF_nuisance_param_"+std::to_string( ++_ipdf )+":T",(preds[imember]-preds[0])/preds[0]); - } - } - else if ( errorType == "hessian") { - for (int imember = first; imember<=last; imember += 2) { - addSystematics("PDF_nuisance_param_"+std::to_string( ++_ipdf )+":T" - ,(preds[imember]-preds[0])/preds[0] - ,(preds[imember+1]-preds[0])/preds[0]); - } - } - else if ( errorType == "replicas") { - // construct average - for (size_t i=0; i<preds[0].size(); i++) { - preds[0][i] = 0; - } - for ( int i=first; i<=last; i++) { - preds[0] += preds[i]; - } - preds[0] /= (last-first+1); - - // convert replicas to deviations from average: - for ( int i=first; i<=last; i++) { - preds[i] -= preds[0]; - } - - // convert to eigenvectors, add to list of systematics - addReplicas(pName,preds); - } - else { - hf_errlog(2018082441,"S: Profiler Unsupported PDF error type : "+errorType); - } - - } - else { - hf_errlog(2018082404,"S: Profiler: sets must be sequence of length 3"); // XXXXXXXXXXXXXXXX - } - } + + size_t endi = sets.size(); + for (size_t i=0; i< endi; i++) { + std::string pName = sets[i].as<string>(); + + if ( members[i].IsSequence() ) { + int msize = members[i].size(); + + if ( msize != 3) { + hf_errlog(2018082404,"S: Profiler: sets must be sequence of length 3"); // XXXXXXXXXXXXXXXX + } + + int central = members[i][0].as<int>(); + int first = members[i][1].as<int>(); + int last = 0; + + try { + last = members[i][2].as<int>(); + } + catch (...) { + last = 0; /// auto-decodez + } + + // save original + auto oSet =gNode["set"]; + auto oMember=gNode["member"]; + + if ( ! oSet || ! oMember ) { + hf_errlog(2018082410,"S: No set or member variables for evolution : "+evolName); + } + + // all predictions + std::vector< std::valarray<double> > preds; + + // Set central PDF and init + gNode["set"] =pName; + gNode["member"]=central; + preds.push_back(evaluatePredictions() ); + + + // now we can get set properties: is it hessian asymmetric, MC or symmetric hessian + std::string errorType = evol->getPropertyS("ErrorType"); + std::cout << errorType << std::endl; + + if ( last == 0) { + // auto determine XXXXXXXXXXXXXXXX + last = evol->getPropertyI("NumMembers")-1; + } + + if ( last > evol->getPropertyI("NumMembers") ) { + hf_errlog(2018082431,"W: Profiler: too many members requested, use max instead"); + }; + + if (errorType == "hessian") { + if ( (last - first + 1) % 2 != 0 ) { + hf_errlog(2018082431,"S: Profiler: hessian sets need even number of members. Check your inputs"); + } + if ( first % 2 == 0 ) { + hf_errlog(2018082432,"S: Profiler: hessian error members should start from odd number. Check your inputs"); + } + } + + // loop over all + for (int imember = first; imember<=last; imember++) { + gNode["member"] = imember; + preds.push_back( evaluatePredictions() ); + // for ( double th : preds[imember] ) { + //std::cout << th << std::endl; + //} + //std::cout << imember << std::endl; + } + + // Restore original + + gNode["set"] =oSet; + gNode["member"]=oMember; + + + // Depending on error type, do nuisance parameters addition + if ( errorType == "symmhessian" ) { + for (int imember = first; imember<=last; imember++) { + addSystematics("PDF_nuisance_param_"+std::to_string( ++_ipdf )+":T",(preds[imember]-preds[0])/preds[0]); + } + } + else if ( errorType == "hessian") { + for (int imember = first; imember<=last; imember += 2) { + addSystematics("PDF_nuisance_param_"+std::to_string( ++_ipdf )+":T" + ,(preds[imember]-preds[0])/preds[0] + ,(preds[imember+1]-preds[0])/preds[0]); + } + } + else if ( errorType == "replicas") { + // construct average + for (size_t i=0; i<preds[0].size(); i++) { + preds[0][i] = 0; + } + for ( int i=first; i<=last; i++) { + preds[0] += preds[i]; + } + preds[0] /= (last-first+1); + + // convert replicas to deviations from average: + for ( int i=first; i<=last; i++) { + preds[i] -= preds[0]; + } + + // convert to eigenvectors, add to list of systematics + addReplicas(pName,preds); + } + else { + hf_errlog(2018082441,"S: Profiler Unsupported PDF error type : "+errorType); + } + } else { - hf_errlog(2018082403,"S: Profiler: sets and members must be the same length"); // XXXXXXXXXXXXXXXX + hf_errlog(2018082404,"S: Profiler: sets must be sequence of length 3"); // XXXXXXXXXXXXXXXX } } - else { - hf_errlog(2018082402,"S: Profiler: sets and members must be sequence"); // XXXXXXXXXXXXXXXX - } } void Profiler::addReplicas(std::string const& pdfName, std::vector< std::valarray<double> > const& uncertainties) { @@ -284,35 +280,35 @@ namespace xfitter for ( int i=0; i<ndata; i++) { for ( int j=i; j<ndata; j++) { - int id = i*ndata + j; - covar[id] = 0; - - for ( int k=1; k<=nrep; k++) { - covar[id] += uncertainties[k][i]*uncertainties[k][j] ; - } - - covar[id] /= nrep; - - if ( i != j ) { - int id2 = j*ndata + i; - covar[id2] = covar[id]; - } + int id = i*ndata + j; + covar[id] = 0; + + for ( int k=1; k<=nrep; k++) { + covar[id] += uncertainties[k][i]*uncertainties[k][j] ; + } + + covar[id] /= nrep; + + if ( i != j ) { + int id2 = j*ndata + i; + covar[id2] = covar[id]; + } } } double *beta = new double[ndata*ndata]; double alpha[ndata]; int ncorr = 0; getnuisancefromcovar_(ndata,ndata,ndata, - covar,beta,0, - ncorr,alpha,0); + covar,beta,0, + ncorr,alpha,0); // std::cout << "NCorr = " << beta[0] << " " << ncorr << std::endl; // ready to add systematic sources for ( int i=0; i<ncorr; i++) { valarray<double> unc(ndata); for ( int j=0; j<ndata; j++) { - unc[j] = beta[ndata*j + i]; - // std::cout << unc[j] << std::endl; + unc[j] = beta[ndata*j + i]; + // std::cout << unc[j] << std::endl; } addSystematics("PDF_nuisance_param_"+std::to_string(++_ipdf)+":T", unc/uncertainties[0] ); } diff --git a/src/ReactionTheory.cc b/src/ReactionTheory.cc index 53b7c8e49c97162bb1ef2b5d1a6013ad0ef88278..25ecd397deb1f5f83522fa27d7098e7154d0d0f6 100644 --- a/src/ReactionTheory.cc +++ b/src/ReactionTheory.cc @@ -17,7 +17,7 @@ using std::string; // Global variable to hold current function -std::function<void(double const& x, double const& Q, double* pdfs)> gProtonPdf; +std::function<void(double const& x, double const& Q, double* pdfs)> gProtonPdf; //The name is misleading, this doesn't actually have to be a proton void protonPDF(double const& x, double const& Q, double* pdfs) { gProtonPdf(x,Q,pdfs); @@ -52,16 +52,25 @@ ReactionTheory::operator=(const ReactionTheory &rt) } void ReactionTheory::initAtIteration() { - // do some magic - std::string evolName = getEvolution(); - gProtonPdf = XFITTER_PARS::retrieveXfxQArray(evolName+":p"); + try{ + gProtonPdf = XFITTER_PARS::retrieveXfxQArray(_evolution); + }catch(std::out_of_range&ex){ + std::cerr<<"Exception in "<<__func__<<": index \""<<_evolution<<"\" not present in gXfxQArrays map\n"; + std::cerr<<ex.what(); + hf_errlog(18091400,"F: Exception in retrieveXfxQArray, details written to stderr"); + } } const pXFXlike ReactionTheory::getXFX(const string& type) { - - // XXXXXXXXXXXXXXXXXXXXXX default evolution - std::string evName = XFITTER_PARS::getParameterS("Evolution"); - gProtonPdf = XFITTER_PARS::retrieveXfxQArray(evName+":p"); + //How is this different from atIteration? + + try{ + gProtonPdf = XFITTER_PARS::retrieveXfxQArray(_evolution); + }catch(std::out_of_range&ex){ + std::cerr<<"Exception in "<<__func__<<": index \""<<_evolution<<"\" not present in gXfxQArrays map\n"; + std::cerr<<ex.what(); + hf_errlog(18091400,"F: Exception in retrieveXfxQArray, details written to stderr"); + } // return _xfx[type]; // double dd[13]; diff --git a/src/TheorEval.cc b/src/TheorEval.cc index 5c72815ea4b1c2365e9eee408be964763ea9a1ec..04563b804f68537190b2c8f60948e7752dd53cfe 100644 --- a/src/TheorEval.cc +++ b/src/TheorEval.cc @@ -265,7 +265,38 @@ TheorEval::initTerm(int iterm, valarray<double> *val) } } - +//Temporary solution, for string parameters only, pending discussion +//Allows to provide dataset-specific reaction parameters in addition to dataset-provided reaction parameters +//using byReaction node in parameters.yaml +//When a parameters is given both in parameters.yaml and datafile, parameters.yaml has priority and a warning iss issued +void LoadParametersFromYAML(std::map<std::string,std::string>&pars,const std::string&reactionName){ + using std::string; + const char*BY_REACTION="byReaction"; + auto it=XFITTER_PARS::gParametersY.find(BY_REACTION); + if(it==XFITTER_PARS::gParametersY.end())return;//No overwrites given, nothing to do + YAML::Node&overwritesNode=it->second; + YAML::Node reactionNode=overwritesNode[reactionName]; + if(reactionNode.IsNull())return;//No overwrite for this reaction, nothing to do + if(!reactionNode.IsMap()){ + cerr<<"[ERROR] In "<<__func__<<"(pars,reactionName="<<reactionName<<"): expected reaction node to be a YAML Map:\n" + <<reactionNode<<"\n[/ERROR]"<<endl; + hf_errlog(18090301,"F: YAML error while loading reaction parameters, details written to stderr"); + } + try{ + for(YAML::const_iterator it=reactionNode.begin();it!=reactionNode.end();++it){ + string key=it->first.as<string>(); + auto pit=pars.find(key); + if(pit!=pars.end()){ + hf_errlog(18091700,"W: Reaction parameter in parameters.yaml overwrites dataset parameter"); + } + pars[key]=it->second.as<string>(); + } + }catch(YAML::TypedBadConversion<string>ex){ + cerr<<"[ERROR] In "<<__func__<<"(pars,reactionName="<<reactionName<<"): YAML failed to convert to string while parsing node:\n" + <<reactionNode<<"\n[/ERROR]"<<endl; + hf_errlog(18090301,"F: YAML error while loading reaction parameters, details written to stderr"); + } +} int TheorEval::initReactionTerm(int iterm, valarray<double> *val) { @@ -327,8 +358,8 @@ TheorEval::initReactionTerm(int iterm, valarray<double> *val) rt->resetParameters(XFITTER_PARS::gParametersY[term_source]); } + std::string evoName =XFITTER_PARS::getDefaultEvolutionName(); // Set the evolution: - std::string evoName = XFITTER_PARS::gParametersS.at("Evolution"); rt->setEvolution(evoName); //Retrieve evolution @@ -363,6 +394,7 @@ TheorEval::initReactionTerm(int iterm, valarray<double> *val) // split term_info into map<string, string> according to key1=value1:key2=value2:key=value3... map<string, string> pars = SplitTermInfo(term_info); + LoadParametersFromYAML(pars,rt->getReactionName()); // and transfer to the module rt->setDatasetParameters(_dsId*1000+iterm, pars, _dsPars); diff --git a/src/ftheor_eval.cc b/src/ftheor_eval.cc index 240f08e2bfc788e750db647e38324c84da1f7c7d..cf9a2667920a2cf4b782dd044a43951b07b02c14 100644 --- a/src/ftheor_eval.cc +++ b/src/ftheor_eval.cc @@ -311,19 +311,18 @@ void init_func_map_() { void init_at_iteration_() { for ( auto pdfdecomposition : XFITTER_PARS::gPdfDecompositions) { - pdfdecomposition.second->initAtIteration(); + pdfdecomposition.second->atIteration(); } - for ( auto evolution : XFITTER_PARS::gEvolutions) { - evolution.second->initAtIteration(); + for(auto it:XFITTER_PARS::gEvolutions) { + xfitter::BaseEvolution*evolution=it.second; + evolution->atIteration(); // register updated PDF XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX - const auto f = evolution.second->xfxQDouble(); - // std::cout << "Gluon(1) = " << f(0, 0.00001, 100) << std::endl; - const std::string evolName = evolution.second->getName() +":p"; + //Wait, do they even get updated between iterations? Is this here even necessary? --Ivan - XFITTER_PARS::registerXfxQArray(evolName,evolution.second->xfxQArray()); + XFITTER_PARS::registerXfxQArray(evolution->_name,evolution->xfxQArray()); } diff --git a/src/main.f b/src/main.f index 204163edf302c453ac0096b0391bf0c632c5deda..444ee8fd7ae5126fc121e695ec2a1ffce71c805e 100644 --- a/src/main.f +++ b/src/main.f @@ -65,12 +65,11 @@ C----------------------------------------------------- * * Read parameters: * - call parse_params() - -* ------------------------------------------------ -* Init new theory modules -* ------------------------------------------------ +*Read the list of dynamically loaded objects from Reactions.txt +*Confusingly, this is used not only for reactions, but also for +*minimizers, decompositions, parameterisations and evolutions call read_reactions() + call parse_params() !read parameters.yaml * * Init evolution @@ -225,16 +224,16 @@ C----------------------------------------------------- #include "steering.inc" double precision fmin, fedm, errdef integer npari, nparx, istat -C---------------------------------------------------------- +C---------------------------------------------------------- call MNSTAT(fmin, fedm, errdef, npari, nparx, istat) if (istat.eq.0) then call hf_errlog(16042801,'I: No minimization has run') else if (istat.eq.3) then call hf_errlog(16042802,'I: Successful run') else if (istat.eq.1) then - call hf_errlog(16042803,'E: Error matrix not accurate') + call hf_errlog(16042803,'E: Error matrix not accurate') else if (istat.eq.2) then - call hf_errlog(16042804,'E: Error matrix forced positive') + call hf_errlog(16042804,'E: Error matrix forced positive') endif C Save in output file open (51,file=Trim(OutDirName)//'/Status.out',status='unknown') diff --git a/src/read_steer.f b/src/read_steer.f index 3e14bbe7095f228a87ba721cd4a6eac28428efef..6d5861f20b736f1ca47eb0ef51578f67f0708d78 100644 --- a/src/read_steer.f +++ b/src/read_steer.f @@ -45,6 +45,7 @@ C call SetPDFStyle endif ! Itheory < 100 + call read_sumrules call read_mcerrorsnml ! MC uncertainties call read_chebnml ! chebyshev parameterisation extra pars call read_polynml @@ -115,8 +116,10 @@ C------------------------------------------------------ iDH_MOD = 0 ! no Dieter Heidt modifications to stat. errros. - PDFStyle = 'HERAPDF' + PDFStyle = 'HERAPDF' PDFType = 'proton' + uvalSum = 2D0 + dvalSum = 1D0 H1QCDFUNC= .False. C================================================= @@ -583,7 +586,22 @@ C--- call HF_stop end - +C +!> Read number of valence up and down quarks for sum rules +C------------------------------------------------------- + subroutine read_sumrules + implicit none +#include "pdfparam.inc" + namelist/sumrule_sums/uvalSum,dvalSum + open(51,file='steering.txt',status='old') + read(51,nml=sumrule_sums,ERR=1718,end=1717) + 1717 continue + close(51) + return + 1718 continue + print '(''Error reading namelist &sumrule_sums, STOP'')' + call HF_stop + end C !> Read MC errors namelist C------------------------------------------------------- diff --git a/src/sumrules.f b/src/sumrules.f index 7387449428d17f8230c4682864f2df696e66f679..65ae8521f706e23a16250cf0946926ff7c46946b 100644 --- a/src/sumrules.f +++ b/src/sumrules.f @@ -37,8 +37,8 @@ double precision tstr,tNoGlue,tPho *add for mixed CTEQHERA double precision SumRuleCTEQ, SumRuleCTEQhera - C----------------------------------------- + kflag=0 zero = 1d-10 @@ -85,7 +85,7 @@ C* -- sum rule : D - Dbar = 1 : gives ADval C* if (pardval(1).eq.0) then - pardval(1) = 1.0d0/CalcIntPdf(pardval) + pardval(1)=dvalSum/CalcIntPdf(pardval) else dv_sum = pardval(1)*CalcIntPdf(pardval) endif @@ -94,10 +94,11 @@ C********************************************************** C* -- sum rule : U - Ubar = 2 : gives AUval C* if (paruval(1).eq.0) then - paruval(1) = 2.0D0/CalcIntPdf(paruval) + paruval(1)=uvalSum/CalcIntPdf(paruval) else uv_sum = paruval(1)*CalcIntPdf(paruval)/2. endif +C* --TODO: cvalSum sumrule here? C Also integrate momenta, for momentum sum rule: tUv = paruval(1)*CalcIntXpdf(paruval) @@ -934,11 +935,11 @@ C--------------------------------------------------------------- C Counting sum-rule for uv: sumUv = SumRuleASpar(-1,asuval) - asuval(1) = 2.0D0 / sumUv + asuval(1) = uvalSum / sumUv C Counting sum-rule for dv: sumDv = SumRuleASpar(-1,asdval) - asdval(1) = 1.0D0 / sumDv + asdval(1) = dvalSum / sumDv C Momentum sum rule: sumMom = 2.D0*asubar(1)*SumRuleASpar(0,asubar) + @@ -1150,11 +1151,11 @@ C--------------------------------------------------------------- C Counting sum-rule for uv: sumUv = SumRuleCTEQ(-1,ctuval) - ctuval(1) = 2.0D0 / sumUv + ctuval(1) = uvalSum / sumUv C Counting sum-rule for dv: sumDv = SumRuleCTEQ(-1,ctdval) - ctdval(1) = 1.0D0 / sumDv + ctdval(1) = dvalSum / sumDv C Momentum sum rule: C---------------- @@ -1167,7 +1168,6 @@ C Sea: tStr = 0 ! Strange already included in Dbar endif - sumMom = 2.D0*ctubar(1)*SumRuleCTEQ(0,ctubar) + $ 2.D0*ctdbar(1)*SumRuleCTEQ(0,ctdbar) + $ ctuval(1)*SumRuleCTEQ(0,ctuval) + @@ -1228,11 +1228,11 @@ C--------------------------------------------------------------- C Counting sum-rule for uv: sumUv = SumRuleCTEQhera(-1,ctuval) - ctuval(1) = 2.0D0 / sumUv + ctuval(1) = uvalSum / sumUv C Counting sum-rule for dv: sumDv = SumRuleCTEQhera(-1,ctdval) - ctdval(1) = 1.0D0 / sumDv + ctdval(1) = dvalSum / sumDv C Momentum sum rule: C---------------- diff --git a/src/xfitter_pars.cc b/src/xfitter_pars.cc index 7e0e63ae3c7ee82a82e20681fc539a3e0aeb318f..bfcf0765d6353f5b8a95e27c8a6ab3db70577841 100644 --- a/src/xfitter_pars.cc +++ b/src/xfitter_pars.cc @@ -8,10 +8,12 @@ */ #include "xfitter_pars.h" +#include "xfitter_steer.h" #include "xfitter_cpp.h" #include "xfitter_cpp_base.h" #include <fstream> #include <string.h> +#include <cmath> #include "BaseEvolution.h" #include "BasePdfDecomposition.h" #include "BaseMinimizer.h" @@ -21,12 +23,12 @@ extern "C" { void parse_params_(); //!< Fortran callable parsing /// 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); + 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); void add_to_param_map_(map<std::string,double*>* map,double &value, int& global, char *name, int len); // Get parameters in fortran, for backward compatibility: double getparamd_(const char* name, int len); @@ -39,6 +41,7 @@ namespace XFITTER_PARS { // Global vars: xfitter::BaseMinimizer* gMinimizer(nullptr); + YAML::Node rootNode; map <string, double*> gParameters; map <string, int> gParametersI; map <string, string> gParametersS; @@ -49,9 +52,70 @@ namespace XFITTER_PARS { // Also keep list of loaded evolutions here: map<string,xfitter::BaseEvolution*> gEvolutions; - // Also keep list of loaded evolutions here: + // Also keep list of loaded decompositions here: map<string,xfitter::BasePdfDecomposition*> gPdfDecompositions; - + // Also keep list of loaded parameterisations here: + map<string,xfitter::BasePdfParam*> gParameterisations; + + using namespace xfitter; + xfitter::InitialPDFfunction getInputFunctionFromYaml(const YAML::Node&rootNode){ + YAML::Node node=rootNode["decomposition"]; + string name; + try{ + name=node.as<string>(); + }catch(YAML::TypedBadConversion<string>ex){ + ostringstream s;s<<"W: YAML exception: "<<ex.what()<<"; while trying to extract decomposition name from node: "<<node<<"; using default decomposition name"; + hf_errlog(18082930,s.str()); + name=getDefaultDecompositionName(); + } + return xfitter::get_pdfDecomposition(name)->f0(); + } + YAML::Node getEvolutionNode(const std::string&name){ + auto it=gParametersY.find("Evolutions"); + if(it==gParametersY.end()){ + hf_errlog(18082900,"F:Failed to get evolution "+name+": Evolutions node not found in parameters.yaml"); + } + YAML::Node instanceNode=it->second[name]; + if(!instanceNode.IsMap()){ + std::ostringstream s; + s<<"F:Failed to get evolution \""<<name<<"\": "; + if(!instanceNode)s<<"no subnode with this name under the node Evolutions"; + else s<<"corresponding subnode is not of type Map"; + hf_errlog(18082901,s.str().c_str()); + } + return instanceNode; + } + YAML::Node getDecompositionNode(const std::string&name){ + auto it=gParametersY.find("Decompositions"); + if(it==gParametersY.end()){ + hf_errlog(18082902,"F:Failed to get decomposition "+name+": \"Decompositions\" node not found"); + } + YAML::Node instanceNode=it->second[name]; + if(!instanceNode.IsMap()){ + std::ostringstream s; + s<<"F:Failed to get decomposition \""<<name<<"\": "; + if(!instanceNode)s<<"no subnode with this name under the node Decompositions"; + else s<<"corresponding subnode is not of type Map"; + hf_errlog(18082903,s.str().c_str()); + } + return instanceNode; + } + YAML::Node getParameterisationNode(const std::string&name){ + auto it=gParametersY.find("Parameterisations"); + if(it==gParametersY.end()){ + hf_errlog(18082904,"F:Failed to get parameterisation "+name+": \"Parameterisations\" node not found"); + } + YAML::Node instanceNode=it->second[name]; + if(!instanceNode.IsMap()){ + std::ostringstream s; + s<<"F:Failed to get parameterisation \""<<name<<"\": "; + if(!instanceNode)s<<"no subnode with this name under the node Parameterisations"; + else s<<"corresponding subnode is not of type Map"; + hf_errlog(18082905,s.str().c_str()); + } + return instanceNode; + } + std::string getParameterS(std::string name) { auto search = gParametersS.find(name); if ( search != gParametersS.end() ) { @@ -62,7 +126,17 @@ namespace XFITTER_PARS { return ""; // not found } } - + + string getDefaultEvolutionName(){ + auto it=gParametersS.find("DefaultEvolution"); + if(it!=gParametersS.end()) return it->second; + return "default"; + } + string getDefaultDecompositionName(){ + auto it=gParametersS.find("DefaultDecomposition"); + if(it!=gParametersS.end())return it->second; + return "default"; + } // Helper function bool is_file_exist(const char *fileName) { @@ -74,11 +148,13 @@ namespace XFITTER_PARS { { try { if ( ! std::ifstream(name).good()) { - string text = "F: Problems opening parameters file " + name; - hf_errlog_(18032001,text.c_str(), text.size()); + string text = "F: Problems opening parameters file " + name; + hf_errlog_(18032001,text.c_str(), text.size()); } YAML::Node node = YAML::LoadFile(name); parse_node(node, gParameters, gParametersI, gParametersS, gParametersV, gParametersY); + //HACKY way to get rootNode, pending include rewrite + if(rootNode.IsNull())rootNode=node; } catch (const std::exception& e) { std::cout << e.what() << std::endl; @@ -93,11 +169,11 @@ namespace XFITTER_PARS { // Parse @param node and return maps void parse_node(const YAML::Node& node, - std::map<string,double*>& dMap, - std::map<string,int>& iMap, - std::map<string,string>& sMap, - std::map<string,vector<double> >& vMap, - std::map<string,YAML::Node> & yMap ) + std::map<string,double*>& dMap, + std::map<string,int>& iMap, + std::map<string,string>& sMap, + std::map<string,vector<double> >& vMap, + std::map<string,YAML::Node> & yMap ) { for ( YAML::const_iterator it = node.begin(); it != node.end(); ++it) { YAML::Node key = it->first; @@ -108,95 +184,95 @@ namespace XFITTER_PARS { // Check if asked to read another file: if ( p_name == "include" ) { - auto fileName = value.as<string>(); - if (is_file_exist(fileName.c_str())) { - parse_file( fileName ); - } - else { - // Now try default location: - if (is_file_exist((PREFIX +string("/")+fileName).c_str())) { - parse_file( PREFIX +string("/")+fileName ); - } - else { - string msg = "F: Include Yaml parameters file "+fileName+" not found"; - hf_errlog_(17041601,msg.c_str(), msg.size()); - } - } + auto fileName = value.as<string>(); + if (is_file_exist(fileName.c_str())) { + parse_file( fileName ); + } + else { + // Now try default location: + if (is_file_exist((PREFIX +string("/")+fileName).c_str())) { + parse_file( PREFIX +string("/")+fileName ); + } + else { + string msg = "F: Include Yaml parameters file "+fileName+" not found"; + hf_errlog_(17041601,msg.c_str(), msg.size()); + } + } } if (value.IsScalar()) { - // Alright, store directly - // Try to read as int, float, string: - try { - int i = value.as<int>(); - iMap[p_name] = i; - continue; - } - catch (const std::exception& e) { - } - - try { - double f = value.as<double>(); - dMap[p_name] = new double(f); - continue; - } - catch (const std::exception& e) { - } - - try { - std::string s = value.as<std::string>(); - sMap[p_name] = s; - continue; - } - catch (const std::exception& e) { - } + // Alright, store directly + // Try to read as int, float, string: + try { + int i = value.as<int>(); + iMap[p_name] = i; + continue; + } + catch (const std::exception& e) { + } + + try { + double f = value.as<double>(); + dMap[p_name] = new double(f); + continue; + } + catch (const std::exception& e) { + } + + try { + std::string s = value.as<std::string>(); + sMap[p_name] = s; + continue; + } + catch (const std::exception& e) { + } } else { // Potentially this may go to minuit, if step is not zero. - - if (value.IsMap()) { + + if (value.IsMap()) { - // Check if this is a minimisation block, true if step is present - if (value["step"] || value["value"]) { - // Defaults - double val = 0; - double step = 0; - double minv = 0; - double maxv = 0; - double priorVal = 0; - double priorUnc = 0; - int add = true; - - if (value["value"]) { - val = value["value"].as<double>(); - } - else { - string text = "F: missing value field for parameter " + p_name; - hf_errlog_(17032401,text.c_str(),text.size()); - } - if (value["step"]) step = value["step"].as<double>(); - if (value["prior"]) priorVal = value["prior"].as<double>(); - if (value["priorUnc"]) priorUnc = value["priorUnc"].as<double>(); - if (value["min"]) minv = value["min"].as<double>(); - if (value["max"]) maxv = value["max"].as<double>(); - // Goes to fortran - addexternalparam_(p_name.c_str(), val, step, minv, maxv, - priorVal, priorUnc, add, &dMap, p_name.size()); - } - else { - // no step or value, store as it is as a yaml node: - yMap[p_name] = value; - } - } - - else if (value.IsSequence() ) { - size_t len = value.size(); - vector<double> v(len); - for (size_t i=0; i<len; i++) { - v[i] = value[i].as<double>(); - } - vMap[p_name] = v; - } - } + // Check if this is a minimisation block, true if step is present + if (value["step"] || value["value"]) { + // Defaults + double val = 0; + double step = 0; + double minv = 0; + double maxv = 0; + double priorVal = 0; + double priorUnc = 0; + int add = true; + + if (value["value"]) { + val = value["value"].as<double>(); + } + else { + string text = "F: missing value field for parameter " + p_name; + hf_errlog_(17032401,text.c_str(),text.size()); + } + if (value["step"]) step = value["step"].as<double>(); + if (value["prior"]) priorVal = value["prior"].as<double>(); + if (value["priorUnc"]) priorUnc = value["priorUnc"].as<double>(); + if (value["min"]) minv = value["min"].as<double>(); + if (value["max"]) maxv = value["max"].as<double>(); + // Goes to fortran + addexternalparam_(p_name.c_str(), val, step, minv, maxv, + priorVal, priorUnc, add, &dMap, p_name.size()); + } + else { + // no step or value, store as it is as a yaml node: + yMap[p_name] = value; + } + } + + else if (value.IsSequence() ) { + size_t len = value.size(); + vector<double> v(len); + for (size_t i=0; i<len; i++) { + v[i] = value[i].as<double>(); + } + vMap[p_name] = v; + } + } } } @@ -264,10 +340,116 @@ namespace XFITTER_PARS { const std::function<void(double const& x, double const& Q, double* pdfs)> retrieveXfxQArray(const std::string& name) { return gXfxQArrays.at(name); } + void createParameters(){ + using namespace std; + try{ + YAML::Node parsNode=gParametersY["Parameters"]; + if(!parsNode.IsMap()){ + hf_errlog(18091710,"F: Failed to create parameters: bad \"Parameters\" YAML node"); + } + xfitter::BaseMinimizer*minimizer=xfitter::get_minimizer(); + for(YAML::const_iterator it=parsNode.begin();it!=parsNode.end();++it){ + string parameterName=it->first.as<string>(); + double value=nan(""); + double step=nan(""); + double min=nan(""); + double max=nan(""); + double pr_mean=nan(""); + double pr_sigma=nan(""); + //TODO: bounds and priors + YAML::Node pNode=it->second; + switch(pNode.Type()){ + case YAML::NodeType::Scalar:{//Should be a special string DEPENDENT + string key=pNode.as<string>(); + if(key=="DEPENDENT"){ + value=1; + step=0;//This means that this parameter will be calculated using sum rules + } + else hf_errlog(18091712,"F: Bad parameter definition"); + break;} + case YAML::NodeType::Sequence:{// interpret sequence as [value,step,min,max] + int size=pNode.size(); + if(size>0)value=pNode[0].as<double>(); + else break; + if(size>1)step=pNode[1].as<double>(); + else break; + if(size>2)min=pNode[2].as<double>(); + else break; + if(size>3)max=pNode[3].as<double>(); + else break; + if(size>4)pr_mean=pNode[4].as<double>(); + else break; + if(size>5)pr_sigma=pNode[5].as<double>(); + else break; + if(size>6){ + cerr<<"[ERROR] Too many arguments("<<size<<") in definition of parameter "<<parameterName<<"; expected at most 6"<<endl; + hf_errlog(18091716,"F: Too many arguments in parameter definition, see stderr"); + } + break;} + case YAML::NodeType::Map:{ + for(YAML::const_iterator it=pNode.begin();it!=pNode.end();++it){ + string key=it->first.as<string>(); + if (key=="value")value=it->second.as<double>(); + else if(key=="step")step=it->second.as<double>(); + else if(key=="min")min=it->second.as<double>(); + else if(key=="max")max=it->second.as<double>(); + else if(key=="pr_mean")pr_mean=it->second.as<double>(); + else if(key=="pr_sigma")pr_sigma=it->second.as<double>(); + else{ + cerr<<"[ERROR] Unknown key "<<key<<" in definition of parameter "<<parameterName<<endl; + hf_errlog(18091711,"F: Unknown key in parameter definition, see stderr"); + } + } + break;} + default: + hf_errlog(18091710,"F: Failed to parse parameters definition node, expected Sequence or Map"); + } + if(std::isnan(value)){ + value=1; + hf_errlog(18091713,"W: Initial value not given for a parameter, assuming 1.0"); + } + if(std::isnan(step)||step<0){ + step=fabs(value)*0.01; + hf_errlog(18091714,"W: Step not given for a parameter, assuming 1%"); + } + if(std::isnan(min)^std::isnan(max)){ + cerr<<"[ERROR] Only one of two bounds ("<<(std::isnan(min)?"upper":"lower")<<") is given for parameter "<<parameterName<<endl; + hf_errlog(18091715,"F: Only one of two bounds is given for one of parameters, see stderr"); + } + if(std::isnan(pr_mean)^std::isnan(pr_sigma)){ + if(std::isnan(pr_sigma)){ + pr_sigma=1; + hf_errlog(18091716,"W: Prior mean value, but not sigma given for a parameter, assuming sigma=1.0"); + }else{ + cerr<<"[ERROR] Prior sigma value, but not prior mean fiven for parameter "<<parameterName<<endl; + hf_errlog(18091717,"F: Prior sigma value, but not mean given for a parameter, see stderr"); + } + } + //Note that if we allocate memory here, it will never be freed + double*bounds=nullptr; + double*priors=nullptr; + if(!std::isnan(min)){ + bounds=new double[2]; + bounds[0]=min; + bounds[1]=max; + } + if(!std::isnan(pr_mean)){ + priors=new double[2]; + priors[0]=pr_mean; + priors[1]=pr_sigma; + } + minimizer->addParameter(value,parameterName,step,bounds,priors); + } + }catch(YAML::Exception&ex){ + cerr<<"[ERROR] YAML exception:\n"<<ex.what()<<"\n[/ERROR]"<<endl; + hf_errlog(18091713,"F: YAML exception while creating parameters, details written to stderr"); + } + } } void parse_params_(){ XFITTER_PARS::parse_file("parameters.yaml"); + XFITTER_PARS::createParameters(); XFITTER_PARS::ParsToFortran(); } diff --git a/src/xfitter_steer.cc b/src/xfitter_steer.cc index 2bb75cca468eea062bad9e8399c453bb6929eed5..22892010b601e6c03eeb5498067f05d6b22c1f67 100644 --- a/src/xfitter_steer.cc +++ b/src/xfitter_steer.cc @@ -9,105 +9,127 @@ #include <iostream> #include <yaml-cpp/yaml.h> #include <Profiler.h> +using std::string; extern std::map<string,string> gReactionLibs; -namespace xfitter { - - BaseEvolution* get_evolution(std::string name) { - if (name == "") { - // get the name from the map - name = XFITTER_PARS::getParameterS("Evolution"); - } - - // Check if already present - if (XFITTER_PARS::gEvolutions.count(name) == 1) { - return XFITTER_PARS::gEvolutions[name]; //already loaded - } - - - // Load corresponding shared library: - string libname = gReactionLibs[name]; - if ( libname == "") { - hf_errlog(18071302,"F: Shared library for evolution "+name+" not found"); +void*createDynamicObject(const string&classname,const string&instanceName){ + //instantiate an object from a shared library + //Used to create evolution, decomposition, parameterisation, could be used to create minimizer + string libpath; + try{ + libpath=PREFIX+string("/lib/")+gReactionLibs.at(classname); + }catch(const std::out_of_range&ex){ + std::cerr<<"[ERROR] out_of_range in function "<<__func__<<":\n"<<ex.what()<<"\n[/ERROR]\n"; + std::ostringstream s; + if(gReactionLibs.count(classname)==0){ + s<<"F: Unknown dynamically loaded class \""<<classname<<"\""; + hf_errlog(18091901,s.str().c_str()); } + s<<"F: Unknown out_of_range exception in "<<__func__; + hf_errlog(18091902,s.str().c_str()); + } + void*shared_library=dlopen(libpath.c_str(),RTLD_NOW); + //By the way, do we ever call dlclose? I don't think so... Maybe we should call it eventually. --Ivan Novikov + if(shared_library==NULL){ + std::cerr<<"[ERROR] dlopen() error while trying to open shared library for class \""<<classname<<"\":\n"<<dlerror()<<"\n[/ERROR]"<<std::endl; + hf_errlog(18091900,"F:dlopen() error, see stderr"); + } + //reset errors + dlerror(); + void*create=dlsym(shared_library,"create"); + if(create==NULL){ + std::cerr<<"[ERROR] dlsym() failed to find \"create\" function for class \""<<classname<<"\":\n"<<dlerror()<<"\n[/ERROR]"<<std::endl; + hf_errlog(18091902,"F:dlsym() error, see stderr"); + } + return((void*(*)(const char*))create)(instanceName.c_str()); +} - // load the library: - void *evolution_handler = dlopen((PREFIX+string("/lib/")+libname).c_str(), RTLD_NOW); - if (evolution_handler == NULL) { - std::cout << dlerror() << std::endl; - hf_errlog(18071303,"F: Evolution shared library ./lib/" + libname + " not present for evolution" + name + ". Check Reactions.txt file"); +namespace xfitter { + BaseEvolution*get_evolution(string name){ + if(name=="")name=XFITTER_PARS::getDefaultEvolutionName(); + // Check if already present + if(XFITTER_PARS::gEvolutions.count(name)==1){ + return XFITTER_PARS::gEvolutions.at(name); } - - // reset errors - dlerror(); - - create_evolution *dispatch_ev = (create_evolution*) dlsym(evolution_handler, "create"); - BaseEvolution *evolution = dispatch_ev(); - - // Now we attach corresponding PDFdecomposition. First try specific, next: global - std::string pdfDecomp = XFITTER_PARS::getParameterS("PDFDecomposition"); - - // XXXXXXXXXXXXXXXXXXXXXXXX - if ( XFITTER_PARS::gParametersY.count(name) > 0 ) { - auto evolNode = XFITTER_PARS::gParametersY[name]; - if ( evolNode["PDFDecomposition"] ) { - pdfDecomp = evolNode["PDFDecomposition"].as<std::string>(); - std::cout << " here here \n"; - } - else { - std::cout << " ho here \n"; - } + // Else create a new instance of evolution + YAML::Node instanceNode=XFITTER_PARS::getEvolutionNode(name); + YAML::Node classnameNode=instanceNode["class"]; + if(!classnameNode.IsScalar()){ + std::ostringstream s; + s<<"F:Failed to get evolution \""<<name<<"\": evolution must have a node \"class\" with the class name as a string"; + hf_errlog(18082902,s.str().c_str()); } - - std::cout << "PDF decomp=" << pdfDecomp << "\n"; - - // Get corresponding PDF decomposition: - BasePdfDecomposition* pdfD = get_pdfDecomposition(pdfDecomp); - evolution->SetPdfDecomposition( pdfD->f0() ); - - // Init it: - evolution->initAtStart(); - - // Store on the map + string classname=classnameNode.as<string>(); + BaseEvolution*evolution=(BaseEvolution*)createDynamicObject(classname,name); + //Note that unlike in the pervious version of this function, we do not set decompositions for evolutions + //Evolution objects are expected to get their decomposition themselves based on YAML parameters, during atStart + evolution->atStart(); + // Store the newly created evolution on the global map XFITTER_PARS::gEvolutions[name] = evolution; return evolution; } - - BasePdfDecomposition* get_pdfDecomposition(std::string name) { - if (name == "") { - // get the name from the map - name = XFITTER_PARS::getParameterS("PDFDecomposition"); - } - // Check if already present - if (XFITTER_PARS::gPdfDecompositions.count(name) == 1) { - return XFITTER_PARS::gPdfDecompositions[name]; //already loaded - } - - // Load corresponding shared library: - string libname = gReactionLibs[name]; - if ( libname == "") { - hf_errlog(18072302,"F: Shared library for pdf decomposition "+name+" not found"); + BasePdfDecomposition*get_pdfDecomposition(string name){ + try{ + if(name=="")name=XFITTER_PARS::getDefaultDecompositionName(); + auto it=XFITTER_PARS::gPdfDecompositions.find(name); + if(it!=XFITTER_PARS::gPdfDecompositions.end())return it->second; + string classname=XFITTER_PARS::getDecompositionNode(name)["class"].as<string>(); + BasePdfDecomposition*ret=(BasePdfDecomposition*)createDynamicObject(classname,name); + ret->atStart(); + XFITTER_PARS::gPdfDecompositions[name]=ret; + return ret; + }catch(YAML::InvalidNode&ex){ + const int errcode=18092401; + const char*errmsg="F: YAML::InvalidNode exception while creating decomposition, details written to stderr"; + using namespace std; + cerr<<"[ERROR]"<<__func__<<'('<<name<<')'<<endl; + YAML::Node node=XFITTER_PARS::getDecompositionNode(name); + if(!node.IsMap()){ + cerr<<"Invalid node Decompositions/"<<name<<"\nnode is not a map\n[/ERROR]"<<endl; + hf_errlog(errcode,errmsg); + } + YAML::Node node_class=node["class"]; + if(!node_class.IsScalar()){ + if(node_class.IsNull())cerr<<"Missing node Decompositions/"<<name<<"/class"; + else cerr<<"Invalid node Decompositions/"<<name<<"/class\nnode is not a scalar"; + cerr<<"\n[/ERROR]"<<endl; + hf_errlog(errcode,errmsg); + } + cerr<<"Unexpected YAML exception\nNode:\n"<<node<<"\n[/ERROR]"<<endl; + hf_errlog(errcode,errmsg); } - - // load the library: - void *pdfDecomposition_handler = dlopen((PREFIX+string("/lib/")+libname).c_str(), RTLD_NOW); - if (pdfDecomposition_handler == NULL) { - std::cout << dlerror() << std::endl; - hf_errlog(18072303,"F: PdfDecomposition shared library ./lib/" + libname + " not present for pdfDecomposition" + name + ". Check Reactions.txt file"); + } + BasePdfParam*getParameterisation(const string&name){ + try{ + auto it=XFITTER_PARS::gParameterisations.find(name); + if(it!=XFITTER_PARS::gParameterisations.end())return it->second; + //Else create a new instance + string classname=XFITTER_PARS::getParameterisationNode(name)["class"].as<string>(); + BasePdfParam*ret=(BasePdfParam*)createDynamicObject(classname,name); + ret->atStart(); + XFITTER_PARS::gParameterisations[name]=ret; + return ret; + }catch(YAML::InvalidNode&ex){ + const int errcode=18092400; + const char*errmsg="F: YAML::InvalidNode exception while creating parameterisation, details written to stderr"; + using namespace std; + cerr<<"[ERROR]"<<__func__<<'('<<name<<')'<<endl; + YAML::Node node=XFITTER_PARS::getParameterisationNode(name); + if(!node.IsMap()){ + cerr<<"Invalid node Parameterisations/"<<name<<"\nnode is not a map\n[/ERROR]"<<endl; + hf_errlog(errcode,errmsg); + } + YAML::Node node_class=node["class"]; + if(!node_class.IsScalar()){ + if(node_class.IsNull())cerr<<"Missing node Parameterisations/"<<name<<"/class"; + else cerr<<"Invalid node Parameterisations/"<<name<<"/class\nnode is not a scalar"; + cerr<<"\n[/ERROR]"<<endl; + hf_errlog(errcode,errmsg); + } + cerr<<"Unexpected YAML exception\nNode:\n"<<node<<"\n[/ERROR]"<<endl; + hf_errlog(errcode,errmsg); } - - // reset errors - dlerror(); - - create_pdfDecomposition *dispatch_decomp = (create_pdfDecomposition*) dlsym(pdfDecomposition_handler, "create"); - BasePdfDecomposition *pdfDecomp = dispatch_decomp(); - pdfDecomp->initAtStart(""); - - // store on the map - XFITTER_PARS::gPdfDecompositions[name] = pdfDecomp; - - return pdfDecomp; } @@ -137,7 +159,7 @@ namespace xfitter { create_minimizer *dispatch_minimizer = (create_minimizer*) dlsym(lib_handler, "create"); BaseMinimizer *minimizer = dispatch_minimizer(); - minimizer->initAtStart(); + minimizer->atStart(); // store on the map XFITTER_PARS::gMinimizer = minimizer; @@ -157,11 +179,12 @@ extern "C" { } void init_evolution_() { - auto evol = xfitter::get_evolution(); + //TODO: reimplement for new interface with multiple evolutions + //auto evol = xfitter::get_evolution(); } void init_minimizer_() { - /// initAtStart is called inside + /// atStart is called inside auto mini = xfitter::get_minimizer(); } diff --git a/steering.txt b/steering.txt index 2aa8a012e66102045b93ee157a4aeccc46339252..580cc1d21ac3e92db6f860c0fb26dc3e25693190 100644 --- a/steering.txt +++ b/steering.txt @@ -194,7 +194,7 @@ CHI2SettingsName = 'StatScale', 'UncorSysScale', 'CorSysScale', 'UncorChi2Type', 'CorChi2Type' Chi2Settings = 'Poisson' , 'Linear', 'Linear' , 'Diagonal' , 'Hessian' - ! Chi2ExtraParam = 'PoissonCorr' + Chi2ExtraParam = 'PoissonCorr' ! Flag to define if native APPLgrid CKM values should be kept. LUseAPPLgridCKM = True diff --git a/tools/AddPdfParam.py b/tools/AddPdfParam.py index 227e9684fcf7e544ee93289f39fc17933f79b55a..a6f4239b34811d7212b53e79172c388fd056ce39 100644 --- a/tools/AddPdfParam.py +++ b/tools/AddPdfParam.py @@ -45,16 +45,16 @@ with open(hFile,"w+") as f: class {name:s}PdfParam:public BasePdfParam{{ public: - {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; + {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,date=datetime.date.today().isoformat()) ) diff --git a/tools/draw/src/xfitter-draw.cc b/tools/draw/src/xfitter-draw.cc index c98b6200488601ed326a222b41100bae7f5ca320..ea5d5cef9588c5e9ef175c85a9fbe7a997096d7a 100644 --- a/tools/draw/src/xfitter-draw.cc +++ b/tools/draw/src/xfitter-draw.cc @@ -35,13 +35,13 @@ int main(int argc, char **argv) for (vector<string>::iterator it1 = opts.labels.begin(); it1 != opts.labels.end(); it1++) for (vector<string>::iterator it2 = it1+1; it2 != opts.labels.end(); it2++) if (*it1 == *it2) - { - cout << endl; - cout << "Error: label (or directory) " << *it1 << " can appear only once in labels list" << endl; - cout << "Specify different labels" << endl; - cout << endl; - exit(-1); - } + { + cout << endl; + cout << "Error: label (or directory) " << *it1 << " can appear only once in labels list" << endl; + cout << "Specify different labels" << endl; + cout << endl; + exit(-1); + } //Associate colors and styles to labels for (vector<string>::iterator itl = opts.labels.begin(); itl != opts.labels.end(); itl++) @@ -67,19 +67,19 @@ int main(int argc, char **argv) //loop on Q2 bins vector <float> ql = q2list(); for (vector<float>::iterator qit = ql.begin(); qit != ql.end(); qit++) - { - //loop on pdf types - for (vector <pdftype>::iterator pit = pdfs.begin(); pit != pdfs.end(); pit++) - { - vector <TCanvas*> pdfcnv = PdfsPainter(*qit, *pit); - pdfscanvaslist.push_back(pdfcnv[0]); - if (pdfcnv.size() > 1) - pdfscanvasratiolist.push_back(pdfcnv[1]); - } - - if (!opts.q2all) - break; - } + { + //loop on pdf types + for (vector <pdftype>::iterator pit = pdfs.begin(); pit != pdfs.end(); pit++) + { + vector <TCanvas*> pdfcnv = PdfsPainter(*qit, *pit); + pdfscanvaslist.push_back(pdfcnv[0]); + if (pdfcnv.size() > 1) + pdfscanvasratiolist.push_back(pdfcnv[1]); + } + + if (!opts.q2all) + break; + } } //-------------------------------------------------- @@ -91,14 +91,14 @@ int main(int argc, char **argv) //loop on datasets vector <int> dl = datalist(); for (vector<int>::iterator dit = dl.begin(); dit != dl.end(); dit++) - { - //extract dataset index and subplot index - int dataindex = (int)(*dit) / 100; - int subplotindex = *dit - dataindex * 100; - TCanvas *dataplot = DataPainter(dataindex, subplotindex); - if (dataplot != 0) - datapullscanvaslist.push_back(dataplot); - } + { + //extract dataset index and subplot index + int dataindex = (int)(*dit) / 100; + int subplotindex = *dit - dataindex * 100; + TCanvas *dataplot = DataPainter(dataindex, subplotindex); + if (dataplot != 0) + datapullscanvaslist.push_back(dataplot); + } } //-------------------------------------------------- @@ -141,11 +141,11 @@ int main(int argc, char **argv) TCanvas * pagecnv = new TCanvas(numb, "", 0, 0, opts.resolution * opts.plotsperpage, opts.resolution * opts.plotsperpage); pagecnv->Divide(opts.plotsperpage, opts.plotsperpage); for (int i = 1; i <= opts.plotsperpage*opts.plotsperpage && it != pdfscanvaslist.end(); i++) - { - pagecnv->cd(i); - (*it)->DrawClonePad(); - it++; - } + { + pagecnv->cd(i); + (*it)->DrawClonePad(); + it++; + } pgn++; sprintf(pgnum, "%d", pgn); @@ -158,11 +158,11 @@ int main(int argc, char **argv) TCanvas * pagecnv = new TCanvas(numb, "", 0, 0, opts.resolution *opts.plotsperpage, opts.resolution * opts.plotsperpage); pagecnv->Divide(opts.plotsperpage, opts.plotsperpage); for (int i = 1; i <= opts.plotsperpage*opts.plotsperpage && it != pdfscanvasratiolist.end(); i++) - { - pagecnv->cd(i); - (*it)->DrawClonePad(); - it++; - } + { + pagecnv->cd(i); + (*it)->DrawClonePad(); + it++; + } pgn++; sprintf(pgnum, "%d", pgn); pagecnv->Print((opts.outdir + "plots_" + pgnum + ".eps").c_str()); @@ -179,29 +179,29 @@ int main(int argc, char **argv) sprintf(numb, "data_%d", it - datapullscanvaslist.begin()); TCanvas * pagecnv; if (opts.twopanels || opts.threepanels) - { - pagecnv = new TCanvas(numb, "", 0, 0, opts.resolution * 2, opts.resolution * 2); - pagecnv->Divide(1, 2); - for (int i = 1; i <= 2; i++) - if (it != datapullscanvaslist.end()) - { - pagecnv->cd(i); - (*it)->DrawClonePad(); - it++; - } - } + { + pagecnv = new TCanvas(numb, "", 0, 0, opts.resolution * 2, opts.resolution * 2); + pagecnv->Divide(1, 2); + for (int i = 1; i <= 2; i++) + if (it != datapullscanvaslist.end()) + { + pagecnv->cd(i); + (*it)->DrawClonePad(); + it++; + } + } else - { - pagecnv = new TCanvas(numb, "", 0, 0, opts.resolution * 2, opts.resolution * 2); - pagecnv->Divide(opts.plotsperpage, opts.plotsperpage); - for (int i = 1; i <= opts.plotsperpage*opts.plotsperpage; i++) - if (it != datapullscanvaslist.end()) - { - pagecnv->cd(i); - (*it)->DrawClonePad(); - it++; - } - } + { + pagecnv = new TCanvas(numb, "", 0, 0, opts.resolution * 2, opts.resolution * 2); + pagecnv->Divide(opts.plotsperpage, opts.plotsperpage); + for (int i = 1; i <= opts.plotsperpage*opts.plotsperpage; i++) + if (it != datapullscanvaslist.end()) + { + pagecnv->cd(i); + (*it)->DrawClonePad(); + it++; + } + } pgn++; sprintf(pgnum, "%d", pgn); pagecnv->Print((opts.outdir + "plots_" + pgnum + ".eps").c_str()); @@ -230,12 +230,12 @@ int main(int argc, char **argv) pagecnv = new TCanvas(numb, "", 0, 0, opts.resolution * 2, opts.resolution * 2); pagecnv->Divide(2, 2); for (int i = 1; i <= 4; i++) - if (it != chi2scancanvaslist.end()) - { - pagecnv->cd(i); - (*it)->DrawClonePad(); - it++; - } + if (it != chi2scancanvaslist.end()) + { + pagecnv->cd(i); + (*it)->DrawClonePad(); + it++; + } pgn++; sprintf(pgnum, "%d", pgn); pagecnv->Print((opts.outdir + "plots_" + pgnum + ".eps").c_str()); @@ -245,36 +245,36 @@ int main(int argc, char **argv) { string ext = opts.ext; if (opts.ext == "pdf") - ext = "eps"; + ext = "eps"; gStyle->SetPaperSize(opts.pagewidth / 2., opts.pagewidth / 2.); for (vector <TCanvas*>::iterator it = pdfscanvaslist.begin(); it != pdfscanvaslist.end(); it++) - (*it)->Print((opts.outdir + (*it)->GetName() + "." + ext).c_str()); + (*it)->Print((opts.outdir + (*it)->GetName() + "." + ext).c_str()); for (vector <TCanvas*>::iterator it = pdfscanvasratiolist.begin(); it != pdfscanvasratiolist.end(); it++) - (*it)->Print((opts.outdir + (*it)->GetName() + "." + ext).c_str()); + (*it)->Print((opts.outdir + (*it)->GetName() + "." + ext).c_str()); if (opts.twopanels || opts.threepanels) - gStyle->SetPaperSize(opts.pagewidth, opts.pagewidth / 2.); + gStyle->SetPaperSize(opts.pagewidth, opts.pagewidth / 2.); for (vector <TCanvas*>::iterator it = datapullscanvaslist.begin(); it != datapullscanvaslist.end(); it++) - (*it)->Print((opts.outdir + (*it)->GetName() + "." + ext).c_str()); + (*it)->Print((opts.outdir + (*it)->GetName() + "." + ext).c_str()); gStyle->SetPaperSize(opts.pagewidth, opts.pagewidth); for (vector <TCanvas*>::iterator it = shiftcanvaslist.begin(); it != shiftcanvaslist.end(); it++) - (*it)->Print((opts.outdir + (*it)->GetName() + "." + ext).c_str()); + (*it)->Print((opts.outdir + (*it)->GetName() + "." + ext).c_str()); for (vector <TCanvas*>::iterator it = chi2scancanvaslist.begin(); it != chi2scancanvaslist.end(); it++) - (*it)->Print((opts.outdir + (*it)->GetName() + "." + ext).c_str()); + (*it)->Print((opts.outdir + (*it)->GetName() + "." + ext).c_str()); if (opts.ext == "pdf") - { - for (vector <TCanvas*>::iterator it = pdfscanvaslist.begin(); it != pdfscanvaslist.end(); it++) - system(((string)"ps2pdf -dEPSCrop " + opts.outdir + (*it)->GetName() + ".eps " + opts.outdir + (*it)->GetName() + ".pdf").c_str()); - for (vector <TCanvas*>::iterator it = pdfscanvasratiolist.begin(); it != pdfscanvasratiolist.end(); it++) - system(((string)"ps2pdf -dEPSCrop " + opts.outdir + (*it)->GetName() + ".eps " + opts.outdir + (*it)->GetName() + ".pdf").c_str()); - for (vector <TCanvas*>::iterator it = datapullscanvaslist.begin(); it != datapullscanvaslist.end(); it++) - system(((string)"ps2pdf -dEPSCrop " + opts.outdir + (*it)->GetName() + ".eps " + opts.outdir + (*it)->GetName() + ".pdf").c_str()); - for (vector <TCanvas*>::iterator it = shiftcanvaslist.begin(); it != shiftcanvaslist.end(); it++) - system(((string)"ps2pdf -dEPSCrop " + opts.outdir + (*it)->GetName() + ".eps " + opts.outdir + (*it)->GetName() + ".pdf").c_str()); - for (vector <TCanvas*>::iterator it = chi2scancanvaslist.begin(); it != chi2scancanvaslist.end(); it++) - system(((string)"ps2pdf -dEPSCrop " + opts.outdir + (*it)->GetName() + ".eps " + opts.outdir + (*it)->GetName() + ".pdf").c_str()); - } + { + for (vector <TCanvas*>::iterator it = pdfscanvaslist.begin(); it != pdfscanvaslist.end(); it++) + system(((string)"ps2pdf -dEPSCrop " + opts.outdir + (*it)->GetName() + ".eps " + opts.outdir + (*it)->GetName() + ".pdf").c_str()); + for (vector <TCanvas*>::iterator it = pdfscanvasratiolist.begin(); it != pdfscanvasratiolist.end(); it++) + system(((string)"ps2pdf -dEPSCrop " + opts.outdir + (*it)->GetName() + ".eps " + opts.outdir + (*it)->GetName() + ".pdf").c_str()); + for (vector <TCanvas*>::iterator it = datapullscanvaslist.begin(); it != datapullscanvaslist.end(); it++) + system(((string)"ps2pdf -dEPSCrop " + opts.outdir + (*it)->GetName() + ".eps " + opts.outdir + (*it)->GetName() + ".pdf").c_str()); + for (vector <TCanvas*>::iterator it = shiftcanvaslist.begin(); it != shiftcanvaslist.end(); it++) + system(((string)"ps2pdf -dEPSCrop " + opts.outdir + (*it)->GetName() + ".eps " + opts.outdir + (*it)->GetName() + ".pdf").c_str()); + for (vector <TCanvas*>::iterator it = chi2scancanvaslist.begin(); it != chi2scancanvaslist.end(); it++) + system(((string)"ps2pdf -dEPSCrop " + opts.outdir + (*it)->GetName() + ".eps " + opts.outdir + (*it)->GetName() + ".pdf").c_str()); + } cout << "Multiple " << opts.ext << " plots saved in: " << (opts.outdir + "*." + opts.ext) << endl; } @@ -285,21 +285,21 @@ int main(int argc, char **argv) f->mkdir("Canvas"); f->cd("Canvas"); for (vector <TCanvas*>::iterator it = pdfscanvaslist.begin(); it != pdfscanvaslist.end(); it++) - (*it)->Write(); + (*it)->Write(); for (vector <TCanvas*>::iterator it = pdfscanvasratiolist.begin(); it != pdfscanvasratiolist.end(); it++) - (*it)->Write(); + (*it)->Write(); for (vector <TCanvas*>::iterator it = datapullscanvaslist.begin(); it != datapullscanvaslist.end(); it++) - (*it)->Write(); + (*it)->Write(); for (vector <TCanvas*>::iterator it = shiftcanvaslist.begin(); it != shiftcanvaslist.end(); it++) - (*it)->Write(); + (*it)->Write(); for (vector <TCanvas*>::iterator it = chi2scancanvaslist.begin(); it != chi2scancanvaslist.end(); it++) - (*it)->Write(); + (*it)->Write(); f->cd(""); f->mkdir("Graphs"); f->cd("Graphs"); for (vector <TGraphAsymmErrors*>::iterator git = allgraphs.begin(); git != allgraphs.end(); git++) - (*git)->Write(); + (*git)->Write(); f->Close(); cout << "TCanvas saved in: " << (opts.outdir + "plots.root") << endl;