Merge pull request #28145 from cms-nanoAOD/pr_nanov6_102X

NanoAODv6 updates [102X]

PhysicsTools/NanoAOD/plugins/EventStringOutputBranches.cc

@@ -0,0 +1,37 @@
+#include "PhysicsTools/NanoAOD/plugins/EventStringOutputBranches.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/Registry.h"
+
+#include <iostream>
+
+void 
+EventStringOutputBranches::updateEventStringNames(TTree & tree, const std::string & evstring) 
+{  
+  bool found = false;
+  for (auto & existing : m_evStringBranches) {
+    existing.buffer = false;
+    if (evstring==existing.name) {existing.buffer = true; found=true;}
+  }
+  if (!found && (!evstring.empty())){
+    NamedBranchPtr nb(evstring,"EventString bit");
+    bool backFillValue = false;
+    nb.branch = tree.Branch(nb.name.c_str(), &backFillValue, (nb.name + "/O").c_str());
+    nb.branch->SetTitle(nb.title.c_str());
+    for(size_t i=0;i<m_fills;i++) nb.branch->Fill(); // Back fill
+    nb.buffer = true;
+    m_evStringBranches.push_back(nb);
+    for (auto & existing : m_evStringBranches) existing.branch->SetAddress(&(existing.buffer)); // m_evStringBranches might have been resized
+  }
+}
+
+void EventStringOutputBranches::fill(const edm::EventForOutput &iEvent,TTree & tree) 
+{
+  if ((!m_update_only_at_new_lumi) || m_lastLumi!=iEvent.id().luminosityBlock()) {
+    edm::Handle<std::string> handle;
+    iEvent.getByToken(m_token, handle);
+    const std::string & evstring = *handle;
+    m_lastLumi=iEvent.id().luminosityBlock();
+    updateEventStringNames(tree,evstring);
+  }
+ m_fills++;
+}

PhysicsTools/NanoAOD/plugins/EventStringOutputBranches.h

@@ -0,0 +1,40 @@
+#ifndef PhysicsTools_NanoAOD_EventStringOutputBranches_h
+#define PhysicsTools_NanoAOD_EventStringOutputBranches_h
+
+#include <string>
+#include <vector>
+#include <TTree.h>
+#include "FWCore/Framework/interface/EventForOutput.h"
+#include "DataFormats/Provenance/interface/BranchDescription.h"
+#include "FWCore/Utilities/interface/EDGetToken.h"
+
+class EventStringOutputBranches {
+ public:
+ EventStringOutputBranches(const edm::BranchDescription *desc, const edm::EDGetToken & token, bool update_only_at_new_lumi = false ) :
+  m_token(token), m_lastLumi(-1), m_fills(0), m_update_only_at_new_lumi(update_only_at_new_lumi)
+    {
+        if (desc->className() != "std::basic_string<char,std::char_traits<char> >") throw cms::Exception("Configuration", "NanoAODOutputModule/EventStringOutputBranches can only write out std::string objects");
+    }
+
+	void updateEventStringNames(TTree &, const std::string &);
+	void fill(const edm::EventForOutput &iEvent,TTree & tree) ;
+
+ private:
+
+    edm::EDGetToken m_token;
+    struct NamedBranchPtr {
+        std::string name, title;
+        TBranch * branch;
+	bool buffer;
+        NamedBranchPtr(const std::string & aname, const std::string & atitle, TBranch *branchptr = nullptr) : 
+            name(aname), title(atitle), branch(branchptr), buffer(false) {}
+    };
+    std::vector<NamedBranchPtr> m_evStringBranches;
+    long m_lastLumi;
+    unsigned long m_fills;
+    bool m_update_only_at_new_lumi;
+
+};
+
+#endif
+

PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc

@@ -10,6 +10,7 @@
 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
+#include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoHeader.h"
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "boost/algorithm/string.hpp"
 
@@ -85,6 +86,21 @@ namespace {
         }
     };
 
