FWLite/bin/FWLiteWithPythonConfig.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 "DataFormats/MuonReco/interface/Muon.h"
0016 #include "DataFormats/PatCandidates/interface/Muon.h"
0017 #include "PhysicsTools/FWLite/interface/TFileService.h"
0018
0019 int main(int argc, char* argv[]) {
0020   // define what muon you are using; this is necessary as FWLite is not
0021   // capable of reading edm::Views
0022   using reco::Muon;
0023
0024   // ----------------------------------------------------------------------
0025   // First Part:
0026   //
0027   //  * enable FWLite
0028   //  * book the histograms of interest
0029   //  * open the input file
0030   // ----------------------------------------------------------------------
0031
0032   // load framework libraries
0033   gSystem->Load("libFWCoreFWLite");
0034   FWLiteEnabler::enable();
0035
0036   // parse arguments
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   // now get each parameter
0052   const edm::ParameterSet& ana = process.getParameter<edm::ParameterSet>("muonAnalyzer");
0053   edm::InputTag muons_(ana.getParameter<edm::InputTag>("muons"));
0054
0055   // book a set of histograms
0056   fwlite::TFileService fs = fwlite::TFileService(outputHandler_.file());
0057   TFileDirectory dir = fs.mkdir("analyzeBasicPat");
0058   TH1F* muonPt_ = dir.make<TH1F>("muonPt", "pt", 100, 0., 300.);
0059   TH1F* muonEta_ = dir.make<TH1F>("muonEta", "eta", 100, -3., 3.);
0060   TH1F* muonPhi_ = dir.make<TH1F>("muonPhi", "phi", 100, -5., 5.);
0061   TH1F* mumuMass_ = dir.make<TH1F>("mumuMass", "mass", 90, 30., 120.);
0062
0063   // loop the events
0064   int ievt = 0;
0065   int maxEvents_(inputHandler_.maxEvents());
0066   for (unsigned int iFile = 0; iFile < inputHandler_.files().size(); ++iFile) {
0067     // open input file (can be located on castor)
0068     TFile* inFile = TFile::Open(inputHandler_.files()[iFile].c_str());
0069     if (inFile) {
0070       // ----------------------------------------------------------------------
0071       // Second Part:
0072       //
0073       //  * loop the events in the input file
0074       //  * receive the collections of interest via fwlite::Handle
0075       //  * fill the histograms
0076       //  * after the loop close the input file
0077       // ----------------------------------------------------------------------
0078       fwlite::Event ev(inFile);
0079       for (ev.toBegin(); !ev.atEnd(); ++ev, ++ievt) {
0080         edm::EventBase const& event = ev;
0081         // break loop if maximal number of events is reached
0082         if (maxEvents_ > 0 ? ievt + 1 > maxEvents_ : false)
0083           break;
0084         // simple event counter
0085         if (inputHandler_.reportAfter() != 0 ? (ievt > 0 && ievt % inputHandler_.reportAfter() == 0) : false)
0086           std::cout << "  processing event: " << ievt << std::endl;
0087
0088         // Handle to the muon collection
0089         edm::Handle<std::vector<Muon> > muons;
0090         event.getByLabel(muons_, muons);
0091
0092         // loop muon collection and fill histograms
0093         for (std::vector<Muon>::const_iterator mu1 = muons->begin(); mu1 != muons->end(); ++mu1) {
0094           muonPt_->Fill(mu1->pt());
0095           muonEta_->Fill(mu1->eta());
0096           muonPhi_->Fill(mu1->phi());
0097           if (mu1->pt() > 20 && fabs(mu1->eta()) < 2.1) {
0098             for (std::vector<Muon>::const_iterator mu2 = muons->begin(); mu2 != muons->end(); ++mu2) {
0099               if (mu2 > mu1) {                            // prevent double conting
0100                 if (mu1->charge() * mu2->charge() < 0) {  // check only muon pairs of unequal charge
0101                   if (mu2->pt() > 20 && fabs(mu2->eta()) < 2.1) {
0102                     mumuMass_->Fill((mu1->p4() + mu2->p4()).mass());
0103                   }
0104                 }
0105               }
0106             }
0107           }
0108         }
0109       }
0110       // close input file
0111       inFile->Close();
0112     }
0113     // break loop if maximal number of events is reached:
0114     // this has to be done twice to stop the file loop as well
0115     if (maxEvents_ > 0 ? ievt + 1 > maxEvents_ : false)
0116       break;
0117   }
0118   return 0;
0119 }