/ElectroWeakAnalysis/Utilities/src/PdfWeightProducer.cc

  1#include <memory>
  2
  3#include "FWCore/Framework/interface/Frameworkfwd.h"
  4#include "FWCore/Framework/interface/EDProducer.h"
  5
  6#include "FWCore/Framework/interface/Event.h"
  7#include "FWCore/Framework/interface/MakerMacros.h"
  8
  9#include "FWCore/ParameterSet/interface/ParameterSet.h"
 10
 11#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
 12#include "DataFormats/Candidate/interface/CompositeCandidate.h"
 13#include "CommonTools/CandUtils/interface/AddFourMomenta.h"
 14#include "CommonTools/CandUtils/interface/Booster.h"
 15#include <Math/VectorUtil.h>
 16
 17//
 18// class declaration
 19//
 20class PdfWeightProducer : public edm::EDProducer {
 21   public:
 22      explicit PdfWeightProducer(const edm::ParameterSet&);
 23      ~PdfWeightProducer();
 24
 25   private:
 26      virtual void beginJob() override ;
 27      virtual void produce(edm::Event&, const edm::EventSetup&) override;
 28      virtual void endJob() override ;
 29
 30      std::string fixPOWHEG_;
 31      bool useFirstAsDefault_;
 32      edm::InputTag genTag_;
 33      edm::InputTag pdfInfoTag_;
 34      std::vector<std::string> pdfSetNames_;
 35      std::vector<std::string> pdfShortNames_;
 36};
 37
 38#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
 39namespace LHAPDF {
 40      void initPDFSet(int nset, const std::string& filename, int member=0);
 41      int numberPDF(int nset);
 42      void usePDFMember(int nset, int member);
 43      double xfx(int nset, double x, double Q, int fl);
 44      double getXmin(int nset, int member);
 45      double getXmax(int nset, int member);
 46      double getQ2min(int nset, int member);
 47      double getQ2max(int nset, int member);
 48      void extrapolate(bool extrapolate=true);
 49}
 50
 51/////////////////////////////////////////////////////////////////////////////////////
 52PdfWeightProducer::PdfWeightProducer(const edm::ParameterSet& pset) :
 53 fixPOWHEG_(pset.getUntrackedParameter<std::string> ("FixPOWHEG", "")),
 54 useFirstAsDefault_(pset.getUntrackedParameter<bool>("useFirstAsDefault",false)),
 55 genTag_(pset.getUntrackedParameter<edm::InputTag> ("GenTag", edm::InputTag("genParticles"))),
 56 pdfInfoTag_(pset.getUntrackedParameter<edm::InputTag> ("PdfInfoTag", edm::InputTag("generator"))),
 57 pdfSetNames_(pset.getUntrackedParameter<std::vector<std::string> > ("PdfSetNames"))
 58{
 59      if (fixPOWHEG_ != "") pdfSetNames_.insert(pdfSetNames_.begin(),fixPOWHEG_);
 60
 61      if (pdfSetNames_.size()>3) {
 62            edm::LogWarning("") << pdfSetNames_.size() << " PDF sets requested on input. Using only the first 3 sets and ignoring the rest!!";
 63            pdfSetNames_.erase(pdfSetNames_.begin()+3,pdfSetNames_.end());
 64      }
 65
 66      for (unsigned int k=0; k<pdfSetNames_.size(); k++) {
 67            size_t dot = pdfSetNames_[k].find_first_of('.');
 68            size_t underscore = pdfSetNames_[k].find_first_of('_');
 69            if (underscore<dot) {
 70                  pdfShortNames_.push_back(pdfSetNames_[k].substr(0,underscore));
 71            } else {
 72                  pdfShortNames_.push_back(pdfSetNames_[k].substr(0,dot));
 73            }
 74            produces<std::vector<double> >(pdfShortNames_[k].data());
 75      }
 76} 
 77
 78/////////////////////////////////////////////////////////////////////////////////////
 79PdfWeightProducer::~PdfWeightProducer(){}
 80
 81/////////////////////////////////////////////////////////////////////////////////////
 82void PdfWeightProducer::beginJob() {
 83      for (unsigned int k=1; k<=pdfSetNames_.size(); k++) {
 84            LHAPDF::initPDFSet(k,pdfSetNames_[k-1]);
 85      }
 86}
 87
 88/////////////////////////////////////////////////////////////////////////////////////
 89void PdfWeightProducer::endJob(){}
 90
 91/////////////////////////////////////////////////////////////////////////////////////
 92void PdfWeightProducer::produce(edm::Event& iEvent, const edm::EventSetup&) {
 93
 94      if (iEvent.isRealData()) return;
 95
 96      edm::Handle<GenEventInfoProduct> pdfstuff;
 97      if (!iEvent.getByLabel(pdfInfoTag_, pdfstuff)) {
 98            edm::LogError("PDFWeightProducer") << ">>> PdfInfo not found: " << pdfInfoTag_.encode() << " !!!";
 99            return;
100      }
101
102      float Q = pdfstuff->pdf()->scalePDF;
103
104      int id1 = pdfstuff->pdf()->id.first;
105      double x1 = pdfstuff->pdf()->x.first;
106      double pdf1 = pdfstuff->pdf()->xPDF.first;
107
108      int id2 = pdfstuff->pdf()->id.second;
109      double x2 = pdfstuff->pdf()->x.second;
110      double pdf2 = pdfstuff->pdf()->xPDF.second; 
111      if (useFirstAsDefault_ && pdf1 == -1. && pdf2 == -1. ) {
112         LHAPDF::usePDFMember(1,0);
113         pdf1 = LHAPDF::xfx(1, x1, Q, id1)/x1;
114         pdf2 = LHAPDF::xfx(1, x2, Q, id2)/x2;
115      }
116
117      // Ad-hoc fix for POWHEG
118      if (fixPOWHEG_!="") {
119            edm::Handle<reco::GenParticleCollection> genParticles;
120            if (!iEvent.getByLabel(genTag_, genParticles)) {
121                  edm::LogError("PDFWeightProducer") << ">>> genParticles  not found: " << genTag_.encode() << " !!!";
122                  return;
123            }
124            unsigned int gensize = genParticles->size();
125            double mboson = 0.;
126            for(unsigned int i = 0; i<gensize; ++i) {
127                  const reco::GenParticle& part = (*genParticles)[i];
128                  int status = part.status();
129                  if (status!=3) continue;
130                  int id = part.pdgId();
131                  if (id!=23 && abs(id)!=24) continue;
132                  mboson = part.mass();
133                  break;
134            }
135            Q = sqrt(mboson*mboson+Q*Q);
136            LHAPDF::usePDFMember(1,0);
137            pdf1 = LHAPDF::xfx(1, x1, Q, id1)/x1;
138            pdf2 = LHAPDF::xfx(1, x2, Q, id2)/x2;
139      }
140
141      // Put PDF weights in the event
142      for (unsigned int k=1; k<=pdfSetNames_.size(); ++k) {
143            std::auto_ptr<std::vector<double> > weights (new std::vector<double>);
144            unsigned int nweights = 1;
145            if (LHAPDF::numberPDF(k)>1) nweights += LHAPDF::numberPDF(k);
146            weights->reserve(nweights);
147        
148            for (unsigned int i=0; i<nweights; ++i) {
149                  LHAPDF::usePDFMember(k,i);
150                  double newpdf1 = LHAPDF::xfx(k, x1, Q, id1)/x1;
151                  double newpdf2 = LHAPDF::xfx(k, x2, Q, id2)/x2;
152                  weights->push_back(newpdf1/pdf1*newpdf2/pdf2);
153            }
154            iEvent.put(weights,pdfShortNames_[k-1]);
155      }
156}
157
158DEFINE_FWK_MODULE(PdfWeightProducer);