PdfWeightProducerTool.cc

Go to the documentation of this file.

#include "PhysicsTools/Heppy/interface/PdfWeightProducerTool.h"
#include <cassert>

namespace LHAPDF {
  void initPDFSet(int nset, const std::string &filename, int member = 0);
  int numberPDF(int nset);
  void usePDFMember(int nset, int member);
  double xfx(int nset, double x, double Q, int fl);
  double getXmin(int nset, int member);
  double getXmax(int nset, int member);
  double getQ2min(int nset, int member);
  double getQ2max(int nset, int member);
  void extrapolate(bool extrapolate = true);
}  // namespace LHAPDF

namespace heppy {

  void PdfWeightProducerTool::addPdfSet(const std::string &name) {
    pdfs_.push_back(name);
    weights_[name] = std::vector<double>();
  }

  void PdfWeightProducerTool::beginJob() {
    for (unsigned int i = 0, n = pdfs_.size(); i < n; ++i) {
      LHAPDF::initPDFSet(i + 1, pdfs_[i]);
    }
  }

  void PdfWeightProducerTool::processEvent(const GenEventInfoProduct &pdfstuff) {
    float Q = pdfstuff.pdf()->scalePDF;

    int id1 = pdfstuff.pdf()->id.first;
    double x1 = pdfstuff.pdf()->x.first;
    double pdf1 = pdfstuff.pdf()->xPDF.first;

    int id2 = pdfstuff.pdf()->id.second;
    double x2 = pdfstuff.pdf()->x.second;
    double pdf2 = pdfstuff.pdf()->xPDF.second;

    for (unsigned int i = 0, n = pdfs_.size(); i < n; ++i) {
      std::vector<double> &weights = weights_[pdfs_[i]];
      unsigned int nweights = 1;
      if (LHAPDF::numberPDF(i + 1) > 1)
        nweights += LHAPDF::numberPDF(i + 1);
      weights.resize(nweights);

      for (unsigned int j = 0; j < nweights; ++j) {
        LHAPDF::usePDFMember(i + 1, j);
        double newpdf1 = LHAPDF::xfx(i + 1, x1, Q, id1) / x1;
        double newpdf2 = LHAPDF::xfx(i + 1, x2, Q, id2) / x2;
        weights[j] = newpdf1 / pdf1 * newpdf2 / pdf2;
      }
    }
  }

  const std::vector<double> &PdfWeightProducerTool::getWeights(const std::string &name) const {
    std::map<std::string, std::vector<double> >::const_iterator match = weights_.find(name);
    assert(match != weights_.end());
    return match->second;
  }

}  // namespace heppy