Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-13 13:14:13

0001 // CMSDAS11DijetAnalyzer.cc
0002 // Description: A basic dijet analyzer for the CMSDAS 2011
0003 // Author: John Paul Chou
0004 // Date: January 12, 2011
0005 
0006 #include "RecoJets/JetAnalyzers/interface/CMSDAS11DijetAnalyzer.h"
0007 
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0012 
0013 #include "DataFormats/VertexReco/interface/Vertex.h"
0014 #include "JetMETCorrections/Objects/interface/JetCorrector.h"
0015 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0016 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0017 
0018 #include <TH1D.h>
0019 
0020 CMSDAS11DijetAnalyzer::CMSDAS11DijetAnalyzer(edm::ParameterSet const& params)
0021     : jetSrc(params.getParameter<edm::InputTag>("jetSrc")),
0022       vertexSrc(params.getParameter<edm::InputTag>("vertexSrc")),
0023       jetCorrections(params.getParameter<std::string>("jetCorrections")),
0024       innerDeltaEta(params.getParameter<double>("innerDeltaEta")),
0025       outerDeltaEta(params.getParameter<double>("outerDeltaEta")),
0026       JESbias(params.getParameter<double>("JESbias")) {
0027   // setup file service
0028   usesResource(TFileService::kSharedResource);
0029   edm::Service<TFileService> fs;
0030 
0031   const int NBINS = 36;
0032   Double_t BOUNDARIES[NBINS] = {220,  244,  270,  296,  325,  354,  386,  419,  453,  489,  526,  565,
0033                                 606,  649,  693,  740,  788,  838,  890,  944,  1000, 1058, 1118, 1181,
0034                                 1246, 1313, 1383, 1455, 1530, 1607, 1687, 1770, 1856, 1945, 2037, 2132};
0035 
0036   // setup histograms
0037   hVertexZ = fs->make<TH1D>("hVertexZ", "Z position of the Vertex", 50, -20, 20);
0038   hJetRawPt = fs->make<TH1D>("hJetRawPt", "Raw Jet Pt", 50, 0, 1000);
0039   hJetCorrPt = fs->make<TH1D>("hJetCorrPt", "Corrected Jet Pt", 50, 0, 1000);
0040   hJet1Pt = fs->make<TH1D>("hJet1Pt", "Corrected Jet1 Pt", 50, 0, 1000);
0041   hJet2Pt = fs->make<TH1D>("hJet2Pt", "Corrected Jet2 Pt", 50, 0, 1000);
0042 
0043   hJetEta = fs->make<TH1D>("hJetEta", "Corrected Jet Eta", 10, -5, 5);
0044   hJet1Eta = fs->make<TH1D>("hJet1Eta", "Corrected Jet1 Eta", 10, -5, 5);
0045   hJet2Eta = fs->make<TH1D>("hJet2Eta", "Corrected Jet2 Eta", 10, -5, 5);
0046 
0047   hJetPhi = fs->make<TH1D>("hJetPhi", "Corrected Jet Phi", 10, -3.1415, 3.1415);
0048   hJet1Phi = fs->make<TH1D>("hJet1Phi", "Corrected Jet1 Phi", 10, -3.1415, 3.1415);
0049   hJet2Phi = fs->make<TH1D>("hJet2Phi", "Corrected Jet2 Phi", 10, -3.1415, 3.1415);
0050 
0051   hJetEMF = fs->make<TH1D>("hJetEMF", "EM Fraction of Jets", 50, 0, 1);
0052   hJet1EMF = fs->make<TH1D>("hJet1EMF", "EM Fraction of Jet1", 50, 0, 1);
0053   hJet2EMF = fs->make<TH1D>("hJet2EMF", "EM Fraction of Jet2", 50, 0, 1);
0054 
0055   hCorDijetMass = fs->make<TH1D>("hCorDijetMass", "Corrected Dijet Mass", NBINS - 1, BOUNDARIES);
0056   hDijetDeltaPhi = fs->make<TH1D>("hDijetDeltaPhi", "Dijet |#Delta #phi|", 50, 0, 3.1415);
0057   hDijetDeltaEta = fs->make<TH1D>("hDijetDeltaEta", "Dijet |#Delta #eta|", 50, 0, 6);
0058 
0059   hInnerDijetMass = fs->make<TH1D>("hInnerDijetMass", "Corrected Inner Dijet Mass", NBINS - 1, BOUNDARIES);
0060   hOuterDijetMass = fs->make<TH1D>("hOuterDijetMass", "Corrected Outer Dijet Mass", NBINS - 1, BOUNDARIES);
0061 }
0062 
0063 void CMSDAS11DijetAnalyzer::endJob(void) {}
0064 
0065 void CMSDAS11DijetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0066   ////////Get Weight, if this is MC//////////////
0067   double mWeight;
0068   edm::Handle<GenEventInfoProduct> hEventInfo;
0069   iEvent.getByLabel("generator", hEventInfo);
0070   if (hEventInfo.isValid())
0071     mWeight = hEventInfo->weight();
0072   else
0073     mWeight = 1.;
0074 
0075   ////////////////////////////////////////////
0076   // Get event ID information
0077   ////////////////////////////////////////////
0078   //int nrun=iEvent.id().run();
0079   //int nlumi=iEvent.luminosityBlock();
0080   //int nevent=iEvent.id().event();
0081 
0082   ////////////////////////////////////////////
0083   // Get Primary Vertex Information
0084   ////////////////////////////////////////////
0085 
0086   // magic to get the vertices from EDM
0087   edm::Handle<std::vector<reco::Vertex> > vertices_h;
0088   iEvent.getByLabel(vertexSrc, vertices_h);
0089   if (!vertices_h.isValid()) {
0090     std::cout << "Didja hear the one about the empty vertex collection?\n";
0091     return;
0092   }
0093 
0094   // require in the event that there is at least one reconstructed vertex
0095   if (vertices_h->empty())
0096     return;
0097 
0098   // pick the first (i.e. highest sum pt) verte
0099   const reco::Vertex* theVertex = &(vertices_h->front());
0100 
0101   // require that the vertex meets certain criteria
0102   if (theVertex->ndof() < 5)
0103     return;
0104   if (fabs(theVertex->z()) > 24.0)
0105     return;
0106   if (fabs(theVertex->position().rho()) > 2.0)
0107     return;
0108 
0109   ////////////////////////////////////////////
0110   // Get Jet Information
0111   ////////////////////////////////////////////
0112 
0113   // magic to get the jets from EDM
0114   edm::Handle<reco::CaloJetCollection> jets_h;
0115   iEvent.getByLabel(jetSrc, jets_h);
0116   if (!jets_h.isValid()) {
0117     std::cout << "Didja hear the one about the empty jet collection?\n";
0118     return;
0119   }
0120 
0121   // magic to get the jet energy corrections
0122   const JetCorrector* corrector = JetCorrector::getJetCorrector(jetCorrections, iSetup);
0123 
0124   // collection of selected jets
0125   std::vector<reco::CaloJet> selectedJets;
0126 
0127   // loop over the jet collection
0128   for (reco::CaloJetCollection::const_iterator j_it = jets_h->begin(); j_it != jets_h->end(); j_it++) {
0129     reco::CaloJet jet = *j_it;
0130 
0131     // calculate and apply the correction
0132     double scale = corrector->correction(jet.p4());
0133 
0134     // Introduce a purposeful bias to the correction, to show what happens
0135     scale *= JESbias;
0136 
0137     // fill the histograms
0138     hJetRawPt->Fill(jet.pt(), mWeight);
0139     hJetEta->Fill(jet.eta(), mWeight);
0140     hJetPhi->Fill(jet.phi(), mWeight);
0141     hJetEMF->Fill(jet.emEnergyFraction(), mWeight);
0142 
0143     // correct the jet energy
0144     jet.scaleEnergy(scale);
0145 
0146     // now fill the correct jet pt after the correction
0147     hJetCorrPt->Fill(jet.pt(), mWeight);
0148 
0149     // put the selected jets into a collection
0150     selectedJets.push_back(jet);
0151   }
0152 
0153   // require at least two jets to continue
0154   if (selectedJets.size() < 2)
0155     return;
0156 
0157   //sort by corrected pt (not the same order as raw pt, sometimes)
0158   sort(selectedJets.begin(), selectedJets.end(), compare_JetPt);
0159 
0160   // select high pt, central, non-noise-like jets
0161   if (selectedJets[0].pt() < 50.0)
0162     return;
0163   if (fabs(selectedJets[0].eta()) > 2.5)
0164     return;
0165   if (selectedJets[0].emEnergyFraction() < 0.01)
0166     return;
0167   if (selectedJets[1].pt() < 50.0)
0168     return;
0169   if (fabs(selectedJets[1].eta()) > 2.5)
0170     return;
0171   if (selectedJets[1].emEnergyFraction() < 0.01)
0172     return;
0173 
0174   // fill histograms for the jets in our dijets, only
0175   hJet1Pt->Fill(selectedJets[0].pt(), mWeight);
0176   hJet1Eta->Fill(selectedJets[0].eta(), mWeight);
0177   hJet1Phi->Fill(selectedJets[0].phi(), mWeight);
0178   hJet1EMF->Fill(selectedJets[0].emEnergyFraction(), mWeight);
0179   hJet2Pt->Fill(selectedJets[1].pt(), mWeight);
0180   hJet2Eta->Fill(selectedJets[1].eta(), mWeight);
0181   hJet2Phi->Fill(selectedJets[1].phi(), mWeight);
0182   hJet2EMF->Fill(selectedJets[1].emEnergyFraction(), mWeight);
0183 
0184   //Get the mass of the two leading jets.  Needs their 4-vectors...
0185   double corMass = (selectedJets[0].p4() + selectedJets[1].p4()).M();
0186   double deltaEta = fabs(selectedJets[0].eta() - selectedJets[1].eta());
0187   if (corMass < 489)
0188     return;
0189   if (deltaEta > 1.3)
0190     return;
0191   hCorDijetMass->Fill(corMass, mWeight);
0192   hDijetDeltaPhi->Fill(fabs(selectedJets[0].phi() - selectedJets[1].phi()), mWeight);
0193   hDijetDeltaEta->Fill(deltaEta, mWeight);
0194 
0195   //Fill the inner and outer dijet mass spectra, to make the ratio from
0196   if (deltaEta < innerDeltaEta)
0197     hInnerDijetMass->Fill(corMass, mWeight);
0198   else if (deltaEta < outerDeltaEta)
0199     hOuterDijetMass->Fill(corMass, mWeight);
0200 
0201   hVertexZ->Fill(theVertex->z(), mWeight);
0202   return;
0203 }
0204 
0205 DEFINE_FWK_MODULE(CMSDAS11DijetAnalyzer);