NanoAOD/plugins/NPUTablesProducer.cc

0001 #include "FWCore/Framework/interface/global/EDProducer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0006 #include "DataFormats/NanoAOD/interface/FlatTable.h"
0007 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0008 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0009 #include "DataFormats/VertexReco/interface/Vertex.h"
0010
0011 #include <vector>
0012 #include <iostream>
0013 #include <algorithm>
0014
0015 class NPUTablesProducer : public edm::global::EDProducer<> {
0016 public:
0017   NPUTablesProducer(edm::ParameterSet const& params)
0018       : npuTag_(consumes<std::vector<PileupSummaryInfo>>(params.getParameter<edm::InputTag>("src"))),
0019         pvTag_(consumes<std::vector<reco::Vertex>>(params.getParameter<edm::InputTag>("pvsrc"))),
0020         vz_(params.getParameter<std::vector<double>>("zbins")),
0021         savePtHatMax_(params.getParameter<bool>("savePtHatMax")) {
0022     produces<nanoaod::FlatTable>();
0023   }
0024
0025   ~NPUTablesProducer() override {}
0026
0027   void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
0028     auto npuTab = std::make_unique<nanoaod::FlatTable>(1, "Pileup", true);
0029
0030     const auto& pvProd = iEvent.get(pvTag_);
0031     const double refpvz = pvProd.at(0).position().z();
0032
0033     edm::Handle<std::vector<PileupSummaryInfo>> npuInfo;
0034     if (iEvent.getByToken(npuTag_, npuInfo)) {
0035       fillNPUObjectTable(*npuInfo, *npuTab, refpvz);
0036     }
0037
0038     iEvent.put(std::move(npuTab));
0039   }
0040
0041   void fillNPUObjectTable(const std::vector<PileupSummaryInfo>& npuProd, nanoaod::FlatTable& out, double refpvz) const {
0042     // Get BX 0
0043     unsigned int bx0 = 0;
0044     float nt = 0;
0045     unsigned int npu = 0;
0046
0047     auto zbin = std::lower_bound(vz_.begin(), vz_.end() - 1, std::abs(refpvz));
0048     float pudensity = 0;
0049     float gpudensity = 0;
0050
0051     float pthatmax = 0;
0052
0053     for (unsigned int ibx = 0; ibx < npuProd.size(); ibx++) {
0054       if (npuProd[ibx].getBunchCrossing() == 0) {
0055         bx0 = ibx;
0056         nt = npuProd[ibx].getTrueNumInteractions();
0057         npu = npuProd[ibx].getPU_NumInteractions();
0058
0059         std::vector<float> zpositions;
0060         unsigned int nzpositions = npuProd[ibx].getPU_zpositions().size();
0061         for (unsigned int j = 0; j < nzpositions; ++j) {
0062           zpositions.push_back(npuProd[ibx].getPU_zpositions()[j]);
0063           if (std::abs(zpositions.back() - refpvz) < 0.1)
0064             pudensity++;  //N_PU/mm
0065           auto bin = std::lower_bound(vz_.begin(), vz_.end() - 1, std::abs(zpositions.back()));
0066           if (bin != vz_.end() && bin == zbin)
0067             gpudensity++;
0068         }
0069         gpudensity /= (20.0 * (*(zbin) - *(zbin - 1)));
0070
0071         if (savePtHatMax_) {
0072           if (!npuProd[ibx].getPU_pT_hats().empty()) {
0073             pthatmax = *max_element(npuProd[ibx].getPU_pT_hats().begin(), npuProd[ibx].getPU_pT_hats().end());
0074           }
0075         }
0076       }
0077     }
0078     unsigned int eoot = 0;
0079     for (unsigned int ipu = 0; ipu < bx0; ipu++) {
0080       eoot += npuProd[ipu].getPU_NumInteractions();
0081     }
0082     unsigned int loot = 0;
0083     for (unsigned int ipu = npuProd.size() - 1; ipu > bx0; ipu--) {
0084       loot += npuProd[ipu].getPU_NumInteractions();
0085     }
0086     out.addColumnValue<float>("nTrueInt",
0087                               nt,
0088                               "the true mean number of the poisson distribution for this event from which the number "
0089                               "of interactions each bunch crossing has been sampled",
0090                               10);
0091     out.addColumnValue<int>(
0092         "nPU",
0093         npu,
0094         "the number of pileup interactions that have been added to the event in the current bunch crossing");
0095     out.addColumnValue<int>("sumEOOT", eoot, "number of early out of time pileup");
0096     out.addColumnValue<int>("sumLOOT", loot, "number of late out of time pileup");
0097     out.addColumnValue<float>("pudensity", pudensity, "PU vertices / mm");
0098     out.addColumnValue<float>("gpudensity", gpudensity, "Generator-level PU vertices / mm");
0099     if (savePtHatMax_) {
0100       out.addColumnValue<float>("pthatmax", pthatmax, "Maximum pt-hat", 10);
0101     }
0102   }
0103
0104   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0105     edm::ParameterSetDescription desc;
0106     desc.add<edm::InputTag>("src", edm::InputTag("slimmedAddPileupInfo"))
0107         ->setComment("tag for the PU information (vector<PileupSummaryInfo>)");
0108     desc.add<edm::InputTag>("pvsrc", edm::InputTag("offlineSlimmedPrimaryVertices"))->setComment("tag for the PVs");
0109     desc.add<std::vector<double>>("zbins", {})
0110         ->setComment("Z bins to compute the generator-level number of PU vertices per mm");
0111     desc.add<bool>("savePtHatMax", false)->setComment("Store maximum pt-hat of PU");
0112     descriptions.add("puTable", desc);
0113   }
0114
0115 protected:
0116   const edm::EDGetTokenT<std::vector<PileupSummaryInfo>> npuTag_;
0117   const edm::EDGetTokenT<std::vector<reco::Vertex>> pvTag_;
0118
0119   const std::vector<double> vz_;
0120
0121   bool savePtHatMax_;
0122 };
0123
0124 #include "FWCore/Framework/interface/MakerMacros.h"
0125 DEFINE_FWK_MODULE(NPUTablesProducer);