+  struct CounterMap {
+    std::map<std::string, Counter> countermap;
+    Counter* active_el = nullptr;
+    std::string active_label = "";
+    void merge(const CounterMap & other) {
+      for (const auto &y : other.countermap) countermap[y.first].merge(y.second);
+      active_el = nullptr;
+    }
+    void clear() {for (auto x : countermap) x.second.clear();}
+    void setLabel(std::string label) { active_el = &(countermap[label]); active_label = label;}
+    void checkLabelSet() { if (!active_el) throw cms::Exception("LogicError", "Called CounterMap::get() before setting the active label\n"); }
+    Counter* get() { checkLabelSet(); return active_el; }
+    std::string& getLabel() { checkLabelSet(); return active_label; }
+  };
+
     ///  ---- RunCache object for dynamic choice of LHE IDs ----
     struct DynamicWeightChoice {
         // choice of LHE weights
@@ -138,13 +154,14 @@ namespace {
     };
 }
 
-class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<Counter>, edm::RunCache<DynamicWeightChoice>, edm::RunSummaryCache<Counter>, edm::EndRunProducer> {
+class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<CounterMap>, edm::RunCache<DynamicWeightChoice>, edm::RunSummaryCache<CounterMap>, edm::EndRunProducer> {
     public:
         GenWeightsTableProducer( edm::ParameterSet const & params ) :
             genTag_(consumes<GenEventInfoProduct>(params.getParameter<edm::InputTag>("genEvent"))),
             lheLabel_(params.getParameter<std::vector<edm::InputTag>>("lheInfo")),
             lheTag_(edm::vector_transform(lheLabel_, [this](const edm::InputTag & tag) { return mayConsume<LHEEventProduct>(tag); })),
             lheRunTag_(edm::vector_transform(lheLabel_, [this](const edm::InputTag & tag) { return mayConsume<LHERunInfoProduct, edm::InRun>(tag); })),
+            genLumiInfoHeadTag_(mayConsume<GenLumiInfoHeader,edm::InLumi>(params.getParameter<edm::InputTag>("genLumiInfoHeader"))),
             namedWeightIDs_(params.getParameter<std::vector<std::string>>("namedWeightIDs")),
             namedWeightLabels_(params.getParameter<std::vector<std::string>>("namedWeightLabels")),
             lheWeightPrecision_(params.getParameter<int32_t>("lheWeightPrecision")),
@@ -153,6 +170,7 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
             hasIssuedWarning_(false)
         {
             produces<nanoaod::FlatTable>();
+            produces<std::string>("genModel");
             produces<nanoaod::FlatTable>("LHEScale");
             produces<nanoaod::FlatTable>("LHEPdf");
             produces<nanoaod::FlatTable>("LHEReweighting");
@@ -175,7 +193,7 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
 
         void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
             // get my counter for weights
-            Counter * counter = streamCache(id);
+            Counter * counter = streamCache(id)->get();
 
             // generator information (always available)
             edm::Handle<GenEventInfoProduct> genInfo;
@@ -188,6 +206,10 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
             out->addColumnValue<float>("", weight, "generator weight", nanoaod::FlatTable::FloatColumn);
             iEvent.put(std::move(out));
 
+	    std::string model_label = streamCache(id)->getLabel();
+	    auto outM = std::make_unique<std::string>((!model_label.empty()) ? std::string("GenModel_") + model_label : "");
+	    iEvent.put(std::move(outM),"genModel");
+
             // tables for LHE weights, may not be filled
             std::unique_ptr<nanoaod::FlatTable> lheScaleTab, lhePdfTab, lheRwgtTab, lheNamedTab;
             std::unique_ptr<nanoaod::FlatTable> genPSTab;
@@ -565,45 +587,61 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
 
 
         // create an empty counter
-        std::unique_ptr<Counter> beginStream(edm::StreamID) const override { 
-            return std::make_unique<Counter>(); 
+        std::unique_ptr<CounterMap> beginStream(edm::StreamID) const override {
+            return std::make_unique<CounterMap>();
         }
         // inizialize to zero at begin run
         void streamBeginRun(edm::StreamID id, edm::Run const&, edm::EventSetup const&) const override { 
             streamCache(id)->clear(); 
         }
+        void streamBeginLuminosityBlock(edm::StreamID id, edm::LuminosityBlock const& lumiBlock, edm::EventSetup const& eventSetup) const override {
+	  auto counterMap = streamCache(id);
+	  edm::Handle<GenLumiInfoHeader> genLumiInfoHead;
+	  lumiBlock.getByToken(genLumiInfoHeadTag_,genLumiInfoHead);
+	  if (!genLumiInfoHead.isValid()) edm::LogWarning("LHETablesProducer") << "No GenLumiInfoHeader product found, will not fill generator model string.\n";
+	  counterMap->setLabel(genLumiInfoHead.isValid() ? genLumiInfoHead->configDescription() : "");
+	}
         // create an empty counter
-        std::shared_ptr<Counter> globalBeginRunSummary(edm::Run const&, edm::EventSetup const&) const override { 
-            return std::make_shared<Counter>(); 
+        std::shared_ptr<CounterMap> globalBeginRunSummary(edm::Run const&, edm::EventSetup const&) const override {
+            return std::make_shared<CounterMap>();
         }
         // add this stream to the summary
-        void streamEndRunSummary(edm::StreamID id, edm::Run const&, edm::EventSetup const&, Counter* runCounter) const override { 
-            runCounter->merge(*streamCache(id)); 
+        void streamEndRunSummary(edm::StreamID id, edm::Run const&, edm::EventSetup const&, CounterMap* runCounterMap) const override {
+            runCounterMap->merge(*streamCache(id));
         }
         // nothing to do per se
-        void globalEndRunSummary(edm::Run const&, edm::EventSetup const&, Counter* runCounter) const override { 
+        void globalEndRunSummary(edm::Run const&, edm::EventSetup const&, CounterMap* runCounterMap) const override {
         }
         // write the total to the run 
-        void globalEndRunProduce(edm::Run& iRun, edm::EventSetup const&, Counter const* runCounter) const override {
+        void globalEndRunProduce(edm::Run& iRun, edm::EventSetup const&, CounterMap const* runCounterMap) const override {
             auto out = std::make_unique<nanoaod::MergeableCounterTable>();
-            out->addInt("genEventCount", "event count", runCounter->num);
-            out->addFloat("genEventSumw", "sum of gen weights", runCounter->sumw);
-            out->addFloat("genEventSumw2", "sum of gen (weight^2)", runCounter->sumw2);
+
+	    for (auto x : runCounterMap->countermap) {
+
+	      auto runCounter = &(x.second);
+	      std::string label = std::string("_") + x.first;
+	      std::string doclabel = (!x.first.empty()) ? (std::string(", for model label ") + x.first) : "";
+
+            out->addInt("genEventCount"+label, "event count"+doclabel, runCounter->num);
+            out->addFloat("genEventSumw"+label, "sum of gen weights"+doclabel, runCounter->sumw);
+            out->addFloat("genEventSumw2"+label, "sum of gen (weight^2)"+doclabel, runCounter->sumw2);
 
             double norm = runCounter->sumw ? 1.0/runCounter->sumw : 1;
             auto sumScales = runCounter->sumScale; for (auto & val : sumScales) val *= norm;
-            out->addVFloat("LHEScaleSumw", "Sum of genEventWeight * LHEScaleWeight[i], divided by genEventSumw", sumScales);
+            out->addVFloat("LHEScaleSumw"+label, "Sum of genEventWeight * LHEScaleWeight[i], divided by genEventSumw"+doclabel, sumScales);
             auto sumPDFs = runCounter->sumPDF; for (auto & val : sumPDFs) val *= norm;
-            out->addVFloat("LHEPdfSumw", "Sum of genEventWeight * LHEPdfWeight[i], divided by genEventSumw", sumPDFs);
+            out->addVFloat("LHEPdfSumw"+label, "Sum of genEventWeight * LHEPdfWeight[i], divided by genEventSumw"+doclabel, sumPDFs);
             if (!runCounter->sumRwgt.empty()) {
                 auto sumRwgts = runCounter->sumRwgt; for (auto & val : sumRwgts) val *= norm;
-                out->addVFloat("LHEReweightingSumw", "Sum of genEventWeight * LHEReweightingWeight[i], divided by genEventSumw", sumRwgts);
+                out->addVFloat("LHEReweightingSumw"+label, "Sum of genEventWeight * LHEReweightingWeight[i], divided by genEventSumw"+doclabel, sumRwgts);
             }
             if (!runCounter->sumNamed.empty()) { // it could be empty if there's no LHE info in the sample
                 for (unsigned int i = 0, n = namedWeightLabels_.size(); i < n; ++i) {
-                    out->addFloat("LHESumw_"+namedWeightLabels_[i], "Sum of genEventWeight * LHEWeight_"+namedWeightLabels_[i]+", divided by genEventSumw", runCounter->sumNamed[i] * norm);
+                    out->addFloat("LHESumw_"+namedWeightLabels_[i]+label, "Sum of genEventWeight * LHEWeight_"+namedWeightLabels_[i]+", divided by genEventSumw"+doclabel, runCounter->sumNamed[i] * norm);
                 }
             }
+
+	    }
             iRun.put(std::move(out));
         }
         // nothing to do here
@@ -612,6 +650,7 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
         static void fillDescriptions(edm::ConfigurationDescriptions & descriptions) {
             edm::ParameterSetDescription desc;
             desc.add<edm::InputTag>("genEvent", edm::InputTag("generator"))->setComment("tag for the GenEventInfoProduct, to get the main weight");
+            desc.add<edm::InputTag>("genLumiInfoHeader", edm::InputTag("generator"))->setComment("tag for the GenLumiInfoProduct, to get the model string");
             desc.add<std::vector<edm::InputTag>>("lheInfo", std::vector<edm::InputTag>{{"externalLHEProducer"},{"source"}})->setComment("tag(s) for the LHE information (LHEEventProduct and LHERunInfoProduct)");
 
             edm::ParameterSetDescription prefpdf;
@@ -632,6 +671,7 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
         const std::vector<edm::InputTag> lheLabel_;
         const std::vector<edm::EDGetTokenT<LHEEventProduct>> lheTag_;
         const std::vector<edm::EDGetTokenT<LHERunInfoProduct>> lheRunTag_;
+        const edm::EDGetTokenT<GenLumiInfoHeader> genLumiInfoHeadTag_;
 
         std::vector<uint32_t> preferredPDFLHAIDs_;
         std::unordered_map<std::string,uint32_t> lhaNameToID_;

PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc

@@ -38,6 +38,7 @@
 #include "DataFormats/NanoAOD/interface/UniqueString.h"
 #include "PhysicsTools/NanoAOD/plugins/TableOutputBranches.h"
 #include "PhysicsTools/NanoAOD/plugins/TriggerOutputBranches.h"
+#include "PhysicsTools/NanoAOD/plugins/EventStringOutputBranches.h"
 #include "PhysicsTools/NanoAOD/plugins/SummaryTableOutputBranches.h"
 
 #include <iostream>
@@ -115,6 +116,7 @@ class NanoAODOutputModule : public edm::one::OutputModule<> {
 
   std::vector<TableOutputBranches> m_tables;
   std::vector<TriggerOutputBranches> m_triggers;
+  std::vector<EventStringOutputBranches> m_evstrings;
 
   std::vector<SummaryTableOutputBranches> m_runTables;
 
@@ -198,6 +200,8 @@ NanoAODOutputModule::write(edm::EventForOutput const& iEvent) {
   }
   // fill triggers
   for (auto & t : m_triggers) t.fill(iEvent,*m_tree);
+  // fill event branches
+  for (auto & t : m_evstrings) t.fill(iEvent,*m_tree);
   m_tree->Fill();
 
   m_processHistoryRegistry.registerProcessHistory(iEvent.processHistory());
@@ -272,6 +276,7 @@ NanoAODOutputModule::openFile(edm::FileBlock const&) {
   /* Setup file structure here */
   m_tables.clear();
   m_triggers.clear();
+  m_evstrings.clear();
   m_runTables.clear();
   const auto & keeps = keptProducts();
   for (const auto & keep : keeps[edm::InEvent]) {
@@ -281,6 +286,9 @@ NanoAODOutputModule::openFile(edm::FileBlock const&) {
 	  {
 	      m_triggers.emplace_back(keep.first, keep.second);
 	  }
+      else if(keep.first->className() == "std::basic_string<char,std::char_traits<char> >" && keep.first->productInstanceName()=="genModel") { // friendlyClassName == "String"
+	m_evstrings.emplace_back(keep.first, keep.second, true); // update only at lumiBlock transitions
+      }
       else throw cms::Exception("Configuration", "NanoAODOutputModule cannot handle class " + keep.first->className());     
   }
 
@@ -358,7 +366,7 @@ NanoAODOutputModule::fillDescriptions(edm::ConfigurationDescriptions& descriptio
         ->setComment("Autoflush parameter for ROOT file");
 
   //replace with whatever you want to get from the EDM by default
-  const std::vector<std::string> keep = {"drop *", "keep nanoaodFlatTable_*Table_*_*", "keep edmTriggerResults_*_*_*", "keep nanoaodMergeableCounterTable_*Table_*_*", "keep nanoaodUniqueString_nanoMetadata_*_*"};
+  const std::vector<std::string> keep = {"drop *", "keep nanoaodFlatTable_*Table_*_*", "keep edmTriggerResults_*_*_*", "keep String_*_genModel_*", "keep nanoaodMergeableCounterTable_*Table_*_*", "keep nanoaodUniqueString_nanoMetadata_*_*"};
   edm::OutputModule::fillDescription(desc, keep);
 
   //Used by Workflow management for their own meta data

PhysicsTools/NanoAOD/plugins/VertexTableProducer.cc

@@ -37,6 +37,7 @@
 
 #include "DataFormats/NanoAOD/interface/FlatTable.h"
 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
+#include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
 #include "DataFormats/Common/interface/ValueMap.h"
@@ -149,8 +150,9 @@ VertexTableProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
     edm::Handle<edm::View<reco::VertexCompositePtrCandidate> > svsIn;
     iEvent.getByToken(svs_, svsIn);
     auto selCandSv = std::make_unique<PtrVector<reco::Candidate>>();
-    std::vector<float> dlen,dlenSig,pAngle;
+    std::vector<float> dlen,dlenSig,pAngle,dxy,dxySig;
     VertexDistance3D vdist;
+    VertexDistanceXY vdistXY;
 
     size_t i=0;
     const auto & PV0 = pvsIn->front();
@@ -163,8 +165,11 @@ VertexTableProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
                edm::Ptr<reco::Candidate> c =  svsIn->ptrAt(i);
                selCandSv->push_back(c);
                double dx = (PV0.x() - sv.vx()), dy = (PV0.y() - sv.vy()), dz = (PV0.z() - sv.vz());
-               double pdotv = (dx * sv.px() + dy*sv.py() + dz*sv.pz())/sv.p();
+               double pdotv = (dx * sv.px() + dy*sv.py() + dz*sv.pz())/sv.p()/sqrt(dx*dx + dy*dy + dz*dz);
                pAngle.push_back(std::acos(pdotv));
+	       Measurement1D d2d = vdistXY.distance(PV0, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
+	       dxy.push_back(d2d.value());
+	       dxySig.push_back(d2d.significance());
            }
        }
        i++;
@@ -175,6 +180,8 @@ VertexTableProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
     // For SV we fill from here only stuff that cannot be created with the SimpleFlatTableProducer 
     svsTable->addColumn<float>("dlen",dlen,"decay length in cm",nanoaod::FlatTable::FloatColumn,10);
     svsTable->addColumn<float>("dlenSig",dlenSig,"decay length significance",nanoaod::FlatTable::FloatColumn, 10);
+    svsTable->addColumn<float>("dxy", dxy, "2D decay length in cm", nanoaod::FlatTable::FloatColumn, 10);
+    svsTable->addColumn<float>("dxySig", dxySig, "2D decay length significance", nanoaod::FlatTable::FloatColumn, 10);
     svsTable->addColumn<float>("pAngle",pAngle,"pointing angle, i.e. acos(p_SV * (SV - PV)) ",nanoaod::FlatTable::FloatColumn,10);
 
 

PhysicsTools/NanoAOD/python/NanoAODEDMEventContent_cff.py

@@ -5,6 +5,7 @@
         'drop *',
         "keep nanoaodFlatTable_*Table_*_*",     # event data
         "keep edmTriggerResults_*_*_*",  # event data
+        "keep String_*_genModel_*",  # generator model data
         "keep nanoaodMergeableCounterTable_*Table_*_*", # accumulated per/run or per/lumi data
         "keep nanoaodUniqueString_nanoMetadata_*_*",   # basic metadata
     )