PatExamples/bin/PatBasicFWLiteJetAnalyzer.cc

0001 #include <memory>
0002 #include <string>
0003 #include <vector>
0004 #include <sstream>
0005 #include <fstream>
0006 #include <iostream>
0007
0008 #include <TH1F.h>
0009 #include <TROOT.h>
0010 #include <TFile.h>
0011 #include <TSystem.h>
0012
0013 #include "DataFormats/FWLite/interface/Event.h"
0014 #include "DataFormats/Common/interface/Handle.h"
0015 #include "DataFormats/PatCandidates/interface/Jet.h"
0016 #include "FWCore/ParameterSet/interface/ProcessDesc.h"
0017 #include "FWCore/FWLite/interface/FWLiteEnabler.h"
0018 #include "PhysicsTools/FWLite/interface/TFileService.h"
0019 #include "FWCore/ParameterSetReader/interface/ProcessDescImpl.h"
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   edm::InputTag jets_(fwliteParameters.getParameter<edm::InputTag>("jets"));
0050
0051   // book a set of histograms
0052   fwlite::TFileService fs = fwlite::TFileService(output_);
0053   TFileDirectory theDir = fs.mkdir("analyzeBasicPat");
0054   TH1F* jetPt_ = theDir.make<TH1F>("jetPt", "pt", 100, 0., 300.);
0055   TH1F* jetEta_ = theDir.make<TH1F>("jetEta", "eta", 100, -3., 3.);
0056   TH1F* jetPhi_ = theDir.make<TH1F>("jetPhi", "phi", 100, -5., 5.);
0057   TH1F* disc_ = theDir.make<TH1F>("disc", "Discriminant", 100, 0.0, 10.0);
0058   TH1F* constituentPt_ = theDir.make<TH1F>("constituentPt", "Constituent pT", 100, 0, 300.0);
0059
0060   // open input file (can be located on castor)
0061   TFile* inFile = TFile::Open(input_.c_str());
0062
0063   // ----------------------------------------------------------------------
0064   // Second Part:
0065   //
0066   //  * loop the events in the input file
0067   //  * receive the collections of interest via fwlite::Handle
0068   //  * fill the histograms
0069   //  * after the loop close the input file
0070   // ----------------------------------------------------------------------
0071
0072   // loop the events
0073   unsigned int iEvent = 0;
0074   fwlite::Event ev(inFile);
0075   for (ev.toBegin(); !ev.atEnd(); ++ev, ++iEvent) {
0076     edm::EventBase const& event = ev;
0077
0078     // break loop after end of file is reached
0079     // or after 1000 events have been processed
0080     if (iEvent == 1000)
0081       break;
0082
0083     // simple event counter
0084     if (iEvent > 0 && iEvent % 1 == 0) {
0085       std::cout << "  processing event: " << iEvent << std::endl;
0086     }
0087
0088     // Handle to the jet collection
0089     edm::Handle<std::vector<pat::Jet> > jets;
0090     edm::InputTag jetLabel(argv[3]);
0091     event.getByLabel(jets_, jets);
0092
0093     // loop jet collection and fill histograms
0094     for (unsigned i = 0; i < jets->size(); ++i) {
0095       // basic kinematics
0096       jetPt_->Fill((*jets)[i].pt());
0097       jetEta_->Fill((*jets)[i].eta());
0098       jetPhi_->Fill((*jets)[i].phi());
0099       // access tag infos
0100       reco::SecondaryVertexTagInfo const* svTagInfos = (*jets)[i].tagInfoSecondaryVertex("secondaryVertex");
0101       if (svTagInfos != nullptr) {
0102         if (svTagInfos->nVertices() > 0) {
0103           disc_->Fill(svTagInfos->flightDistance(0).value());
0104         }
0105       }
0106       // access calo towers
0107       std::vector<reco::PFCandidatePtr> const& pfConstituents = (*jets)[i].getPFConstituents();
0108       for (std::vector<reco::PFCandidatePtr>::const_iterator ibegin = pfConstituents.begin(),
0109                                                              iend = pfConstituents.end(),
0110                                                              iconstituent = ibegin;
0111            iconstituent != iend;
0112            ++iconstituent) {
0113         constituentPt_->Fill((*iconstituent)->pt());
0114       }
0115     }
0116   }
0117   // close input file
0118   inFile->Close();
0119
0120   // ----------------------------------------------------------------------
0121   // Third Part:
0122   //
0123   //  * never forget to free the memory of objects you created
0124   // ----------------------------------------------------------------------
0125
0126   // in this example there is nothing to do
0127
0128   // that's it!
0129   return 0;
0130 }