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