CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/PhysicsTools/PatExamples/bin/PatBasicFWLiteJetAnalyzer.cc

Go to the documentation of this file.
00001 #include <memory>
00002 #include <string>
00003 #include <vector>
00004 #include <sstream>
00005 #include <fstream>
00006 #include <iostream>
00007 
00008 #include <TH1F.h>
00009 #include <TROOT.h>
00010 #include <TFile.h>
00011 #include <TSystem.h>
00012 
00013 #include "DataFormats/FWLite/interface/Event.h"
00014 #include "DataFormats/Common/interface/Handle.h"
00015 #include "DataFormats/PatCandidates/interface/Jet.h"
00016 #include "FWCore/ParameterSet/interface/ProcessDesc.h"
00017 #include "FWCore/FWLite/interface/AutoLibraryLoader.h"
00018 #include "PhysicsTools/FWLite/interface/TFileService.h"
00019 #include "FWCore/PythonParameterSet/interface/PythonProcessDesc.h"
00020 
00021 
00022 int main(int argc, char* argv[])
00023 {
00024   // ----------------------------------------------------------------------
00025   // First Part:
00026   //
00027   //  * enable the AutoLibraryLoader
00028   //  * book the histograms of interest
00029   //  * open the input file
00030   // ----------------------------------------------------------------------
00031 
00032   // load framework libraries
00033   gSystem->Load( "libFWCoreFWLite" );
00034   AutoLibraryLoader::enable();
00035 
00036   // only allow one argument for this simple example which should be the
00037   // the python cfg file
00038   if ( argc < 2 ) {
00039     std::cout << "Usage : " << argv[0] << " [parameters.py]" << std::endl;
00040     return 0;
00041   }
00042  
00043   // get the python configuration
00044   PythonProcessDesc builder(argv[1]);
00045   const edm::ParameterSet& fwliteParameters = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("FWLiteParams");
00046 
00047   // now get each parameter
00048   std::string   input_ ( fwliteParameters.getParameter<std::string  >("inputFile"  ) );
00049   std::string   output_( fwliteParameters.getParameter<std::string  >("outputFile" ) );
00050   edm::InputTag jets_  ( fwliteParameters.getParameter<edm::InputTag>("jets") );
00051 
00052   // book a set of histograms
00053   fwlite::TFileService fs = fwlite::TFileService(output_.c_str());
00054   TFileDirectory theDir = fs.mkdir("analyzeBasicPat");
00055   TH1F* jetPt_  = theDir.make<TH1F>("jetPt", "pt",    100,  0.,300.);
00056   TH1F* jetEta_ = theDir.make<TH1F>("jetEta","eta",   100, -3.,  3.);
00057   TH1F* jetPhi_ = theDir.make<TH1F>("jetPhi","phi",   100, -5.,  5.);
00058   TH1F* disc_   = theDir.make<TH1F>("disc", "Discriminant", 100, 0.0, 10.0);
00059   TH1F* constituentPt_ = theDir.make<TH1F>("constituentPt", "Constituent pT", 100, 0, 300.0);
00060  
00061   // open input file (can be located on castor)
00062   TFile* inFile = TFile::Open(input_.c_str());
00063 
00064   // ----------------------------------------------------------------------
00065   // Second Part:
00066   //
00067   //  * loop the events in the input file
00068   //  * receive the collections of interest via fwlite::Handle
00069   //  * fill the histograms
00070   //  * after the loop close the input file
00071   // ----------------------------------------------------------------------
00072 
00073   // loop the events
00074   unsigned int iEvent=0;
00075   fwlite::Event ev(inFile);
00076   for(ev.toBegin(); !ev.atEnd(); ++ev, ++iEvent){
00077     edm::EventBase const & event = ev;
00078 
00079     // break loop after end of file is reached
00080     // or after 1000 events have been processed
00081     if( iEvent==1000 ) break;
00082    
00083     // simple event counter
00084     if(iEvent>0 && iEvent%1==0){
00085       std::cout << "  processing event: " << iEvent << std::endl;
00086     }
00087 
00088     // Handle to the jet collection
00089     edm::Handle<std::vector<pat::Jet> > jets;
00090     edm::InputTag jetLabel( argv[3] );
00091     event.getByLabel(jets_, jets);
00092    
00093     // loop jet collection and fill histograms
00094     for(unsigned i=0; i<jets->size(); ++i){
00095       // basic kinematics
00096       jetPt_ ->Fill( (*jets)[i].pt()  );
00097       jetEta_->Fill( (*jets)[i].eta() );
00098       jetPhi_->Fill( (*jets)[i].phi() );
00099       // access tag infos
00100       reco::SecondaryVertexTagInfo const *svTagInfos = (*jets)[i].tagInfoSecondaryVertex("secondaryVertex");
00101       if( svTagInfos != 0 ) {
00102         if( svTagInfos->nVertices() > 0 ){
00103           disc_->Fill( svTagInfos->flightDistance(0).value() );
00104         }
00105       }
00106       // access calo towers
00107       std::vector<reco::PFCandidatePtr> const & pfConstituents =  (*jets)[i].getPFConstituents();
00108       for( std::vector<reco::PFCandidatePtr>::const_iterator ibegin=pfConstituents.begin(), iend=pfConstituents.end(), iconstituent=ibegin; iconstituent!=iend; ++iconstituent){
00109         constituentPt_->Fill( (*iconstituent)->pt() );
00110       }
00111     }
00112   } 
00113   // close input file
00114   inFile->Close();
00115   
00116   // ----------------------------------------------------------------------
00117   // Third Part:
00118   //
00119   //  * never forget to free the memory of objects you created
00120   // ----------------------------------------------------------------------
00121 
00122   // in this example there is nothing to do
00123  
00124   // that's it!
00125   return 0;
00126 }
00127