PatExamples/bin/PatBasicFWLiteAnalyzer.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 "FWCore/FWLite/interface/FWLiteEnabler.h"
0016
0017 #include "DataFormats/FWLite/interface/InputSource.h"
0018 #include "DataFormats/FWLite/interface/OutputFiles.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/ParameterSetReader/interface/ParameterSetReader.h"
0021
0022 #include "DataFormats/PatCandidates/interface/Muon.h"
0023 #include "PhysicsTools/FWLite/interface/TFileService.h"
0024
0025 int main(int argc, char* argv[]) {
0026   // ----------------------------------------------------------------------
0027   // First Part:
0028   //
0029   //  * enable FWLite
0030   //  * book the histograms of interest
0031   //  * open the input file
0032   // ----------------------------------------------------------------------
0033
0034   // load framework libraries
0035   gSystem->Load("libFWCoreFWLite");
0036   FWLiteEnabler::enable();
0037
0038   // only allow one argument for this simple example which should be the
0039   // the python cfg file
0040   if (argc < 2) {
0041     std::cout << "Usage : " << argv[0] << " [parameters.py]" << std::endl;
0042     return 0;
0043   }
0044   if (!edm::readPSetsFrom(argv[1])->existsAs<edm::ParameterSet>("process")) {
0045     std::cout << " ERROR: ParametersSet 'process' is missing in your configuration file" << std::endl;
0046     exit(0);
0047   }
0048   // get the python configuration
0049   const edm::ParameterSet& process = edm::readPSetsFrom(argv[1])->getParameter<edm::ParameterSet>("process");
0050   fwlite::InputSource inputHandler_(process);
0051   fwlite::OutputFiles outputHandler_(process);
0052
0053   // now get each parameter
0054   const edm::ParameterSet& ana = process.getParameter<edm::ParameterSet>("muonAnalyzer");
0055   edm::InputTag muons_(ana.getParameter<edm::InputTag>("muons"));
0056
0057   // book a set of histograms
0058   fwlite::TFileService fs = fwlite::TFileService(outputHandler_.file());
0059   TFileDirectory dir = fs.mkdir("analyzeBasicPat");
0060   TH1F* muonPt_ = dir.make<TH1F>("muonPt", "pt", 100, 0., 300.);
0061   TH1F* muonEta_ = dir.make<TH1F>("muonEta", "eta", 100, -3., 3.);
0062   TH1F* muonPhi_ = dir.make<TH1F>("muonPhi", "phi", 100, -5., 5.);
0063
0064   // loop the events
0065   int ievt = 0;
0066   int maxEvents_(inputHandler_.maxEvents());
0067   for (unsigned int iFile = 0; iFile < inputHandler_.files().size(); ++iFile) {
0068     // open input file (can be located on castor)
0069     TFile* inFile = TFile::Open(inputHandler_.files()[iFile].c_str());
0070     if (inFile) {
0071       // ----------------------------------------------------------------------
0072       // Second Part:
0073       //
0074       //  * loop the events in the input file
0075       //  * receive the collections of interest via fwlite::Handle
0076       //  * fill the histograms
0077       //  * after the loop close the input file
0078       // ----------------------------------------------------------------------
0079       fwlite::Event ev(inFile);
0080       for (ev.toBegin(); !ev.atEnd(); ++ev, ++ievt) {
0081         edm::EventBase const& event = ev;
0082         // break loop if maximal number of events is reached
0083         if (maxEvents_ > 0 ? ievt + 1 > maxEvents_ : false)
0084           break;
0085         // simple event counter
0086         if (inputHandler_.reportAfter() != 0 ? (ievt > 0 && ievt % inputHandler_.reportAfter() == 0) : false)
0087           std::cout << "  processing event: " << ievt << std::endl;
0088
0089         // Handle to the muon collection
0090         edm::Handle<std::vector<pat::Muon> > muons;
0091         event.getByLabel(muons_, muons);
0092
0093         // loop muon collection and fill histograms
0094         for (unsigned i = 0; i < muons->size(); ++i) {
0095           muonPt_->Fill((*muons)[i].pt());
0096           muonEta_->Fill((*muons)[i].eta());
0097           muonPhi_->Fill((*muons)[i].phi());
0098         }
0099       }
0100       // close input file
0101       inFile->Close();
0102     }
0103     // break loop if maximal number of events is reached:
0104     // this has to be done twice to stop the file loop as well
0105     if (maxEvents_ > 0 ? ievt + 1 > maxEvents_ : false)
0106       break;
0107   }
0108   return 0;
0109 }