PatExamples/bin/PatCleaningExercise.cc

0001 #include <vector>
0002
0003 #include "TH1.h"
0004 #include "TFile.h"
0005 #include <TROOT.h>
0006 #include <TSystem.h>
0007
0008 #include "DataFormats/Math/interface/deltaR.h"
0009 #include "DataFormats/FWLite/interface/Event.h"
0010 #include "DataFormats/Common/interface/Handle.h"
0011 #include "DataFormats/PatCandidates/interface/Jet.h"
0012 #include "FWCore/ParameterSet/interface/ProcessDesc.h"
0013 #include "FWCore/FWLite/interface/FWLiteEnabler.h"
0014 #include "PhysicsTools/FWLite/interface/TFileService.h"
0015 #include "FWCore/ParameterSetReader/interface/ProcessDescImpl.h"
0016
0017 //using namespace std;
0018 //using namespace reco;
0019 //using namespace pat;
0020
0021 int main(int argc, char* argv[]) {
0022   // ----------------------------------------------------------------------
0023   // First Part:
0024   //
0025   //  * enable FWLite
0026   //  * book the histograms of interest
0027   //  * open the input file
0028   // ----------------------------------------------------------------------
0029
0030   // load framework libraries
0031   gSystem->Load("libFWCoreFWLite");
0032   FWLiteEnabler::enable();
0033
0034   // only allow one argument for this simple example which should be the
0035   // the python cfg file
0036   if (argc < 2) {
0037     std::cout << "Usage : " << argv[0] << " [parameters.py]" << std::endl;
0038     return 0;
0039   }
0040
0041   // get the python configuration
0042   ProcessDescImpl builder(argv[1], true);
0043   const edm::ParameterSet& fwliteParameters =
0044       builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("FWLiteParams");
0045
0046   // now get each parameter
0047   std::string input_(fwliteParameters.getParameter<std::string>("inputFile"));
0048   std::string output_(fwliteParameters.getParameter<std::string>("outputFile"));
0049   std::string overlaps_(fwliteParameters.getParameter<std::string>("overlaps"));
0050   edm::InputTag jets_(fwliteParameters.getParameter<edm::InputTag>("jets"));
0051
0052   // book a set of histograms
0053   fwlite::TFileService fs = fwlite::TFileService(output_);
0054   TFileDirectory theDir = fs.mkdir("analyzePatCleaning");
0055   TH1F* emfAllJets_ = theDir.make<TH1F>("emfAllJets", "f_{emf}(All Jets)", 20, 0., 1.);
0056   TH1F* emfCleanJets_ = theDir.make<TH1F>("emfCleanJets", "f_{emf}(Clean Jets)", 20, 0., 1.);
0057   TH1F* emfOverlapJets_ = theDir.make<TH1F>("emfOverlapJets", "f_{emf}(Overlap Jets)", 20, 0., 1.);
0058   TH1F* deltaRElecJet_ = theDir.make<TH1F>("deltaRElecJet", "#DeltaR (elec, jet)", 10, 0., 0.5);
0059   TH1F* elecOverJet_ = theDir.make<TH1F>("elecOverJet", "E_{elec}/E_{jet}", 100, 0., 2.);
0060   TH1F* nOverlaps_ = theDir.make<TH1F>("nOverlaps", "Number of overlaps", 5, 0., 5.);
0061
0062   // open input file (can be located on castor)
0063   TFile* inFile = TFile::Open(input_.c_str());
0064
0065   // ----------------------------------------------------------------------
0066   // Second Part:
0067   //
0068   //  * loop the events in the input file
0069   //  * receive the collections of interest via fwlite::Handle
0070   //  * fill the histograms
0071   //  * after the loop close the input file
0072   // ----------------------------------------------------------------------
0073
0074   // loop the events
0075   unsigned int iEvent = 0;
0076   fwlite::Event ev(inFile);
0077   for (ev.toBegin(); !ev.atEnd(); ++ev, ++iEvent) {
0078     edm::EventBase const& event = ev;
0079
0080     // break loop after end of file is reached
0081     // or after 1000 events have been processed
0082     if (iEvent == 1000)
0083       break;
0084
0085     // simple event counter
0086     if (iEvent > 0 && iEvent % 1 == 0) {
0087       std::cout << "  processing event: " << iEvent << std::endl;
0088     }
0089
0090     // handle to jet collection
0091     edm::Handle<std::vector<pat::Jet> > jets;
0092     event.getByLabel(jets_, jets);
0093
0094     // loop over the jets in the event
0095     for (std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet != jets->end(); jet++) {
0096       if (jet->pt() > 20 && jet == jets->begin()) {
0097         emfAllJets_->Fill(jet->emEnergyFraction());
0098         if (!jet->hasOverlaps(overlaps_)) {
0099           emfCleanJets_->Fill(jet->emEnergyFraction());
0100         } else {
0101           //get all overlaps
0102           const reco::CandidatePtrVector overlaps = jet->overlaps(overlaps_);
0103           nOverlaps_->Fill(overlaps.size());
0104           emfOverlapJets_->Fill(jet->emEnergyFraction());
0105           //loop over the overlaps
0106           for (reco::CandidatePtrVector::const_iterator overlap = overlaps.begin(); overlap != overlaps.end();
0107                overlap++) {
0108             float deltaR = reco::deltaR((*overlap)->eta(), (*overlap)->phi(), jet->eta(), jet->phi());
0109             deltaRElecJet_->Fill(deltaR);
0110             elecOverJet_->Fill((*overlap)->energy() / jet->energy());
0111           }
0112         }
0113       }
0114     }
0115   }
0116   inFile->Close();
0117   return 0;
0118 }