PatExamples/test/TestCorrections_FWLiteMacro.C

0001 // -*- C++ -*-
0002 //****Author -  Sudhir Malik - malik@fnal.gov *****//
0003
0004   //*****NOTE: To successfully exeute this macro, include the following lines in your root logon file ************************//
0005     /*
0006       {
0007       gSystem->Load("libFWCoreFWLite.so");
0008       FWLiteEnabler::enable();
0009       gSystem->Load("libDataFormatsFWLite.so");
0010       gSystem->Load("libDataFormatsPatCandidates.so");
0011       gSystem->Load("libRooFit") ;
0012       using namespace RooFit ;
0013       cout << "loaded" << endl;
0014       }
0015     */
0016 //*****************************//
0017
0018   // CMS includes
0019 #include "DataFormats/FWLite/interface/Handle.h"
0020 #include "DataFormats/FWLite/interface/Event.h"
0021 #include <TH2.h>
0022
0023
0024 #include "TFile.h"
0025 #include "TH1.h"
0026 #include "TCanvas.h"
0027 #include "TLegend.h"
0028
0029 #if !defined(__CINT__) && !defined(__MAKECINT__)
0030
0031 #include "DataFormats/PatCandidates/interface/Jet.h"
0032 #include "CondFormats/JetMETObjects/interface/CombinedJetCorrector.h"
0033
0034   // these includes are needed to make the  "gSystem" commands below work
0035 #include "TSystem.h"
0036 #include "TROOT.h"
0037 #include "FWCore/FWLite/interface/FWLiteEnabler.h"
0038 #endif
0039
0040
0041 #include <iostream>
0042 #include <vector>
0043 #include <memory>
0044
0045
0046 void JetCorrFWLite()
0047 {
0048
0049
0050   // The input file MYCOPY.root  is available at
0051   //https://cms-service-sdtweb.web.cern.ch/cms-service-sdtweb/validation/physicstools/PATVALIDATION/MYCOPY.root
0052
0053   TFile  *file = new TFile("MYCOPY.root");
0054   TFile  *outputfile    = new TFile("JetCorr.root","RECREATE");
0055
0056
0057
0058   // Book those histograms for Jet stuff!
0059   TH1F * hist_jetMult = new TH1F( "jetMult", "jetMult", 100, 0, 50) ;
0060   TH1F * hist_jetPt = new TH1F( "jetPt", "jetPt", 100, 0, 200) ;
0061   TH1F * hist_jetEta = new TH1F( "jetEta", "jetEta", 100,-5, 5) ;
0062   TH1F * hist_jetPhi = new TH1F( "jetPhi", "jetPhi", 100, -3.5, 3.5) ;
0063   TH2F * hist_Mapping = new TH2F( "Mapping", "Mapping", 500, 0, 500, 500, 0, 500) ;
0064   TH2F * hist_Ratio = new TH2F( "Ratio", "Ratio", 500, 0, 500, 500, 0.8, 1.2) ;
0065
0066
0067   fwlite::Event ev(file);
0068   for( ev.toBegin();
0069        ! ev.atEnd();
0070        ++ev) {
0071
0072
0073     fwlite::Handle<vector<reco::CaloJet> > jetHandle;
0074     jetHandle.getByLabel(ev,"ak5CaloJets");
0075
0076
0077
0078     ////////////////////////////////////////////////////////////////////////
0079     ////////////// Defining the L2L3L5L7JetCorrector ///////////////////////
0080     ////////////////////////////////////////////////////////////////////////
0081
0082     string Levels1 = "L2:L3:L5:L7";
0083     string Tags1 = "900GeV_L2Relative_AK5Calo:900GeV_L3Absolute_AK5Calo:L5Flavor_IC5:L7Parton_IC5";
0084     string Options1 = "Flavor:gJ & Parton:gJ";
0085     CombinedJetCorrector *L2L3L5L7JetCorrector = new CombinedJetCorrector(Levels1,Tags1,Options1);
0086
0087     ///////////////////////////////////////////////////////////////////////
0088       ////////////// Defining the L2L3JetCorrector //////////////////////////
0089       ///////////////////////////////////////////////////////////////////////
0090
0091       string Levels = "L2:L3";
0092       string Tags = "900GeV_L2Relative_AK5Calo:900GeV_L3Absolute_AK5Calo";
0093       CombinedJetCorrector *L2L3JetCorrector = new CombinedJetCorrector(Levels,Tags);
0094
0095
0096
0097
0098
0099       // Loop over the jets
0100       hist_jetMult->Fill(jetHandle->size());
0101       const vector<reco::CaloJet>::const_iterator kJetEnd = jetHandle->end();
0102       for (vector<reco::CaloJet>::const_iterator jetIter = jetHandle->begin();
0103            kJetEnd != jetIter;
0104            ++jetIter)
0105     {
0106       hist_jetPt ->Fill (jetIter->pt());
0107       hist_jetEta->Fill (jetIter->eta());
0108       hist_jetPhi->Fill (jetIter->phi());
0109       ////////////////////Jet Correctiion Stuff ////////////
0110           double pt = jetIter->pt();
0111           double eta = jetIter->eta();
0112           double emf = jetIter->emEnergyFraction();
0113
0114           double L2L3scale = L2L3JetCorrector->getCorrection(pt,eta,emf);
0115           double L2L3L5L7scale = L2L3L5L7JetCorrector->getCorrection(pt,eta,emf);
0116
0117           vector<double> L2L3factors = L2L3JetCorrector->getSubCorrections(pt,eta,emf);
0118           vector<double> L2L3L5L7factors = L2L3L5L7JetCorrector->getSubCorrections(pt,eta,emf);
0119
0120           cout<<"Pt = "<<pt<<", Eta = "<<eta<<", EMF = "<<emf<<endl;
0121           cout<<"L2L3correction = "<<L2L3scale<<", L2L3CorPt = "<<L2L3scale*pt<<endl;
0122           for(unsigned int i=0;i < L2L3factors.size();i++)
0123             cout<<L2L3factors[i]<<endl;
0124           cout<<"L2L3L5L7correction = "<<L2L3L5L7scale<<", L2L3L5L7CorPt = "<<L2L3L5L7scale*pt<<endl;
0125           for(unsigned int i=0;i< L2L3L5L7factors.size();i++)
0126             cout<<L2L3L5L7factors[i]<<endl;
0127           hist_Mapping->Fill (L2L3scale*pt,L2L3L5L7scale*pt);
0128           hist_Ratio->Fill (L2L3scale*pt,L2L3L5L7scale/L2L3scale);
0129       cout << "L2L3scale*pt = " << L2L3scale*pt << endl;
0130       cout << "L2L3L5L7scale*pt = " << L2L3L5L7scale*pt << endl;
0131           //////////////////// End of Jet Correctiion Stuff ////////////
0132
0133     } // for jetIter
0134
0135
0136
0137   } // for event loop
0138
0139
0140   outputfile->cd();
0141   outputfile->Write();
0142   outputfile->Close();
0143
0144 }
0145