GenericPDF.h

Go to the documentation of this file.

#pragma once
#include "PDFxTMDLib/Common/ConfigWrapper.h"
#include "PDFxTMDLib/Common/Exception.h"
#include "PDFxTMDLib/Common/PDFUtils.h"
#include "PDFxTMDLib/Common/PartonUtils.h"
#include "PDFxTMDLib/Implementation/Extrapolator/Collinear/DPD/CDPDZeroExtrapolator.h"
#include "PDFxTMDLib/Implementation/Extrapolator/Collinear/SPDF/CContinuationExtrapolator.h"
#include "PDFxTMDLib/Implementation/Extrapolator/TMD/SPDF/TZeroExtrapolator.h"
#include "PDFxTMDLib/Implementation/Interpolator/Collinear/DPD/CPDFxTMDDPDInterpolator.h"
#include "PDFxTMDLib/Implementation/Interpolator/Collinear/SPDF/CLHAPDFBicubicInterpolator.h"
#include "PDFxTMDLib/Implementation/Interpolator/TMD/SPDF/TTrilinearInterpolator.h"
#include "PDFxTMDLib/Implementation/Reader/Collinear/DPD/CDefaultDPDReader.h"
#include "PDFxTMDLib/Implementation/Reader/Collinear/SPDF/CDefaultLHAPDFFileReader.h"
#include "PDFxTMDLib/Implementation/Reader/TMD/SPDF/TDefaultLHAPDF_TMDReader.h"
#include "PDFxTMDLib/Interface/IExtrapolator.h"
#include "PDFxTMDLib/Interface/IReader.h"
#include "PDFxTMDLib/Interface/SPDF/ITMDInterpolator.h"
#include "PDFxTMDLib/Interface/SPDF/IcPDFInterpolator.h"
#include <PDFxTMDLib/Common/YamlMetaInfo/YamlStandardPDFInfo.h>
#include <algorithm>
#include <array>
#include <iostream>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <vector>
 
namespace PDFxTMD
{
struct TMDPDFTag
{
};
struct CollinearPDFTag
{
};
struct CollinearDPDTag
{
};
 
// Type trait to get default implementations based on tag
template <typename Tag> struct DefaultPDFImplementations;
 
// Specialization for TMDPDFTag
template <> struct DefaultPDFImplementations<TMDPDFTag>
{
    using Reader = TDefaultLHAPDF_TMDReader;
    using Interpolator = TTrilinearInterpolator<Reader>;
    using Extrapolator = TZeroExtrapolator;
};
 
// Specialization for CollinearPDFTag
template <> struct DefaultPDFImplementations<CollinearPDFTag>
{
    using Reader = CDefaultLHAPDFFileReader;
    using Interpolator = CLHAPDFBicubicInterpolator<Reader>;
    using Extrapolator = CContinuationExtrapolator<CLHAPDFBicubicInterpolator<Reader>>;
};
 
template <> struct DefaultPDFImplementations<CollinearDPDTag>
{
    using Reader = CDefaultDPDReader;
    using Interpolator = CPDFxTMDDPDInterpolator<Reader>;
    using Extrapolator = CDPDZeroExtrapolator;
};
 
template <typename Tag, typename Reader = typename DefaultPDFImplementations<Tag>::Reader,
          typename Interpolator = typename DefaultPDFImplementations<Tag>::Interpolator,
          typename Extrapolator = typename DefaultPDFImplementations<Tag>::Extrapolator>
class GenericPDF
{
  public:
    GenericPDF(const std::string &pdfName, int setNumber)
        : m_pdfName(pdfName), m_setNumber(setNumber)
    {
        loadStandardInfo();
        loadData();
    }
    ~GenericPDF() = default;
    double pdf(PartonFlavor flavor, double x, double mu2) const
    {
        if constexpr (std::is_same_v<Tag, CollinearPDFTag>)
        {
            if (isInRange(m_reader, x, mu2))
                return m_interpolator.interpolate(flavor, x, mu2);
            return m_extrapolator.extrapolate(flavor, x, mu2);
        }
        else
        {
            throw std::logic_error(
                "pdf(PartonFlavor, double, double) is not supported for this tag.");
        }
    }
    void pdf(double x, double mu2, std::array<double, DEFAULT_TOTAL_PDFS> &output) const
    {
        if constexpr (std::is_same_v<Tag, CollinearPDFTag>)
        {
            if (isInRange(m_reader, x, mu2))
                return m_interpolator.interpolate(x, mu2, output);
            return m_extrapolator.extrapolate(x, mu2, output);
        }
        else
        {
            throw std::logic_error("pdf(double, double, std::array<double, DEFAULT_TOTAL_PDFS>&) "
                                   "is not supported for this tag.");
        }
    }
    GenericPDF(GenericPDF &&other) noexcept
        : m_pdfName(std::move(other.m_pdfName)), m_setNumber(other.m_setNumber),
          m_reader(std::move(other.m_reader)), m_interpolator(std::move(other.m_interpolator)),
          m_extrapolator(std::move(other.m_extrapolator)), m_stdInfo(std::move(other.m_stdInfo))
    {
        m_interpolator.initialize(&m_reader);
        if constexpr (std::is_base_of_v<IcAdvancedPDFExtrapolator<Extrapolator>, Extrapolator>)
        {
            m_extrapolator.setInterpolator(&m_interpolator);
        }
    }
    GenericPDF &operator=(GenericPDF &&other) noexcept
    {
        if (this == &other)
        {
            return *this;
        }
        m_pdfName = std::move(other.m_pdfName);
        m_setNumber = other.m_setNumber;
        m_reader = std::move(other.m_reader);
        m_interpolator = std::move(other.m_interpolator);
        m_extrapolator = std::move(other.m_extrapolator);
        m_stdInfo = std::move(other.m_stdInfo);
 
        // Re-initialize internal dependencies
        m_interpolator.initialize(&m_reader);
        if constexpr (std::is_base_of_v<IcAdvancedPDFExtrapolator<Extrapolator>, Extrapolator>)
        {
            m_extrapolator.setInterpolator(&m_interpolator);
        }
        return *this;
    }
 
