JetAnalyzers/src/FFTJetPileupAnalyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    FFTJetPileupAnalyzer
0004 // Class:      FFTJetPileupAnalyzer
0005 //
0006 /**\class FFTJetPileupAnalyzer FFTJetPileupAnalyzer.cc RecoJets/JetAnalyzers/src/FFTJetPileupAnalyzer.cc
0007
0008  Description: collects the info produced by FFTJetPileupProcessor and FFTJetPileupEstimator
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Igor Volobouev
0015 //         Created:  Thu Apr 21 15:52:11 CDT 2011
0016 //
0017 //
0018
0019 #include <cassert>
0020 #include <sstream>
0021 #include <numeric>
0022
0023 #include "TNtuple.h"
0024
0025 // user include files
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 #include "FWCore/ServiceRegistry/interface/Service.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032
0033 #include "DataFormats/Common/interface/Handle.h"
0034 #include <TH2D.h>
0035 #include "DataFormats/JetReco/interface/FFTJetPileupSummary.h"
0036 #include "DataFormats/VertexReco/interface/Vertex.h"
0037 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0038
0039 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0040
0041 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0042
0043 #include "DataFormats/JetReco/interface/DiscretizedEnergyFlow.h"
0044
0045 #include "DataFormats/Provenance/interface/RunLumiEventNumber.h"
0046
0047 #define init_param(type, varname) varname(ps.getParameter<type>(#varname))
0048
0049 //
0050 // class declaration
0051 //
0052 class FFTJetPileupAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0053 public:
0054   explicit FFTJetPileupAnalyzer(const edm::ParameterSet&);
0055   FFTJetPileupAnalyzer() = delete;
0056   FFTJetPileupAnalyzer(const FFTJetPileupAnalyzer&) = delete;
0057   FFTJetPileupAnalyzer& operator=(const FFTJetPileupAnalyzer&) = delete;
0058   ~FFTJetPileupAnalyzer() override;
0059
0060 private:
0061   // The following method should take all necessary info from
0062   // PileupSummaryInfo and fill out the ntuple
0063   void analyzePileup(const std::vector<PileupSummaryInfo>& pInfo);
0064
0065   void beginJob() override;
0066   void analyze(const edm::Event&, const edm::EventSetup&) override;
0067   void endJob() override;
0068
0069   edm::InputTag histoLabel;
0070   edm::InputTag summaryLabel;
0071   edm::InputTag fastJetRhoLabel;
0072   edm::InputTag fastJetSigmaLabel;
0073   edm::InputTag gridLabel;
0074   edm::InputTag srcPVs;
0075   std::string pileupLabel;
0076
0077   edm::EDGetTokenT<TH2D> histoToken;
0078   edm::EDGetTokenT<reco::FFTJetPileupSummary> summaryToken;
0079   edm::EDGetTokenT<double> fastJetRhoToken;
0080   edm::EDGetTokenT<double> fastJetSigmaToken;
0081   edm::EDGetTokenT<reco::DiscretizedEnergyFlow> gridToken;
0082   edm::EDGetTokenT<reco::VertexCollection> srcPVsToken;
0083   edm::EDGetTokenT<std::vector<PileupSummaryInfo> > pileupToken;
0084   edm::EDGetTokenT<std::pair<double, double> > etSumToken;
0085
0086   std::string ntupleName;
0087   std::string ntupleTitle;
0088   bool collectHistos;
0089   bool collectSummaries;
0090   bool collectFastJetRho;
0091   bool collectPileup;
0092   bool collectOOTPileup;
0093   bool collectGrids;
0094   bool collectGridDensity;
0095   bool collectVertexInfo;
0096   bool verbosePileupInfo;
0097
0098   double vertexNdofCut;
0099   double crazyEnergyCut;
0100
0101   std::vector<float> ntupleData;
0102   TNtuple* nt;
0103   int totalNpu;
0104   int totalNPV;
0105   unsigned long counter;
0106 };
0107
0108 //
0109 // constructors and destructor
0110 //
0111 FFTJetPileupAnalyzer::FFTJetPileupAnalyzer(const edm::ParameterSet& ps)
0112     : init_param(edm::InputTag, histoLabel),
0113       init_param(edm::InputTag, summaryLabel),
0114       init_param(edm::InputTag, fastJetRhoLabel),
0115       init_param(edm::InputTag, fastJetSigmaLabel),
0116       init_param(edm::InputTag, gridLabel),
0117       init_param(edm::InputTag, srcPVs),
0118       init_param(std::string, pileupLabel),
0119       init_param(std::string, ntupleName),
0120       init_param(std::string, ntupleTitle),
0121       init_param(bool, collectHistos),
0122       init_param(bool, collectSummaries),
0123       init_param(bool, collectFastJetRho),
0124       init_param(bool, collectPileup),
0125       init_param(bool, collectOOTPileup),
0126       init_param(bool, collectGrids),
0127       init_param(bool, collectGridDensity),
0128       init_param(bool, collectVertexInfo),
0129       init_param(bool, verbosePileupInfo),
0130       init_param(double, vertexNdofCut),
0131       init_param(double, crazyEnergyCut),
0132       nt(nullptr),
0133       totalNpu(-1),
0134       totalNPV(-1),
0135       counter(0) {
0136   usesResource(TFileService::kSharedResource);
0137
0138   if (collectPileup || collectOOTPileup)
0139     pileupToken = consumes<std::vector<PileupSummaryInfo> >(pileupLabel);
0140
0141   if (collectHistos)
0142     histoToken = consumes<TH2D>(histoLabel);
0143
0144   if (collectSummaries)
0145     summaryToken = consumes<reco::FFTJetPileupSummary>(summaryLabel);
0146
0147   if (collectFastJetRho) {
0148     fastJetRhoToken = consumes<double>(fastJetRhoLabel);
0149     fastJetSigmaToken = consumes<double>(fastJetSigmaLabel);
0150   }
0151
0152   if (collectGrids)
0153     gridToken = consumes<reco::DiscretizedEnergyFlow>(gridLabel);
0154
0155   if (collectGridDensity)
0156     etSumToken = consumes<std::pair<double, double> >(histoLabel);
0157
0158   if (collectVertexInfo)
0159     srcPVsToken = consumes<reco::VertexCollection>(srcPVs);
0160 }
0161
0162 FFTJetPileupAnalyzer::~FFTJetPileupAnalyzer() {}
0163
0164 //
0165 // member functions
0166 //
0167 void FFTJetPileupAnalyzer::analyzePileup(const std::vector<PileupSummaryInfo>& info) {
0168   const unsigned nBx = info.size();
0169   if (collectPileup)
0170     ntupleData.push_back(static_cast<float>(nBx));
0171
0172   double sumpt_Lo = 0.0, sumpt_Hi = 0.0;
0173   totalNpu = 0;
0174
0175   int npu_by_Bx[3] = {
0176       0,
0177   };
0178   double sumpt_Lo_by_Bx[3] =
0179       {
0180           0.0,
0181       },
0182          sumpt_Hi_by_Bx[3] = {
0183              0.0,
0184          };
0185
0186   if (verbosePileupInfo)
0187     std::cout << "\n**** Pileup info begin" << std::endl;
0188
0189   bool isCrazy = false;
0190   for (unsigned ibx = 0; ibx < nBx; ++ibx) {
0191     const PileupSummaryInfo& puInfo(info[ibx]);
0192
0193     const int bx = puInfo.getBunchCrossing();
0194     const int npu = puInfo.getPU_NumInteractions();
0195     const std::vector<float>& lopt(puInfo.getPU_sumpT_lowpT());
0196     const std::vector<float>& hipt(puInfo.getPU_sumpT_highpT());
0197     const double losum = std::accumulate(lopt.begin(), lopt.end(), 0.0);
0198     const double hisum = std::accumulate(hipt.begin(), hipt.end(), 0.0);
0199
0200     if (losum >= crazyEnergyCut)
0201       isCrazy = true;
0202     if (hisum >= crazyEnergyCut)
0203       isCrazy = true;
0204
0205     totalNpu += npu;
0206     sumpt_Lo += losum;
0207     sumpt_Hi += hisum;
0208
0209     const unsigned idx = bx < 0 ? 0U : (bx == 0 ? 1U : 2U);
0210     npu_by_Bx[idx] += npu;
0211     sumpt_Lo_by_Bx[idx] += losum;
0212     sumpt_Hi_by_Bx[idx] += hisum;
0213
0214     if (verbosePileupInfo)
0215       std::cout << "ibx " << ibx << " bx " << bx << " npu " << npu << " losum " << losum << " hisum " << hisum
0216                 << std::endl;
0217   }
0218
0219   if (verbosePileupInfo)
0220     std::cout << "**** Pileup info end\n" << std::endl;
0221
0222   if (isCrazy) {
0223     totalNpu = -1;
0224     sumpt_Lo = 0.0;
0225     sumpt_Hi = 0.0;
0226     for (unsigned ibx = 0; ibx < 3; ++ibx) {
0227       npu_by_Bx[ibx] = -1;
0228       sumpt_Lo_by_Bx[ibx] = 0.0;
0229       sumpt_Hi_by_Bx[ibx] = 0.0;
0230     }
0231   }
0232
0233   if (collectPileup) {
0234     ntupleData.push_back(totalNpu);
0235     ntupleData.push_back(sumpt_Lo);
0236     ntupleData.push_back(sumpt_Hi);
0237   }
0238
0239   if (collectOOTPileup)
0240     for (unsigned ibx = 0; ibx < 3; ++ibx) {
0241       ntupleData.push_back(npu_by_Bx[ibx]);
0242       ntupleData.push_back(sumpt_Lo_by_Bx[ibx]);
0243       ntupleData.push_back(sumpt_Hi_by_Bx[ibx]);
0244     }
0245 }
0246
0247 // ------------ method called once each job just before starting event loop
0248 void FFTJetPileupAnalyzer::beginJob() {
0249   // Come up with the list of variables
0250   std::string vars = "cnt:run:event";
0251   if (collectPileup)
0252     vars += ":nbx:npu:sumptLowCut:sumptHiCut";
0253   if (collectOOTPileup) {
0254     vars += ":npu_negbx:sumptLowCut_negbx:sumptHiCut_negbx";
0255     vars += ":npu_0bx:sumptLowCut_0bx:sumptHiCut_0bx";
0256     vars += ":npu_posbx:sumptLowCut_posbx:sumptHiCut_posbx";
0257   }
0258   if (collectSummaries)
0259     vars += ":estimate:pileup:uncert:uncertCode";
0260   if (collectFastJetRho)
0261     vars += ":fjrho:fjsigma";
0262   if (collectGridDensity)
0263     vars += ":gridEtDensity:gridEtDensityMixed";
0264   if (collectVertexInfo)
0265     vars += ":nPV";
0266
0267   // Book the ntuple
0268   edm::Service<TFileService> fs;
0269   nt = fs->make<TNtuple>(ntupleName.c_str(), ntupleTitle.c_str(), vars.c_str());
0270   ntupleData.reserve(nt->GetNvar());
0271 }
0272
0273 // ------------ method called to for each event  ------------
0274 void FFTJetPileupAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0275   ntupleData.clear();
0276   ntupleData.push_back(counter);
0277   totalNpu = -1;
0278   totalNPV = -1;
0279
0280   edm::RunNumber_t const runnumber = iEvent.id().run();
0281   edm::EventNumber_t const eventnumber = iEvent.id().event();
0282   ntupleData.push_back(runnumber);
0283   ntupleData.push_back(eventnumber);
0284
0285   // Get pileup information from the pile-up information module
0286   if (collectPileup || collectOOTPileup) {
0287     edm::Handle<std::vector<PileupSummaryInfo> > puInfo;
0288     if (iEvent.getByToken(pileupToken, puInfo))
0289       analyzePileup(*puInfo);
0290     else {
0291       if (collectPileup) {
0292         ntupleData.push_back(-1);
0293         ntupleData.push_back(-1);
0294         ntupleData.push_back(0.f);
0295         ntupleData.push_back(0.f);
0296       }
0297       if (collectOOTPileup)
0298         for (unsigned ibx = 0; ibx < 3; ++ibx) {
0299           ntupleData.push_back(-1);
0300           ntupleData.push_back(0.f);
0301           ntupleData.push_back(0.f);
0302         }
0303     }
0304   }
0305
0306   if (collectHistos) {
0307     edm::Handle<TH2D> input;
0308     iEvent.getByToken(histoToken, input);
0309
0310     edm::Service<TFileService> fs;
0311     TH2D* copy = new TH2D(*input);
0312
0313     std::ostringstream os;
0314     os << copy->GetName() << '_' << counter << '_' << totalNpu << '_' << runnumber << '_' << eventnumber;
0315     const std::string& newname(os.str());
0316     copy->SetNameTitle(newname.c_str(), newname.c_str());
0317
0318     copy->SetDirectory(fs->getBareDirectory());
0319   }
0320
0321   if (collectSummaries) {
0322     edm::Handle<reco::FFTJetPileupSummary> summary;
0323     iEvent.getByToken(summaryToken, summary);
0324
0325     ntupleData.push_back(summary->uncalibratedQuantile());
0326     ntupleData.push_back(summary->pileupRho());
0327     ntupleData.push_back(summary->pileupRhoUncertainty());
0328     ntupleData.push_back(summary->uncertaintyCode());
0329   }
0330
0331   if (collectFastJetRho) {
0332     edm::Handle<double> fjrho, fjsigma;
0333     iEvent.getByToken(fastJetRhoToken, fjrho);
0334     iEvent.getByToken(fastJetSigmaToken, fjsigma);
0335
0336     ntupleData.push_back(*fjrho);
0337     ntupleData.push_back(*fjsigma);
0338   }
0339
0340   if (collectGrids) {
0341     edm::Handle<reco::DiscretizedEnergyFlow> input;
0342     iEvent.getByToken(gridToken, input);
0343
0344     // Make sure the input grid is reasonable
0345     const double* data = input->data();
0346     assert(data);
0347     assert(input->phiBin0Edge() == 0.0);
0348     const unsigned nEta = input->nEtaBins();
0349     const unsigned nPhi = input->nPhiBins();
0350
0351     // Generate a name for the output histogram
0352     std::ostringstream os;
0353     os << "FFTJetGrid_" << counter << '_' << totalNpu << '_' << runnumber << '_' << eventnumber;
0354     const std::string& newname(os.str());
0355
0356     // Make a histogram and copy the grid data into it
0357     edm::Service<TFileService> fs;
0358     TH2F* h =
0359         fs->make<TH2F>(newname.c_str(), newname.c_str(), nEta, input->etaMin(), input->etaMax(), nPhi, 0.0, 2.0 * M_PI);
0360     h->GetXaxis()->SetTitle("Eta");
0361     h->GetYaxis()->SetTitle("Phi");
0362     h->GetZaxis()->SetTitle("Transverse Energy");
0363
0364     for (unsigned ieta = 0; ieta < nEta; ++ieta)
0365       for (unsigned iphi = 0; iphi < nPhi; ++iphi)
0366         h->SetBinContent(ieta + 1U, iphi + 1U, data[ieta * nPhi + iphi]);
0367   }
0368
0369   if (collectGridDensity) {
0370     edm::Handle<std::pair<double, double> > etSum;
0371     iEvent.getByToken(etSumToken, etSum);
0372
0373     ntupleData.push_back(etSum->first);
0374     ntupleData.push_back(etSum->second);
0375   }
0376
0377   if (collectVertexInfo) {
0378     edm::Handle<reco::VertexCollection> pvCollection;
0379     iEvent.getByToken(srcPVsToken, pvCollection);
0380     totalNPV = 0;
0381     if (!pvCollection->empty())
0382       for (reco::VertexCollection::const_iterator pv = pvCollection->begin(); pv != pvCollection->end(); ++pv) {
0383         const double ndof = pv->ndof();
0384         if (!pv->isFake() && ndof > vertexNdofCut)
0385           ++totalNPV;
0386       }
0387     ntupleData.push_back(totalNPV);
0388   }
0389
0390   assert(ntupleData.size() == static_cast<unsigned>(nt->GetNvar()));
0391   nt->Fill(&ntupleData[0]);
0392
0393   ++counter;
0394 }
0395
0396 // ------------ method called once each job just after ending the event loop
0397 void FFTJetPileupAnalyzer::endJob() {}
0398
0399 //define this as a plug-in
0400 DEFINE_FWK_MODULE(FFTJetPileupAnalyzer);