PDFSet.h

// -*- C++ -*-
//
// This file is part of LHAPDF
// Copyright (C) 2012-2023 The LHAPDF collaboration (see AUTHORS for details)
//
#pragma once
#ifndef LHAPDF_PDFSet_H
#define LHAPDF_PDFSet_H

#include "LHAPDF/Info.h"
#include "LHAPDF/Factories.h"
#include "LHAPDF/Version.h"
#include "LHAPDF/Config.h"
#include "LHAPDF/Utils.h"

namespace LHAPDF {


  const double CL1SIGMA = 100*erf(1/sqrt(2));


  // Forward declaration
  class PDF;




  struct PDFUncertainty {
    using ErrPairs = std::vector<std::pair<double,double>>;

    PDFUncertainty(double cent=0, double eplus=0, double eminus=0, double esymm=0, double scalefactor=1,
                   double eplus_pdf=0, double eminus_pdf=0, double esymm_pdf=0,
                   double eplus_par=0, double eminus_par=0, double esymm_par=0)
      : central(cent), errplus(eplus), errminus(eminus), errsymm(esymm), scale(scalefactor),
        errplus_pdf(eplus_pdf), errminus_pdf(eminus_pdf), errsymm_pdf(esymm_pdf),
        errplus_par(eplus_par), errminus_par(eminus_par), errsymm_par(esymm_par)
    {    }

    double central, errplus, errminus, errsymm, scale;

    double errplus_pdf, errminus_pdf, errsymm_pdf;
    double errplus_par, errminus_par, errsymm_par;
    double err_par; 

    ErrPairs errparts;
  };



  struct PDFErrInfo {
    using EnvPart = std::pair<std::string, size_t>;
    using EnvParts = std::vector<EnvPart>;
    using QuadParts = std::vector<EnvParts>;

    PDFErrInfo(QuadParts parts, double cl, const std::string& errtypestr="")
      : qparts(parts), conflevel(cl), errtype(errtypestr)
    {    }

    PDFErrInfo() {}

    QuadParts qparts;

    double conflevel;

    std::string errtype;

    std::string coreType() const { return qpartName(0); }

    std::string qpartName(size_t iq) const;
    std::vector<std::string> qpartNames() const;

    size_t nmemCore() const;
    size_t nmemPar() const;

  };




  class PDFSet : public Info {
  public:

    PDFSet() { }

    PDFSet(const std::string& setname);



    std::string name() const {
      return _setname;
    }

    std::string description() const {
      return get_entry("SetDesc");
    }

    int lhapdfID() const {
      return get_entry_as<int>("SetIndex", -1);
    }

    int dataversion() const {
      return get_entry_as<int>("DataVersion", -1);
    }

    std::string errorType() const {
      return to_lower(get_entry("ErrorType", "UNKNOWN"));
    }

    PDFErrInfo errorInfo() const;

    double errorConfLevel() const;

    // int numMembers() const {
    //   return get_entry_as<int>("NumMembers");
    // }
    size_t size() const {
      return get_entry_as<unsigned int>("NumMembers");
    }

    size_t errorSize() const {
      return size()-1;
    }
    size_t errSize() const {
      return errorSize();
    }



    void print(std::ostream& os=std::cout, int verbosity=1) const;



    PDF* mkPDF(size_t member) const {
      return LHAPDF::mkPDF(name(), member);
    }


    //
    template <typename PTR>
    void mkPDFs(std::vector<PTR>& pdfs) const {
      const int v = verbosity();
      if (v > 0) {
        std::cout << "LHAPDF " << version() << " loading all " << size() << " PDFs in set " << name() << std::endl;
        this->print(std::cout, v);
        if (this->has_key("Note")) std::cout << get_entry("Note") << std::endl;
      }
      pdfs.clear();
      pdfs.reserve(size());
      if (v < 2) setVerbosity(0); //< Disable every-member printout unless verbosity level is high
      for (size_t i = 0; i < size(); ++i) {
        pdfs.push_back( PTR(mkPDF(i)) );
      }
      setVerbosity(v);
    }

    std::vector<PDF*> mkPDFs() const {
      std::vector<PDF*> rtn;
      mkPDFs(rtn);
      return rtn;
    }

    // template <typename PTR=PDF*>
    template <typename PTR>
    std::vector<PTR> mkPDFs() const {
      std::vector<PTR> rtn;
      mkPDFs(rtn);
      return rtn;
    }






    std::vector<std::string> keys() const {
      std::vector<std::string> rtn = getConfig().keys();
      for (const std::string& k : keys_local()) {
        if (!contains(rtn, k)) rtn.push_back(k);
      }
      return rtn;
    }

    bool has_key(const std::string& key) const {
      return has_key_local(key) || getConfig().has_key(key);
    }

    const std::string& get_entry(const std::string& key) const {
      if (has_key_local(key)) return get_entry_local(key); //< value is defined locally
      return getConfig().get_entry(key); //< fall back to the global config
    }

    const std::string& get_entry(const std::string& key, const std::string& fallback) const {
      return Info::get_entry(key, fallback);
    }




    PDFUncertainty uncertainty(const std::vector<double>& values,
                               double cl=CL1SIGMA, bool alternative=false) const {
      PDFUncertainty rtn;
      uncertainty(rtn, values, cl, alternative);
      return rtn;
    }


    // // Trick to ensure no calls with implicit type conversion
    // template <typename T1, typename T2>
    // void uncertainty(const std::vector<double>& values, T1, T2) const = delete;

    // /// Alternative form allowing the alternative computation with default CL
    // PDFUncertainty uncertainty(const std::vector<double>& values,
    //                            bool alternative, double cl=CL1SIGMA) const {
    //   return uncertainty(values, cl, alternative);
    // }


    void uncertainty(PDFUncertainty& rtn,
                     const std::vector<double>& values,
                     double cl=CL1SIGMA, bool alternative=false) const;

    // // Trick to ensure no calls with implicit type conversion
    // template <typename T1, typename T2>
    // void uncertainty(PDFUncertainty& rtn, const std::vector<double>& values, T1, T2) const = delete;

    // /// Alternative form allowing the alternative computation with default CL
    // void uncertainty(PDFUncertainty& rtn,
    //                  const std::vector<double>& values,
    //                  bool alternative, double cl=CL1SIGMA) const {
    //   uncertainty(rtn, values, cl, alternative);
    // }


    std::vector<PDFUncertainty> uncertainties(const std::vector<std::vector<double>>& observables_values,
                                              double cl=CL1SIGMA, bool alternative=false) const {
      std::vector<PDFUncertainty> rtn;
      uncertainties(rtn, observables_values, cl, alternative);
      return rtn;
    }


    void uncertainties(std::vector<PDFUncertainty>& rtn,
                       const std::vector<std::vector<double>>& observables_values,
                       double cl=CL1SIGMA, bool alternative=false) const;


    double correlation(const std::vector<double>& valuesA, const std::vector<double>& valuesB) const;

    double randomValueFromHessian(const std::vector<double>& values, const std::vector<double>& randoms, bool symmetrise=true) const;


    void _checkPdfType(const std::vector<string>& pdftypes) const;



  private:

    std::string _setname;

    mutable PDFErrInfo _errinfo;

  };


}
#endif