Merge pull request #13019 from bendavid/pdfweights_80x

Add PDF weight helper tool (80x, updated)

PhysicsTools/HepMCCandAlgos/BuildFile.xml

@@ -2,7 +2,10 @@
 <use   name="DataFormats/Candidate"/>
 <use   name="DataFormats/Common"/>
 <use   name="DataFormats/HepMCCandidate"/>
+<use   name="FWCore/ServiceRegistry"/>
+<use   name="CommonTools/UtilAlgos"/>
 <use   name="hepmc"/>
+<use   name="eigen"/>
 <export>
   <lib   name="1"/>
 </export>

PhysicsTools/HepMCCandAlgos/interface/PDFWeightsHelper.h

@@ -0,0 +1,26 @@
+#ifndef PhysicsTools_HepMCCandAlgos_PDFWeightsHelper_h
+#define PhysicsTools_HepMCCandAlgos_PDFWeightsHelper_h
+
+#include <Eigen/Dense>
+
+#include <iostream>
+
+#include "FWCore/ParameterSet/interface/FileInPath.h"
+
+class PDFWeightsHelper {
+
+public:
+
+  PDFWeightsHelper();
+
+  void Init(unsigned int nreplicas, unsigned int neigenvectors, const edm::FileInPath &incsv);
+  void DoMC2Hessian(double nomweight, const double *inweights, double *outweights) const;
+
+  unsigned int neigenvectors() const { return transformation_.cols(); }
+
+protected:
+
+  Eigen::MatrixXd transformation_;  
+
+};
+#endif

PhysicsTools/HepMCCandAlgos/plugins/PDFWeightsTest.cc

@@ -0,0 +1,114 @@
+#include "FWCore/Framework/interface/one/EDAnalyzer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "DataFormats/Common/interface/Handle.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
+#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
+#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
+#include "PhysicsTools/HepMCCandAlgos/interface/PDFWeightsHelper.h"
+#include "FWCore/Utilities/interface/EDMException.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "CommonTools/UtilAlgos/interface/TFileService.h"
+#include "TTree.h"
+#include <iostream>
+using namespace std;
+using namespace edm;
+using namespace HepMC;
+
+class PDFWeightsTest : public edm::one::EDAnalyzer<> {
+private:
+  PDFWeightsHelper pdfweightshelper_;
+  EDGetTokenT<LHEEventProduct> srcToken_;
+  EDGetTokenT<LHEEventProduct> srcTokenAlt_;
+  EDGetTokenT<GenEventInfoProduct> srcTokenGen_;
+
+  unsigned int pdfWeightOffset_;
+  unsigned int nPdfWeights_;
+  unsigned int nPdfEigWeights_;
+
+  TTree *tree_;
+  std::vector<float> pdfweights_;
+  std::vector<float> pdfeigweights_;
+  float weight_;
+
+public:
+  explicit PDFWeightsTest( const ParameterSet & cfg ) :
+    srcToken_( consumes<LHEEventProduct>(edm::InputTag("externalLHEProducer"))),
+    srcTokenAlt_( consumes<LHEEventProduct>(edm::InputTag("source"))),
+    srcTokenGen_( consumes<GenEventInfoProduct>(edm::InputTag("generator"))),
+    pdfWeightOffset_(cfg.getParameter<unsigned int>("pdfWeightOffset")),
+    nPdfWeights_(cfg.getParameter<unsigned int>("nPdfWeights")),
+    nPdfEigWeights_(cfg.getParameter<unsigned int>("nPdfEigWeights")),
+    pdfweights_(nPdfWeights_),
+    pdfeigweights_(nPdfEigWeights_)
+  {
+
+    edm::Service<TFileService> fs;
+    tree_ = fs->make<TTree>("tree", "");
+
+    tree_->Branch("pdfrep",&pdfweights_);
+    tree_->Branch("pdfeig",&pdfeigweights_);
+    tree_->Branch("weight",&weight_);
+
+    edm::FileInPath mc2hessianCSV = cfg.getParameter<edm::FileInPath>("mc2hessianCSV");
+    pdfweightshelper_.Init(nPdfWeights_,nPdfEigWeights_,mc2hessianCSV);
+
+  }
+
+private:
+  void analyze( const Event & evt, const EventSetup& es ) override {
+    Handle<LHEEventProduct> lheInfo;
+    evt.getByToken( srcToken_, lheInfo );
+
+    if (!lheInfo.isValid()) {
+      evt.getByToken( srcTokenAlt_, lheInfo );
+    }
+
+    double nomlheweight = lheInfo->weights()[0].wgt;
+
+    Handle<GenEventInfoProduct> genInfo;
+    evt.getByToken( srcTokenGen_, genInfo );
+
+    weight_ = genInfo->weight();
+
+    //get the original mc replica weights
+    std::vector<double> inpdfweights(nPdfWeights_);
+    for (unsigned int ipdf=0; ipdf<nPdfWeights_; ++ipdf) {
+      unsigned int iwgt = ipdf + pdfWeightOffset_;
+
+      //this is the weight to be used for evaluating uncertainties with mc replica weights
+      pdfweights_[ipdf] = lheInfo->weights()[iwgt].wgt*weight_/nomlheweight;
+
+      //this is the raw weight to be fed to the mc2hessian convertor
+      inpdfweights[ipdf] = lheInfo->weights()[iwgt].wgt;
+
+    }
+
+    std::vector<double> outpdfweights(nPdfEigWeights_);
+    //do the actual conversion, where the nominal lhe weight is needed as the reference point for the linearization
+    pdfweightshelper_.DoMC2Hessian(nomlheweight,inpdfweights.data(),outpdfweights.data());
+
+    for (unsigned int iwgt=0; iwgt<nPdfEigWeights_; ++iwgt) {
+      double wgtval = outpdfweights[iwgt];
+
+      //the is the weight to be used for evaluating uncertainties with hessian weights
+      pdfeigweights_[iwgt] = wgtval*weight_/nomlheweight;
+    }    
+
+    tree_->Fill();
+
+  }
+
+
+};
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+DEFINE_FWK_MODULE( PDFWeightsTest );
+
+
+

