CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PdfWeightProducer.cc
Go to the documentation of this file.
1 #include <memory>
2 
5 
8 
10 
16 #include <Math/VectorUtil.h>
17 
18 //
19 // class declaration
20 //
22 public:
23  explicit PdfWeightProducer(const edm::ParameterSet&);
24  ~PdfWeightProducer() override;
25 
26 private:
27  void beginJob() override;
28  void produce(edm::Event&, const edm::EventSetup&) override;
29  void endJob() override;
30 
37  std::vector<std::string> pdfSetNames_;
38  std::vector<std::string> pdfShortNames_;
39 };
40 
41 namespace LHAPDF {
42  void initPDFSet(int nset, const std::string& filename, int member = 0);
43  int numberPDF(int nset);
44  void usePDFMember(int nset, int member);
45  double xfx(int nset, double x, double Q, int fl);
46  double getXmin(int nset, int member);
47  double getXmax(int nset, int member);
48  double getQ2min(int nset, int member);
49  double getQ2max(int nset, int member);
50  void extrapolate(bool extrapolate = true);
51 } // namespace LHAPDF
52 
55  : fixPOWHEG_(pset.getUntrackedParameter<std::string>("FixPOWHEG", "")),
56  useFirstAsDefault_(pset.getUntrackedParameter<bool>("useFirstAsDefault", false)),
57  genTag_(pset.getUntrackedParameter<edm::InputTag>("GenTag", edm::InputTag("genParticles"))),
59  pdfInfoTag_(pset.getUntrackedParameter<edm::InputTag>("PdfInfoTag", edm::InputTag("generator"))),
61  pdfSetNames_(pset.getUntrackedParameter<std::vector<std::string> >("PdfSetNames")) {
62  if (!fixPOWHEG_.empty())
63  pdfSetNames_.insert(pdfSetNames_.begin(), fixPOWHEG_);
64 
65  if (pdfSetNames_.size() > 3) {
66  edm::LogWarning("") << pdfSetNames_.size()
67  << " PDF sets requested on input. Using only the first 3 sets and ignoring the rest!!";
68  pdfSetNames_.erase(pdfSetNames_.begin() + 3, pdfSetNames_.end());
69  }
70 
71  for (unsigned int k = 0; k < pdfSetNames_.size(); k++) {
72  size_t dot = pdfSetNames_[k].find_first_of('.');
73  size_t underscore = pdfSetNames_[k].find_first_of('_');
74  if (underscore < dot) {
75  pdfShortNames_.push_back(pdfSetNames_[k].substr(0, underscore));
76  } else {
77  pdfShortNames_.push_back(pdfSetNames_[k].substr(0, dot));
78  }
79  produces<std::vector<double> >(pdfShortNames_[k]);
80  }
81 }
82 
85 
88  for (unsigned int k = 1; k <= pdfSetNames_.size(); k++) {
90  }
91 }
92 
95 
98  if (iEvent.isRealData())
99  return;
100 
102  if (!iEvent.getByToken(pdfInfoToken_, pdfstuff)) {
103  edm::LogError("PDFWeightProducer") << ">>> PdfInfo not found: " << pdfInfoTag_.encode() << " !!!";
104  return;
105  }
106 
107  float Q = pdfstuff->pdf()->scalePDF;
108 
109  int id1 = pdfstuff->pdf()->id.first;
110  double x1 = pdfstuff->pdf()->x.first;
111  double pdf1 = pdfstuff->pdf()->xPDF.first;
112 
113  int id2 = pdfstuff->pdf()->id.second;
114  double x2 = pdfstuff->pdf()->x.second;
115  double pdf2 = pdfstuff->pdf()->xPDF.second;
116  if (useFirstAsDefault_ && pdf1 == -1. && pdf2 == -1.) {
117  LHAPDF::usePDFMember(1, 0);
118  pdf1 = LHAPDF::xfx(1, x1, Q, id1) / x1;
119  pdf2 = LHAPDF::xfx(1, x2, Q, id2) / x2;
120  }
121 
122  // Ad-hoc fix for POWHEG
123  if (!fixPOWHEG_.empty()) {
125  if (!iEvent.getByToken(genToken_, genParticles)) {
126  edm::LogError("PDFWeightProducer") << ">>> genParticles not found: " << genTag_.encode() << " !!!";
127  return;
128  }
129  unsigned int gensize = genParticles->size();
130  double mboson = 0.;
131  for (unsigned int i = 0; i < gensize; ++i) {
132  const reco::GenParticle& part = (*genParticles)[i];
133  int status = part.status();
134  if (status != 3)
135  continue;
136  int id = part.pdgId();
137  if (id != 23 && abs(id) != 24)
138  continue;
139  mboson = part.mass();
140  break;
141  }
142  Q = sqrt(mboson * mboson + Q * Q);
143  LHAPDF::usePDFMember(1, 0);
144  pdf1 = LHAPDF::xfx(1, x1, Q, id1) / x1;
145  pdf2 = LHAPDF::xfx(1, x2, Q, id2) / x2;
146  }
147 
148  // Put PDF weights in the event
149  for (unsigned int k = 1; k <= pdfSetNames_.size(); ++k) {
150  std::unique_ptr<std::vector<double> > weights(new std::vector<double>);
151  unsigned int nweights = 1;
152  if (LHAPDF::numberPDF(k) > 1)
153  nweights += LHAPDF::numberPDF(k);
154  weights->reserve(nweights);
155 
156  for (unsigned int i = 0; i < nweights; ++i) {
158  double newpdf1 = LHAPDF::xfx(k, x1, Q, id1) / x1;
159  double newpdf2 = LHAPDF::xfx(k, x2, Q, id2) / x2;
160  weights->push_back(newpdf1 / pdf1 * newpdf2 / pdf2);
161  }
162  iEvent.put(std::move(weights), pdfShortNames_[k - 1]);
163  }
164 }
165 
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
int pdgId() const final
PDG identifier.
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:131
const PDF * pdf() const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:525
std::pair< double, double > x
Definition: PdfInfo.h:13
edm::EDGetTokenT< reco::GenParticleCollection > genToken_
edm::InputTag genTag_
bool isRealData() const
Definition: EventBase.h:62
std::string encode() const
Definition: InputTag.cc:159
void beginJob() override
void initPDFSet(int nset, const std::string &filename, int member=0)
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
std::vector< std::string > pdfSetNames_
std::pair< double, double > xPDF
Definition: PdfInfo.h:14
double getXmax(int nset, int member)
double getQ2max(int nset, int member)
int iEvent
Definition: GenABIO.cc:224
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
T sqrt(T t)
Definition: SSEVec.h:19
void extrapolate(bool extrapolate=true)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int numberPDF(int nset)
double getXmin(int nset, int member)
PdfWeightProducer(const edm::ParameterSet &)
edm::EDGetTokenT< GenEventInfoProduct > pdfInfoToken_
std::pair< int, int > id
Definition: PdfInfo.h:12
edm::InputTag pdfInfoTag_
part
Definition: HCALResponse.h:20
void produce(edm::Event &, const edm::EventSetup &) override
std::vector< std::string > pdfShortNames_
fixed size matrix
HLT enums.
T dot(const Basic3DVector &v) const
Scalar product, or "dot" product, with a vector of same type.
double xfx(int nset, double x, double Q, int fl)
int status() const final
status word
~PdfWeightProducer() override
double getQ2min(int nset, int member)
EDGetTokenT< ProductType > mayConsume(edm::InputTag const &tag)
double scalePDF
Definition: PdfInfo.h:15
def move(src, dest)
Definition: eostools.py:511
double mass() const final
mass
void usePDFMember(int nset, int member)
void endJob() override