PatExamples/bin/PatBasicFWLiteJetUnitTest.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/Common/interface/Handle.h"
0014 #include "DataFormats/FWLite/interface/Event.h"
0015 #include "DataFormats/PatCandidates/interface/Jet.h"
0016 #include "FWCore/FWLite/interface/FWLiteEnabler.h"
0017 #include "PhysicsTools/FWLite/interface/TFileService.h"
0018 #include "TStopwatch.h"
0019
0020 int main(int argc, char* argv[]) {
0021   // ----------------------------------------------------------------------
0022   // First Part:
0023   //
0024   //  * enable FWLite
0025   //  * book the histograms of interest
0026   //  * open the input file
0027   // ----------------------------------------------------------------------
0028
0029   if (argc < 4)
0030     return 0;
0031
0032   // load framework libraries
0033   gSystem->Load("libFWCoreFWLite");
0034   FWLiteEnabler::enable();
0035
0036   // book a set of histograms
0037   fwlite::TFileService fs = fwlite::TFileService(argv[2]);
0038   TFileDirectory theDir = fs.mkdir("analyzeBasicPat");
0039   TH1F* jetPt_ = theDir.make<TH1F>("jetPt", "pt", 100, 0., 300.);
0040   TH1F* jetEta_ = theDir.make<TH1F>("jetEta", "eta", 100, -3., 3.);
0041   TH1F* jetPhi_ = theDir.make<TH1F>("jetPhi", "phi", 100, -5., 5.);
0042   TH1F* disc_ = theDir.make<TH1F>("disc", "Discriminant", 100, 0.0, 10.0);
0043   TH1F* constituentPt_ = theDir.make<TH1F>("constituentPt", "Constituent pT", 100, 0, 300.0);
0044
0045   // open input file (can be located on castor)
0046   TFile* inFile = TFile::Open(argv[1]);
0047
0048   // ----------------------------------------------------------------------
0049   // Second Part:
0050   //
0051   //  * loop the events in the input file
0052   //  * receive the collections of interest via fwlite::Handle
0053   //  * fill the histograms
0054   //  * after the loop close the input file
0055   // ----------------------------------------------------------------------
0056
0057   // loop the events
0058   unsigned int iEvent = 0;
0059   fwlite::Event ev(inFile);
0060   TStopwatch timer;
0061   timer.Start();
0062
0063   unsigned int nEventsAnalyzed = 0;
0064   for (ev.toBegin(); !ev.atEnd(); ++ev, ++iEvent) {
0065     edm::EventBase const& event = ev;
0066
0067     // Handle to the jet collection
0068     edm::Handle<std::vector<pat::Jet> > jets;
0069     edm::InputTag jetLabel(argv[3]);
0070     event.getByLabel(jetLabel, jets);
0071
0072     // loop jet collection and fill histograms
0073     for (unsigned i = 0; i < jets->size(); ++i) {
0074       jetPt_->Fill((*jets)[i].pt());
0075       jetEta_->Fill((*jets)[i].eta());
0076       jetPhi_->Fill((*jets)[i].phi());
0077       reco::SecondaryVertexTagInfo const* svTagInfos = (*jets)[i].tagInfoSecondaryVertex("secondaryVertex");
0078       if (svTagInfos != nullptr) {
0079         if (svTagInfos->nVertices() > 0)
0080           disc_->Fill(svTagInfos->flightDistance(0).value());
0081       }
0082       std::vector<CaloTowerPtr> const& caloConstituents = (*jets)[i].getCaloConstituents();
0083       for (std::vector<CaloTowerPtr>::const_iterator ibegin = caloConstituents.begin(),
0084                                                      iend = caloConstituents.end(),
0085                                                      iconstituent = ibegin;
0086            iconstituent != iend;
0087            ++iconstituent) {
0088         constituentPt_->Fill((*iconstituent)->pt());
0089       }
0090     }
0091     ++nEventsAnalyzed;
0092   }
0093   // close input file
0094   inFile->Close();
0095
0096   timer.Stop();
0097
0098   // print some timing statistics
0099   Double_t rtime = timer.RealTime();
0100   Double_t ctime = timer.CpuTime();
0101   printf("Analyzed events: %d \n", nEventsAnalyzed);
0102   printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime);
0103   printf("%4.2f events / RealTime second .\n", (double)nEventsAnalyzed / rtime);
0104   printf("%4.2f events / CpuTime second .\n", (double)nEventsAnalyzed / ctime);
0105
0106   // ----------------------------------------------------------------------
0107   // Third Part:
0108   //
0109   //  * never forget to free the memory of objects you created
0110   // ----------------------------------------------------------------------
0111
0112   // in this example there is nothing to do
0113
0114   // that's it!
0115   return 0;
0116 }