PhysicsTools/HepMCCandAlgos/src/PDFWeightsHelper.cc

@@ -0,0 +1,50 @@
+#include "PhysicsTools/HepMCCandAlgos/interface/PDFWeightsHelper.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include <fstream>
+
+PDFWeightsHelper::PDFWeightsHelper() {
+
+}
+
+void PDFWeightsHelper::Init(unsigned int nreplicas, unsigned int neigenvectors, const edm::FileInPath &incsv) {  
+
+  transformation_.resize(nreplicas,neigenvectors);
+
+  std::ifstream instream(incsv.fullPath());
+  if (!instream.is_open()) {
+    throw cms::Exception("PDFWeightsHelper")
+      << "Could not open csv file" << incsv.relativePath(); 
+  }
+
+
+  for (unsigned int ireplica = 0; ireplica<nreplicas; ++ireplica) {
+    std::string linestr;
+    getline(instream,linestr);
+    std::istringstream line(linestr);
+    for (unsigned int ieigen = 0; ieigen<neigenvectors; ++ieigen) {
+      std::string valstr;
+      getline(line,valstr,',');
+      std::istringstream val(valstr);
+      val >> transformation_(ireplica,ieigen);
+    }
+
+  }
+
+}
+
+void PDFWeightsHelper::DoMC2Hessian(double nomweight, const double *inweights, double *outweights) const {
+  const unsigned int nreplicas = transformation_.rows();
+  const unsigned int neigenvectors = transformation_.cols();
+
+  Eigen::VectorXd inweightv(nreplicas);
+  for (unsigned int irep=0; irep<nreplicas; ++irep) {
+    inweightv[irep] = inweights[irep] - nomweight;
+  }
+
+  Eigen::VectorXd outweightv = transformation_.transpose()*inweightv;
+
+  for (unsigned int ieig=0; ieig<neigenvectors; ++ieig) {
+    outweights[ieig] = outweightv[ieig] + nomweight;
+  }  
+
+}

PhysicsTools/HepMCCandAlgos/test/testPDFWeightsTest.py

@@ -0,0 +1,59 @@
+# Auto generated configuration file
+# using: 
+# Revision: 1.19 
+# Source: /local/reps/CMSSW/CMSSW/Configuration/Applications/python/ConfigBuilder.py,v 
+# with command line options: testpdf -s NONE --no_exec --conditions auto:run2_mc -n -1 --filein file:~/work/CMSSWpdfweight/events_orig.lhe
+import FWCore.ParameterSet.Config as cms
+
+process = cms.Process('LHE')
+
+# import of standard configurations
+process.load('FWCore.MessageService.MessageLogger_cfi')
+process.load('Configuration.EventContent.EventContent_cff')
+process.load('SimGeneral.MixingModule.mixNoPU_cfi')
+process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_condDBv2_cff')
+
+process.maxEvents = cms.untracked.PSet(
+    input = cms.untracked.int32(1000)
+)
+
+# Input source
+process.source = cms.Source("PoolSource",
+    fileNames = cms.untracked.vstring('/store/mc/RunIISpring15FSPremix/SMS-T1bbbb_mGluino-1150_mLSP-400to975-1100to1125_TuneCUETP8M1_13TeV-madgraphMLM-pythia8/MINIAODSIM/MCRUN2_74_V9-v1/50000/320B2CE3-BE5D-E511-8C9E-B083FED76C6C.root')
+)
+
+process.options = cms.untracked.PSet(
+
+)
+
+# Production Info
+process.configurationMetadata = cms.untracked.PSet(
+    annotation = cms.untracked.string('testpdf nevts:-1'),
+    name = cms.untracked.string('Applications'),
+    version = cms.untracked.string('$Revision: 1.19 $')
+)
+
+
+#TFileService for output 
+process.TFileService = cms.Service("TFileService", 
+    fileName = cms.string("testpdf.root"),
+    closeFileFast = cms.untracked.bool(True)
+)
+
+# Other statements
+from Configuration.AlCa.GlobalTag_condDBv2 import GlobalTag
+process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_mc', '')
+
+process.testpdf = cms.EDAnalyzer("PDFWeightsTest",
+                                 pdfWeightOffset = cms.uint32(10), #index of first mc replica weight (careful, this should not be the nominal weight, which is repeated in some mc samples).  The majority of run2 LO madgraph_aMC@NLO samples with 5fs matrix element and pdf would use index 10, corresponding to pdf set 263001, the first alternate mc replica for the nominal pdf set 263000 used for these samples
+                                 nPdfWeights = cms.uint32(100), #number of input weights
+                                 nPdfEigWeights = cms.uint32(60), #number of output weights
+                                 mc2hessianCSV = cms.FileInPath('PhysicsTools/HepMCCandAlgos/data/NNPDF30_lo_as_0130_hessian_60.csv'), #MC2Hessian transformation matrix
+                                 )
+
+process.ana = cms.Path(process.testpdf)
+
+# Schedule definition
+process.schedule = cms.Schedule(process.ana)
+
+