/ElectroWeakAnalysis/Utilities/src/PdfWeightProducer.cc

https://github.com/aivanov-cern/cmssw · C++ · 158 lines · 124 code · 24 blank · 10 comment · 23 complexity · a1f4ce0cd1ab0eb26e9b8f5e3aa8e109 MD5 · raw file

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