FWLite/bin/FWLiteWithSelectorUtils.cc

0001 #include <TH1F.h>
0002 #include <TROOT.h>
0003 #include <TFile.h>
0004 #include <TSystem.h>
0005
0006 #include "DataFormats/FWLite/interface/Event.h"
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "FWCore/FWLite/interface/FWLiteEnabler.h"
0009
0010 #include "DataFormats/FWLite/interface/InputSource.h"
0011 #include "DataFormats/FWLite/interface/OutputFiles.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/ParameterSetReader/interface/ParameterSetReader.h"
0014
0015 #include "TStopwatch.h"
0016 #include "PhysicsTools/FWLite/interface/WSelector.h"
0017 #include "PhysicsTools/FWLite/interface/TFileService.h"
0018 //#include "PhysicsTools/FWLite/interface/WSelectorFast.h"
0019
0020 int main(int argc, char* argv[]) {
0021   // define what muon you are using; this is necessary as FWLite is not
0022   // capable of reading edm::Views
0023   using reco::Muon;
0024
0025   // ----------------------------------------------------------------------
0026   // First Part:
0027   //
0028   //  * enable FWLite
0029   //  * book the histograms of interest
0030   //  * open the input file
0031   // ----------------------------------------------------------------------
0032
0033   // load framework libraries
0034   gSystem->Load("libFWCoreFWLite");
0035   FWLiteEnabler::enable();
0036
0037   if (argc < 2) {
0038     std::cout << "Usage : " << argv[0] << " [parameters.py]" << std::endl;
0039     return 0;
0040   }
0041
0042   if (!edm::readPSetsFrom(argv[1])->existsAs<edm::ParameterSet>("process")) {
0043     std::cout << " ERROR: ParametersSet 'process' is missing in your configuration file" << std::endl;
0044     exit(0);
0045   }
0046   // get the python configuration
0047   const edm::ParameterSet& process = edm::readPSetsFrom(argv[1])->getParameter<edm::ParameterSet>("process");
0048   fwlite::InputSource inputHandler_(process);
0049   fwlite::OutputFiles outputHandler_(process);
0050
0051   // initialize the W selector
0052   edm::ParameterSet selection = process.getParameter<edm::ParameterSet>("selection");
0053   WSelector wSelector(selection);
0054   pat::strbitset wSelectorReturns = wSelector.getBitTemplate();
0055
0056   // book a set of histograms
0057   fwlite::TFileService fs = fwlite::TFileService(outputHandler_.file());
0058   TFileDirectory theDir = fs.mkdir("analyzeBasicPat");
0059   TH1F* muonPt_ = theDir.make<TH1F>("muonPt", "pt", 100, 0., 300.);
0060   TH1F* muonEta_ = theDir.make<TH1F>("muonEta", "eta", 100, -3., 3.);
0061   TH1F* muonPhi_ = theDir.make<TH1F>("muonPhi", "phi", 100, -5., 5.);
0062
0063   // start a CPU timer
0064   TStopwatch timer;
0065   timer.Start();
0066
0067   // loop the events
0068   int ievt = 0;
0069   unsigned int nEventsAnalyzed = 0;
0070   int maxEvents_(inputHandler_.maxEvents());
0071   for (unsigned int iFile = 0; iFile < inputHandler_.files().size(); ++iFile) {
0072     // open input file (can be located on castor)
0073     TFile* inFile = TFile::Open(inputHandler_.files()[iFile].c_str());
0074     if (inFile) {
0075       // ----------------------------------------------------------------------
0076       // Second Part:
0077       //
0078       //  * loop the events in the input file
0079       //  * receive the collections of interest via fwlite::Handle
0080       //  * fill the histograms
0081       //  * after the loop close the input file
0082       // ----------------------------------------------------------------------
0083       fwlite::Event ev(inFile);
0084       for (ev.toBegin(); !ev.atEnd(); ++ev, ++ievt) {
0085         edm::EventBase const& event = ev;
0086         // break loop if maximal number of events is reached
0087         if (maxEvents_ > 0 ? ievt + 1 > maxEvents_ : false)
0088           break;
0089         // simple event counter
0090         if (inputHandler_.reportAfter() != 0 ? (ievt > 0 && ievt % inputHandler_.reportAfter() == 0) : false)
0091           std::cout << "  processing event: " << ievt << std::endl;
0092
0093         if (wSelector(event, wSelectorReturns)) {
0094           pat::Muon const& wMuon = wSelector.wMuon();
0095           muonPt_->Fill(wMuon.pt());
0096           muonEta_->Fill(wMuon.eta());
0097           muonPhi_->Fill(wMuon.phi());
0098         }
0099         ++nEventsAnalyzed;
0100       }
0101       // close input file
0102       inFile->Close();
0103     }
0104     // break loop if maximal number of events is reached:
0105     // this has to be done twice to stop the file loop as well
0106     if (maxEvents_ > 0 ? ievt + 1 > maxEvents_ : false)
0107       break;
0108   }
0109   // stop CPU timer
0110   timer.Stop();
0111
0112   // print selector
0113   wSelector.print(std::cout);
0114
0115   // print some timing statistics
0116   double rtime = timer.RealTime();
0117   double ctime = timer.CpuTime();
0118   // timing printouts
0119   printf("Analyzed events: %d \n", nEventsAnalyzed);
0120   printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime);
0121   printf("%4.2f events / RealTime second .\n", (double)nEventsAnalyzed / rtime);
0122   printf("%4.2f events / CpuTime second .\n", (double)nEventsAnalyzed / ctime);
0123   return 0;
0124 }