LHEInterface/plugins/LHECOMWeightProducer.cc

0001 #include <iostream>
0002 #include <memory>
0003
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/one/EDProducer.h"
0006
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/Run.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012
0013 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0014 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
0015 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0016
0017 #include <sstream>
0018 //class that reweights a pure parton level event from the originale COM energy to
0019 //an energy that is < than original COM energy.
0020
0021 //
0022 // class declaration
0023 //
0024 class LHECOMWeightProducer : public edm::one::EDProducer<edm::one::WatchRuns> {
0025 public:
0026   explicit LHECOMWeightProducer(const edm::ParameterSet&);
0027   ~LHECOMWeightProducer() override = default;
0028
0029 private:
0030   void beginJob() override;
0031   void beginRun(edm::Run const& run, const edm::EventSetup& es) override;
0032   void endRun(edm::Run const& run, const edm::EventSetup& es) override {}
0033   void produce(edm::Event&, const edm::EventSetup&) override;
0034
0035   const edm::InputTag lheTag_;
0036   const double _newECMS;
0037   const edm::EDGetTokenT<LHERunInfoProduct> tokenLHERun_;
0038   const edm::EDGetTokenT<LHEEventProduct> tokenLHEEvent_;
0039   int _pdfset;
0040   int _pdfmember;
0041   double _origECMS;
0042   std::string _label;
0043 };
0044
0045 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0046 namespace LHAPDF {
0047   void initPDFSet(int nset, int setid, int member = 0);
0048   int numberPDF(int nset);
0049   void usePDFMember(int nset, int member);
0050   double xfx(int nset, double x, double Q, int fl);
0051   double getXmin(int nset, int member);
0052   double getXmax(int nset, int member);
0053   double getQ2min(int nset, int member);
0054   double getQ2max(int nset, int member);
0055   void extrapolate(bool extrapolate = true);
0056 }  // namespace LHAPDF
0057
0058 /////////////////////////////////////////////////////////////////////////////////////
0059 LHECOMWeightProducer::LHECOMWeightProducer(const edm::ParameterSet& pset)
0060     : lheTag_(pset.getParameter<edm::InputTag>("lheSrc")),
0061       _newECMS(pset.getParameter<double>("NewECMS")),
0062       tokenLHERun_(consumes<LHERunInfoProduct, edm::InRun>(lheTag_)),
0063       tokenLHEEvent_(consumes<LHEEventProduct>(lheTag_)) {
0064   std::stringstream labelStr;
0065   labelStr << "com"
0066            << "To" << _newECMS;
0067   _label = labelStr.str();
0068   produces<GenEventInfoProduct>(_label);
0069 }
0070
0071 /////////////////////////////////////////////////////////////////////////////////////
0072 void LHECOMWeightProducer::beginRun(edm::Run const& run, const edm::EventSetup& es) {
0073   using namespace edm;
0074   const edm::Handle<LHERunInfoProduct>& lheRun = run.getHandle(tokenLHERun_);
0075   //assumes the same pdf is used for both beams
0076   _pdfset = lheRun->heprup().PDFSUP.first;
0077   _pdfmember = lheRun->heprup().PDFGUP.first;
0078   _origECMS = lheRun->heprup().EBMUP.first + lheRun->heprup().EBMUP.second;
0079   edm::LogInfo("LHECOMWeightProducer") << "PDFSET: " << _pdfset << "; member: " << _pdfmember
0080                                        << "; COM energy: " << _origECMS;
0081   if (_newECMS > _origECMS)
0082     throw cms::Exception("LHECOMWeightProducer") << "You cannot reweight COM energy to a higher than original energy ";
0083   LHAPDF::initPDFSet(1, _pdfset, _pdfmember);
0084 }
0085
0086 /////////////////////////////////////////////////////////////////////////////////////
0087 void LHECOMWeightProducer::beginJob() {
0088   //LHAPDF::initPDFSet(1,pdfSetName_);
0089 }
0090
0091 /////////////////////////////////////////////////////////////////////////////////////
0092 void LHECOMWeightProducer::produce(edm::Event& iEvent, const edm::EventSetup&) {
0093   using namespace std;
0094   bool verbose = false;
0095
0096   if (iEvent.isRealData())
0097     return;
0098
0099   const edm::Handle<LHEEventProduct>& lheevent = iEvent.getHandle(tokenLHEEvent_);
0100
0101   float Q = lheevent->hepeup().SCALUP;
0102
0103   int id1 = lheevent->hepeup().IDUP[0];
0104   double x1 = fabs(lheevent->hepeup().PUP[0][2] / (_origECMS / 2));
0105   double x1prime = fabs(lheevent->hepeup().PUP[0][2] / (_newECMS / 2));
0106
0107   int id2 = lheevent->hepeup().IDUP[1];
0108   double x2 = fabs(lheevent->hepeup().PUP[1][2] / (_origECMS / 2));
0109   double x2prime = fabs(lheevent->hepeup().PUP[1][2] / (_newECMS / 2));
0110
0111   LogTrace("LHECOMWeightProducer") << "*******LHECOMWeightProducer*******\n"
0112                                    << " Q  : " << Q << "\n"
0113                                    << " id1: " << id1 << "\n"
0114                                    << " x1 : " << x1 << "\n"
0115                                    << " x1': " << x1prime << "\n"
0116                                    << " id2: " << id2 << "\n"
0117                                    << " x2 : " << x2 << "\n"
0118                                    << " x2': " << x2prime;
0119   //gluon is 0 in the LHAPDF numbering
0120   if (id1 == 21)
0121     id1 = 0;
0122   if (id2 == 21)
0123     id2 = 0;
0124
0125   // Put PDF weights in the event
0126   if (verbose)
0127     cout << " Set : " << _pdfset << "  member : " << _pdfmember << endl;
0128
0129   LHAPDF::usePDFMember(1, _pdfmember);
0130   double oldpdf1 = LHAPDF::xfx(1, x1, Q, id1) / x1;
0131   double oldpdf2 = LHAPDF::xfx(1, x2, Q, id2) / x2;
0132   double newpdf1 = LHAPDF::xfx(1, x1prime, Q, id1) / x1prime;
0133   double newpdf2 = LHAPDF::xfx(1, x2prime, Q, id2) / x2prime;
0134   LogTrace("LHECOMWeightProducer") << "     xfx1 : " << oldpdf1 << "\n"
0135                                    << "     xfx2 : " << oldpdf2 << "\n"
0136                                    << "     xfx1': " << newpdf1 << "\n"
0137                                    << "     xfx2': " << newpdf2 << "\n"
0138                                    << "     weight:" << (newpdf1 / oldpdf1) * (newpdf2 / oldpdf2);
0139   double weight = (newpdf1 / oldpdf1) * (newpdf2 / oldpdf2);
0140   std::vector<double> weights;
0141   weights.push_back(weight);
0142   std::unique_ptr<GenEventInfoProduct> info(new GenEventInfoProduct());
0143   info->setWeights(weights);
0144   iEvent.put(std::move(info), _label);
0145 }
0146
0147 DEFINE_FWK_MODULE(LHECOMWeightProducer);