File indexing completed on 2024-04-06 12:13:50
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
0019
0020
0021
0022
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 }
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
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
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
0120 if (id1 == 21)
0121 id1 = 0;
0122 if (id2 == 21)
0123 id2 = 0;
0124
0125
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);