/data/refman/pasoursint/CMSSW_5_3_10_patch1/src/RecoJets/JetAnalyzers/src/CMSDAS11DijetTestAnalyzer.cc

Go to the documentation of this file.

// CMSDAS11DijetTestAnalyzer.cc
// Description: A basic dijet analyzer for the CMSDAS 2011
// Author: John Paul Chou
// Date: January 12, 2011

#include "RecoJets/JetAnalyzers/interface/CMSDAS11DijetTestAnalyzer.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "DataFormats/VertexReco/interface/Vertex.h"
#include "JetMETCorrections/Objects/interface/JetCorrector.h"

#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"

#include <TH1D.h>
#include <TH2D.h>

CMSDAS11DijetTestAnalyzer::CMSDAS11DijetTestAnalyzer(edm::ParameterSet const& params) :
  edm::EDAnalyzer(),
  jetSrc(params.getParameter<edm::InputTag>("jetSrc")),
  vertexSrc(params.getParameter<edm::InputTag>("vertexSrc")),
  jetCorrections(params.getParameter<std::string>("jetCorrections")),
  innerDeltaEta(params.getParameter<double>("innerDeltaEta")),
  outerDeltaEta(params.getParameter<double>("outerDeltaEta")),
  JESbias(params.getParameter<double>("JESbias"))
{
  // setup file service
  edm::Service<TFileService> fs;

  const int NBINS=36;
  Double_t BOUNDARIES[NBINS] = { 220, 244, 270, 296, 325, 354, 386, 419, 453,
                                 489, 526, 565, 606, 649, 693, 740, 788, 838,
                                 890, 944, 1000, 1058, 1118, 1181, 1246, 1313, 1383,
                                 1455, 1530, 1607, 1687, 1770, 1856, 1945, 2037, 2132 };
  
  // setup histograms
  hVertexZ = fs->make<TH1D>("hVertexZ", "Z position of the Vertex",50,-20,20);
  hJetRawPt    = fs->make<TH1D>("hJetRawPt","Raw Jet Pt",50,0,1000);
  hJetCorrPt = fs->make<TH1D>("hJetCorrPt","Corrected Jet Pt",50,0,1000);
  hJet1Pt      = fs->make<TH1D>("hJet1Pt","Corrected Jet1 Pt",50,0,1000);
  hJet2Pt      = fs->make<TH1D>("hJet2Pt","Corrected Jet2 Pt",50,0,1000);

  hJetEta      = fs->make<TH1D>("hJetEta","Corrected Jet Eta",   50,-5,5);
  hJet1Eta      = fs->make<TH1D>("hJet1Eta","Corrected Jet1 Eta",50,-5,5);
  hJet2Eta      = fs->make<TH1D>("hJet2Eta","Corrected Jet2 Eta",50,-5,5);

  hJetPhi      = fs->make<TH1D>("hJetPhi","Corrected Jet Phi",   50,-3.1415,3.1415);
  hJet1Phi      = fs->make<TH1D>("hJet1Phi","Corrected Jet1 Phi",50,-3.1415,3.1415);
  hJet2Phi      = fs->make<TH1D>("hJet2Phi","Corrected Jet2 Phi",50,-3.1415,3.1415);

  hJetEMF       = fs->make<TH1D>("hJetEMF","EM Fraction of Jets",50,0,1);
  hJet1EMF      = fs->make<TH1D>("hJet1EMF","EM Fraction of Jet1",50,0,1);
  hJet2EMF      = fs->make<TH1D>("hJet2EMF","EM Fraction of Jet2",50,0,1);

  hCorDijetMass = fs->make<TH1D>("hCorDijetMass","Corrected Dijet Mass",NBINS-1,BOUNDARIES);
  hDijetDeltaPhi= fs->make<TH1D>("hDijetDeltaPhi","Dijet |#Delta #phi|",50,0,3.1415);
  hDijetDeltaEta= fs->make<TH1D>("hDijetDeltaEta","Dijet |#Delta #eta|",50,0,1.3);
  hDijetDeltaPhiNJets= fs->make<TH2D>("hDijetDeltaPhiNJets","Dijet |#Delta #phi| v. the number of jets",50,0,3.1415,7,0.5,7.5);
  hDijetEta1Eta2=fs->make<TH2D>("hDijetEta1Eta2","Eta 1 versus Eta 2 of dijet events",50,-5,5,50,-5,5);

  hInnerDijetMass = fs->make<TH1D>("hInnerDijetMass","Corrected Inner Dijet Mass",NBINS-1,BOUNDARIES);
  hOuterDijetMass = fs->make<TH1D>("hOuterDijetMass","Corrected Outer Dijet Mass",NBINS-1,BOUNDARIES);
}

void CMSDAS11DijetTestAnalyzer::endJob(void) {
}

void CMSDAS11DijetTestAnalyzer::analyze( const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  double mWeight;
  edm::Handle<GenEventInfoProduct> hEventInfo;
  iEvent.getByLabel("generator", hEventInfo);
  if (hEventInfo.isValid()) {
    mWeight = hEventInfo->weight();
  }
  else 
    mWeight = 1.;
  
  // Get event ID information
  //int nrun=iEvent.id().run();
  //int nlumi=iEvent.luminosityBlock();
  //int nevent=iEvent.id().event();

  // Get Primary Vertex Information

  // magic to get the vertices from EDM
  edm::Handle< std::vector<reco::Vertex> > vertices_h;
  iEvent.getByLabel(vertexSrc, vertices_h);
  if (!vertices_h.isValid()) {
    std::cout<<"Didja hear the one about the empty vertex collection?\n";
    return;
  }
  
  // require in the event that there is at least one reconstructed vertex
  if(vertices_h->size()<=0) return;

  // pick the first (i.e. highest sum pt) verte
  const reco::Vertex* theVertex=&(vertices_h->front());

  // require that the vertex meets certain criteria
  if(theVertex->ndof()<5) return;
  if(fabs(theVertex->z())>24.0) return;
  if(fabs(theVertex->position().rho())>2.0) return;

  // Get Jet Information

  // magic to get the jets from EDM
  edm::Handle<reco::CaloJetCollection> jets_h;
  iEvent.getByLabel(jetSrc, jets_h);
  if (!jets_h.isValid()) {
    std::cout<<"Didja hear the one about the empty jet collection?\n";
    return;
  }

  // magic to get the jet energy corrections
  const JetCorrector* corrector = JetCorrector::getJetCorrector(jetCorrections,iSetup);
  corrector->vectorialCorrection();
  // collection of selected jets
  std::vector<reco::CaloJet> selectedJets;

  // loop over the jet collection
  for(reco::CaloJetCollection::const_iterator j_it = jets_h->begin(); j_it!=jets_h->end(); j_it++) {
    reco::CaloJet jet = *j_it;

    // put the selected jets into a collection
    selectedJets.push_back(jet);
  }

  hVertexZ->Fill(theVertex->z(), mWeight);
  return;
}


DEFINE_FWK_MODULE(CMSDAS11DijetTestAnalyzer);