    GenericPDF(const GenericPDF &other)
        : m_pdfName(other.m_pdfName), m_setNumber(other.m_setNumber), m_reader(other.m_reader),
          m_interpolator(other.m_interpolator), m_extrapolator(other.m_extrapolator),
          m_stdInfo(other.m_stdInfo)
    {
        m_interpolator.initialize(&m_reader);
        if constexpr (std::is_base_of_v<IcAdvancedPDFExtrapolator<Extrapolator>, Extrapolator>)
        {
            m_extrapolator.setInterpolator(&m_interpolator);
        }
    }
    GenericPDF &operator=(const GenericPDF &other)
    {
        if (this == &other)
        {
            return *this;
        }
        m_pdfName = other.m_pdfName;
        m_setNumber = other.m_setNumber;
        m_reader = other.m_reader;
        m_interpolator = other.m_interpolator;
        m_extrapolator = other.m_extrapolator;
        m_stdInfo = other.m_stdInfo;
 
        // Re-initialize internal dependencies
        m_interpolator.initialize(&m_reader);
        if constexpr (std::is_base_of_v<IcAdvancedPDFExtrapolator<Extrapolator>, Extrapolator>)
        {
            m_extrapolator.setInterpolator(&m_interpolator);
        }
        return *this;
    }
    double tmd(PartonFlavor flavor, double x, double kt2, double mu2)
    {
        if constexpr (std::is_same_v<Tag, TMDPDFTag>)
        {
            if (isInRange(m_reader, x, kt2, mu2))
                return m_interpolator.interpolate(flavor, x, kt2, mu2);
 
            return m_extrapolator.extrapolate(flavor, x, kt2, mu2);
        }
        else
        {
            throw std::logic_error(
                "pdf(double, double, std::array<double, 13>&) is not supported for this tag.");
        }
    }
    void tmd(double x, double kt2, double mu2, std::array<double, DEFAULT_TOTAL_PDFS> &output)
    {
        if constexpr (std::is_same_v<Tag, TMDPDFTag>)
        {
            isInRange(m_reader, x, kt2, mu2) ? m_interpolator.interpolate(x, kt2, mu2, output)
                                             : m_extrapolator.extrapolate(x, kt2, mu2, output);
        }
        else
        {
            throw std::logic_error(
                "pdf(double, double, std::array<double, 13>&) is not supported for this tag.");
        }
    }
    double dpd(PartonFlavor flavor1, PartonFlavor flavor2, double x1, double mu1_2, double x2,
               double mu2_2)
    {
        if constexpr (std::is_same_v<Tag, CollinearDPDTag>)
        {
            // note: my current interpolator supports the case where x is larger than the greatest x
            // in the grid. it should be handeled later.But is enough for most of the analysis.
            // because x > 0.97.. has zero contriution in most cases
            if (isInRange(m_reader, x1, mu1_2) && isInRange(m_reader, x2, mu2_2))
                return m_interpolator.interpolate(flavor1, flavor2, x1, mu1_2, x2, mu2_2);
            return m_extrapolator.extrapolate(flavor1, flavor2, x1, mu1_2, x2, mu2_2);
        }
        else
        {
            throw std::logic_error(
                "pdf(PartonFlavor, double, double) is not supported for this tag.");
        }
    }
    YamlStandardTMDInfo getStdPDFInfo() const
    {
        return m_stdInfo;
    }
 
  private:
    void loadStandardInfo()
    {
        auto infoPathPair = StandardInfoFilePath(m_pdfName);
        if (infoPathPair.second != ErrorType::None)
            throw FileLoadException("Unable to find info file of PDF set " + m_pdfName);
        if constexpr (std::is_same_v<Tag, TMDPDFTag> || std::is_same_v<Tag, CollinearPDFTag> ||
                      std::is_same_v<Tag, CollinearDPDTag>)
        {
            auto pdfStandardInfo = YamlStandardPDFInfoReader(*infoPathPair.first);
            if (pdfStandardInfo.second != ErrorType::None)
                throw InvalidFormatException("Invalid standard info file " + *infoPathPair.first);
            m_stdInfo = *pdfStandardInfo.first;
        }
        else
        {
            static_assert(!std::is_same_v<Tag, Tag>, "Unsupported Tag");
        }
    }
    void loadData()
    {
        m_reader.read(m_pdfName, m_setNumber);
        m_interpolator.initialize(&m_reader);
        if constexpr (std::is_same_v<Tag, CollinearPDFTag>)
        {
            if constexpr (std::is_base_of_v<IcAdvancedPDFExtrapolator<Extrapolator>, Extrapolator>)
            {
                m_extrapolator.setInterpolator(&m_interpolator);
            }
        }
    }
    std::string m_pdfName;
    int m_setNumber;
    Reader m_reader;
    Interpolator m_interpolator;
    Extrapolator m_extrapolator;
    YamlStandardTMDInfo m_stdInfo;
};
 
// Convenient type aliases for common use cases
using TMDPDF = GenericPDF<TMDPDFTag>;
using CollinearPDF = GenericPDF<CollinearPDFTag>;
using CollinearDPD = GenericPDF<CollinearDPDTag>;
 
} // namespace PDFxTMD