Heppy/src/PdfWeightProducerTool.cc

0001 #include "PhysicsTools/Heppy/interface/PdfWeightProducerTool.h"
0002 #include <cassert>
0003
0004 namespace LHAPDF {
0005   void initPDFSet(int nset, const std::string &filename, int member = 0);
0006   int numberPDF(int nset);
0007   void usePDFMember(int nset, int member);
0008   double xfx(int nset, double x, double Q, int fl);
0009   double getXmin(int nset, int member);
0010   double getXmax(int nset, int member);
0011   double getQ2min(int nset, int member);
0012   double getQ2max(int nset, int member);
0013   void extrapolate(bool extrapolate = true);
0014 }  // namespace LHAPDF
0015
0016 namespace heppy {
0017
0018   void PdfWeightProducerTool::addPdfSet(const std::string &name) {
0019     pdfs_.push_back(name);
0020     weights_[name] = std::vector<double>();
0021   }
0022
0023   void PdfWeightProducerTool::beginJob() {
0024     for (unsigned int i = 0, n = pdfs_.size(); i < n; ++i) {
0025       LHAPDF::initPDFSet(i + 1, pdfs_[i]);
0026     }
0027   }
0028
0029   void PdfWeightProducerTool::processEvent(const GenEventInfoProduct &pdfstuff) {
0030     float Q = pdfstuff.pdf()->scalePDF;
0031
0032     int id1 = pdfstuff.pdf()->id.first;
0033     double x1 = pdfstuff.pdf()->x.first;
0034     double pdf1 = pdfstuff.pdf()->xPDF.first;
0035
0036     int id2 = pdfstuff.pdf()->id.second;
0037     double x2 = pdfstuff.pdf()->x.second;
0038     double pdf2 = pdfstuff.pdf()->xPDF.second;
0039
0040     for (unsigned int i = 0, n = pdfs_.size(); i < n; ++i) {
0041       std::vector<double> &weights = weights_[pdfs_[i]];
0042       unsigned int nweights = 1;
0043       if (LHAPDF::numberPDF(i + 1) > 1)
0044         nweights += LHAPDF::numberPDF(i + 1);
0045       weights.resize(nweights);
0046
0047       for (unsigned int j = 0; j < nweights; ++j) {
0048         LHAPDF::usePDFMember(i + 1, j);
0049         double newpdf1 = LHAPDF::xfx(i + 1, x1, Q, id1) / x1;
0050         double newpdf2 = LHAPDF::xfx(i + 1, x2, Q, id2) / x2;
0051         weights[j] = newpdf1 / pdf1 * newpdf2 / pdf2;
0052       }
0053     }
0054   }
0055
0056   const std::vector<double> &PdfWeightProducerTool::getWeights(const std::string &name) const {
0057     std::map<std::string, std::vector<double> >::const_iterator match = weights_.find(name);
0058     assert(match != weights_.end());
0059     return match->second;
0060   }
0061
0062 }  // namespace heppy