HepMCCandAlgos/plugins/PDFWeightsTest.cc

0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/Utilities/interface/InputTag.h"
0007 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0008 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0009 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0010 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0011 #include "PhysicsTools/HepMCCandAlgos/interface/PDFWeightsHelper.h"
0012 #include "FWCore/Utilities/interface/EDMException.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0015 #include "TTree.h"
0016 #include <iostream>
0017 using namespace std;
0018 using namespace edm;
0019 using namespace HepMC;
0020
0021 class PDFWeightsTest : public edm::one::EDAnalyzer<> {
0022 private:
0023   PDFWeightsHelper pdfweightshelper_;
0024   EDGetTokenT<LHEEventProduct> srcToken_;
0025   EDGetTokenT<LHEEventProduct> srcTokenAlt_;
0026   EDGetTokenT<GenEventInfoProduct> srcTokenGen_;
0027
0028   unsigned int pdfWeightOffset_;
0029   unsigned int nPdfWeights_;
0030   unsigned int nPdfEigWeights_;
0031
0032   TTree* tree_;
0033   std::vector<float> pdfweights_;
0034   std::vector<float> pdfeigweights_;
0035   float weight_;
0036
0037 public:
0038   explicit PDFWeightsTest(const ParameterSet& cfg)
0039       : srcToken_(consumes<LHEEventProduct>(edm::InputTag("externalLHEProducer"))),
0040         srcTokenAlt_(consumes<LHEEventProduct>(edm::InputTag("source"))),
0041         srcTokenGen_(consumes<GenEventInfoProduct>(edm::InputTag("generator"))),
0042         pdfWeightOffset_(cfg.getParameter<unsigned int>("pdfWeightOffset")),
0043         nPdfWeights_(cfg.getParameter<unsigned int>("nPdfWeights")),
0044         nPdfEigWeights_(cfg.getParameter<unsigned int>("nPdfEigWeights")),
0045         pdfweights_(nPdfWeights_),
0046         pdfeigweights_(nPdfEigWeights_) {
0047     edm::Service<TFileService> fs;
0048     tree_ = fs->make<TTree>("tree", "");
0049
0050     tree_->Branch("pdfrep", &pdfweights_);
0051     tree_->Branch("pdfeig", &pdfeigweights_);
0052     tree_->Branch("weight", &weight_);
0053
0054     edm::FileInPath mc2hessianCSV = cfg.getParameter<edm::FileInPath>("mc2hessianCSV");
0055     pdfweightshelper_.Init(nPdfWeights_, nPdfEigWeights_, mc2hessianCSV);
0056   }
0057
0058 private:
0059   void analyze(const Event& evt, const EventSetup& es) override {
0060     Handle<LHEEventProduct> lheInfo;
0061     evt.getByToken(srcToken_, lheInfo);
0062
0063     if (!lheInfo.isValid()) {
0064       evt.getByToken(srcTokenAlt_, lheInfo);
0065     }
0066
0067     double nomlheweight = lheInfo->weights()[0].wgt;
0068
0069     Handle<GenEventInfoProduct> genInfo;
0070     evt.getByToken(srcTokenGen_, genInfo);
0071
0072     weight_ = genInfo->weight();
0073
0074     //get the original mc replica weights
0075     std::vector<double> inpdfweights(nPdfWeights_);
0076     for (unsigned int ipdf = 0; ipdf < nPdfWeights_; ++ipdf) {
0077       unsigned int iwgt = ipdf + pdfWeightOffset_;
0078
0079       //this is the weight to be used for evaluating uncertainties with mc replica weights
0080       pdfweights_[ipdf] = lheInfo->weights()[iwgt].wgt * weight_ / nomlheweight;
0081
0082       //this is the raw weight to be fed to the mc2hessian convertor
0083       inpdfweights[ipdf] = lheInfo->weights()[iwgt].wgt;
0084     }
0085
0086     std::vector<double> outpdfweights(nPdfEigWeights_);
0087     //do the actual conversion, where the nominal lhe weight is needed as the reference point for the linearization
0088     pdfweightshelper_.DoMC2Hessian(nomlheweight, inpdfweights.data(), outpdfweights.data());
0089
0090     for (unsigned int iwgt = 0; iwgt < nPdfEigWeights_; ++iwgt) {
0091       double wgtval = outpdfweights[iwgt];
0092
0093       //the is the weight to be used for evaluating uncertainties with hessian weights
0094       pdfeigweights_[iwgt] = wgtval * weight_ / nomlheweight;
0095     }
0096
0097     tree_->Fill();
0098   }
0099 };
0100
0101 #include "FWCore/Framework/interface/MakerMacros.h"
0102
0103 DEFINE_FWK_MODULE(PDFWeightsTest);