Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-19 05:08:37

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     edm::Handle<std::vector<reco::Vertex>> pvsIn;
0031     iEvent.getByToken(pvTag_, pvsIn);
0032     const double refpvz = (*pvsIn)[0].position().z();
0033 
0034     edm::Handle<std::vector<PileupSummaryInfo>> npuInfo;
0035     if (iEvent.getByToken(npuTag_, npuInfo)) {
0036       fillNPUObjectTable(*npuInfo, *npuTab, refpvz);
0037     }
0038 
0039     iEvent.put(std::move(npuTab));
0040   }
0041 
0042   void fillNPUObjectTable(const std::vector<PileupSummaryInfo>& npuProd, nanoaod::FlatTable& out, double refpvz) const {
0043     // Get BX 0
0044     unsigned int bx0 = 0;
0045     float nt = 0;
0046     unsigned int npu = 0;
0047 
0048     auto zbin = std::lower_bound(vz_.begin(), vz_.end() - 1, std::abs(refpvz));
0049     float pudensity = 0;
0050     float gpudensity = 0;
0051 
0052     float pthatmax = 0;
0053 
0054     for (unsigned int ibx = 0; ibx < npuProd.size(); ibx++) {
0055       if (npuProd[ibx].getBunchCrossing() == 0) {
0056         bx0 = ibx;
0057         nt = npuProd[ibx].getTrueNumInteractions();
0058         npu = npuProd[ibx].getPU_NumInteractions();
0059 
0060         std::vector<float> zpositions;
0061         unsigned int nzpositions = npuProd[ibx].getPU_zpositions().size();
0062         for (unsigned int j = 0; j < nzpositions; ++j) {
0063           zpositions.push_back(npuProd[ibx].getPU_zpositions()[j]);
0064           if (std::abs(zpositions.back() - refpvz) < 0.1)
0065             pudensity++;  //N_PU/mm
0066           auto bin = std::lower_bound(vz_.begin(), vz_.end() - 1, std::abs(zpositions.back()));
0067           if (bin != vz_.end() && bin == zbin)
0068             gpudensity++;
0069         }
0070         gpudensity /= (20.0 * (*(zbin) - *(zbin - 1)));
0071 
0072         if (savePtHatMax_) {
0073           if (!npuProd[ibx].getPU_pT_hats().empty()) {
0074             pthatmax = *max_element(npuProd[ibx].getPU_pT_hats().begin(), npuProd[ibx].getPU_pT_hats().end());
0075           }
0076         }
0077       }
0078     }
0079     unsigned int eoot = 0;
0080     for (unsigned int ipu = 0; ipu < bx0; ipu++) {
0081       eoot += npuProd[ipu].getPU_NumInteractions();
0082     }
0083     unsigned int loot = 0;
0084     for (unsigned int ipu = npuProd.size() - 1; ipu > bx0; ipu--) {
0085       loot += npuProd[ipu].getPU_NumInteractions();
0086     }
0087     out.addColumnValue<float>("nTrueInt",
0088                               nt,
0089                               "the true mean number of the poisson distribution for this event from which the number "
0090                               "of interactions each bunch crossing has been sampled",
0091                               10);
0092     out.addColumnValue<int>(
0093         "nPU",
0094         npu,
0095         "the number of pileup interactions that have been added to the event in the current bunch crossing");
0096     out.addColumnValue<int>("sumEOOT", eoot, "number of early out of time pileup");
0097     out.addColumnValue<int>("sumLOOT", loot, "number of late out of time pileup");
0098     out.addColumnValue<float>("pudensity", pudensity, "PU vertices / mm");
0099     out.addColumnValue<float>("gpudensity", gpudensity, "Generator-level PU vertices / mm");
0100     if (savePtHatMax_) {
0101       out.addColumnValue<float>("pthatmax", pthatmax, "Maximum pt-hat");
0102     }
0103   }
0104 
0105   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0106     edm::ParameterSetDescription desc;
0107     desc.add<edm::InputTag>("src", edm::InputTag("slimmedAddPileupInfo"))
0108         ->setComment("tag for the PU information (vector<PileupSummaryInfo>)");
0109     desc.add<edm::InputTag>("pvsrc", edm::InputTag("offlineSlimmedPrimaryVertices"))->setComment("tag for the PVs");
0110     desc.add<std::vector<double>>("zbins", {})
0111         ->setComment("Z bins to compute the generator-level number of PU vertices per mm");
0112     desc.add<bool>("savePtHatMax", false)->setComment("Store maximum pt-hat of PU");
0113     descriptions.add("puTable", desc);
0114   }
0115 
0116 protected:
0117   const edm::EDGetTokenT<std::vector<PileupSummaryInfo>> npuTag_;
0118   const edm::EDGetTokenT<std::vector<reco::Vertex>> pvTag_;
0119 
0120   const std::vector<double> vz_;
0121 
0122   bool savePtHatMax_;
0123 };
0124 
0125 #include "FWCore/Framework/interface/MakerMacros.h"
0126 DEFINE_FWK_MODULE(NPUTablesProducer);