DijetMass.cc

Go to the documentation of this file.

// DijetMass.cc
// Description:  Some Basic validation plots for jets.
// Author: Robert M. Harris
// Date:  30 - August - 2006
// Kalanand Mishra (November 22, 2009):
//          Modified and cleaned up to work in 3.3.X
//
#include "RecoJets/JetAnalyzers/interface/DijetMass.h"

#include "DataFormats/JetReco/interface/CaloJetCollection.h"
#include "DataFormats/JetReco/interface/PFJetCollection.h"
#include "DataFormats/JetReco/interface/GenJetCollection.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/JetReco/interface/PFJet.h"
#include "DataFormats/JetReco/interface/GenJet.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include <TROOT.h>
#include <TSystem.h>
#include <TFile.h>
#include <TCanvas.h>
#include <cmath>
using namespace edm;
using namespace reco;
using namespace std;

template <class Jet>
DijetMass<Jet>::DijetMass(const edm::ParameterSet& cfg) {
  PtHistMax = cfg.getUntrackedParameter<double>("PtHistMax", 3000.0);
  EtaMax = cfg.getUntrackedParameter<double>("EtaMax", 1.3);
  histogramFile = cfg.getUntrackedParameter<std::string>("HistoFileName", "DijetMassHistos.root");

  AKJets = cfg.getParameter<std::string>("AKJets");
  AKCorJets = cfg.getParameter<std::string>("AKCorrectedJets");
  ICJets = cfg.getParameter<std::string>("ICJets");
  ICCorJets = cfg.getParameter<std::string>("ICCorrectedJets");
  SCJets = cfg.getParameter<std::string>("SCJets");
  SCCorJets = cfg.getParameter<std::string>("SCCorrectedJets");
  KTJets = cfg.getParameter<std::string>("KTJets");
  KTCorJets = cfg.getParameter<std::string>("KTCorrectedJets");
}

template <class Jet>
void DijetMass<Jet>::beginJob() {
  cout << "DijetMass: Maximum bin edge for Pt Hists = " << PtHistMax << endl;
  numJets = 2;

  //Initialize some stuff
  evtCount = 0;

  // Open the histogram file and book some associated histograms
  m_file = new TFile(histogramFile.c_str(), "RECREATE");

  //Simple histos

  //AK unc
  ptAKunc = TH1F("ptAKunc", "p_{T} of leading Jets (AK)", 50, 0.0, PtHistMax);
  etaAKunc = TH1F("etaAKunc", "#eta of leading Jets (AK)", 23, -1.0, 1.0);
  phiAKunc = TH1F("phiAKunc", "#phi of leading Jets (AK)", 72, -M_PI, M_PI);
  m2jAKunc = TH1F("m2jAKunc", "Dijet Mass of leading Jets (AK)", 100, 0.0, 2 * PtHistMax);

  //AK cor
  ptAKcor = TH1F("ptAKcor", "p_{T} of leading Corrected Jets (AK)", 50, 0.0, PtHistMax);
  etaAKcor = TH1F("etaAKcor", "#eta of leading Corrected Jets (AK)", 23, -1.0, 1.0);
  phiAKcor = TH1F("phiAKcor", "#phi of leading Corrected Jets (AK)", 72, -M_PI, M_PI);
  m2jAKcor = TH1F("m2jAKcor", "Dijet Mass of leading Corrected Jets (AK)", 100, 0.0, 2 * PtHistMax);

  //IC unc
  ptICunc = TH1F("ptICunc", "p_{T} of leading Jets (IC)", 50, 0.0, PtHistMax);
  etaICunc = TH1F("etaICunc", "#eta of leading Jets (IC)", 23, -1.0, 1.0);
  phiICunc = TH1F("phiICunc", "#phi of leading Jets (IC)", 72, -M_PI, M_PI);
  m2jICunc = TH1F("m2jICunc", "Dijet Mass of leading Jets (IC)", 100, 0.0, 2 * PtHistMax);

  //IC cor
  ptICcor = TH1F("ptICcor", "p_{T} of leading Corrected Jets (IC)", 50, 0.0, PtHistMax);
  etaICcor = TH1F("etaICcor", "#eta of leading Corrected Jets (IC)", 23, -1.0, 1.0);
  phiICcor = TH1F("phiICcor", "#phi of leading Corrected Jets (IC)", 72, -M_PI, M_PI);
  m2jICcor = TH1F("m2jICcor", "Dijet Mass of leading Corrected Jets (IC)", 100, 0.0, 2 * PtHistMax);

  //KT unc
  ptKTunc = TH1F("ptKTunc", "p_{T} of leading Jets (KT)", 50, 0.0, PtHistMax);
  etaKTunc = TH1F("etaKTunc", "#eta of leading Jets (KT)", 23, -1.0, 1.0);
  phiKTunc = TH1F("phiKTunc", "#phi of leading Jets (KT)", 72, -M_PI, M_PI);
  m2jKTunc = TH1F("m2jKTunc", "Dijet Mass of leading Jets (KT)", 100, 0.0, 2 * PtHistMax);

  //KT cor
  ptKTcor = TH1F("ptKTcor", "p_{T} of leading Corrected Jets (KT)", 50, 0.0, PtHistMax);
  etaKTcor = TH1F("etaKTcor", "#eta of leading Corrected Jets (KT)", 23, -1.0, 1.0);
  phiKTcor = TH1F("phiKTcor", "#phi of leading Corrected Jets (KT)", 72, -M_PI, M_PI);
  m2jKTcor = TH1F("m2jKTcor", "Dijet Mass of leading Corrected Jets (KT)", 100, 0.0, 2 * PtHistMax);

  //SC unc
  ptSCunc = TH1F("ptSCunc", "p_{T} of leading Jets (SC)", 50, 0.0, PtHistMax);
  etaSCunc = TH1F("etaSCunc", "#eta of leading Jets (SC)", 23, -1.0, 1.0);
  phiSCunc = TH1F("phiSCunc", "#phi of leading Jets (SC)", 72, -M_PI, M_PI);
  m2jSCunc = TH1F("m2jSCunc", "Dijet Mass of leading Jets (SC)", 100, 0.0, 2 * PtHistMax);

  //SC cor
  ptSCcor = TH1F("ptSCcor", "p_{T} of leading Corrected Jets (SC)", 50, 0.0, PtHistMax);
  etaSCcor = TH1F("etaSCcor", "#eta of leading Corrected Jets (SC)", 23, -1.0, 1.0);
  phiSCcor = TH1F("phiSCcor", "#phi of leading Corrected Jets (SC)", 72, -M_PI, M_PI);
  m2jSCcor = TH1F("m2jSCcor", "Dijet Mass of leading Corrected Jets (SC)", 100, 0.0, 2 * PtHistMax);
}

template <class Jet>
void DijetMass<Jet>::analyze(const Event& evt, const EventSetup& es) {
  evtCount++;
  math::XYZTLorentzVector p4jet[2];
  int jetInd;
  Handle<JetCollection> Jets;

  //Fill Simple Histos
  typename JetCollection::const_iterator i_jet;

  //AK unc
  evt.getByLabel(AKJets, Jets);
  jetInd = 0;
  for (i_jet = Jets->begin(); i_jet != Jets->end() && jetInd < 2; ++i_jet) {
    p4jet[jetInd] = i_jet->p4();
    jetInd++;
  }
  if (jetInd == 2 && abs(p4jet[0].eta()) < EtaMax && abs(p4jet[1].eta()) < EtaMax) {
    m2jAKunc.Fill((p4jet[0] + p4jet[1]).mass());
    ptAKunc.Fill(p4jet[0].Pt());
    ptAKunc.Fill(p4jet[1].Pt());
    etaAKunc.Fill(p4jet[0].eta());
    etaAKunc.Fill(p4jet[1].eta());
    phiAKunc.Fill(p4jet[0].phi());
    phiAKunc.Fill(p4jet[1].phi());
  }

  //AK corrected
  evt.getByLabel(AKCorJets, Jets);
  jetInd = 0;
  for (i_jet = Jets->begin(); i_jet != Jets->end() && jetInd < 2; ++i_jet) {
    p4jet[jetInd] = i_jet->p4();
    jetInd++;
  }
  if (jetInd == 2 && abs(p4jet[0].eta()) < EtaMax && abs(p4jet[1].eta()) < EtaMax) {
    m2jAKcor.Fill((p4jet[0] + p4jet[1]).mass());
    ptAKcor.Fill(p4jet[0].Pt());
    ptAKcor.Fill(p4jet[1].Pt());
    etaAKcor.Fill(p4jet[0].eta());
    etaAKcor.Fill(p4jet[1].eta());
    phiAKcor.Fill(p4jet[0].phi());
    phiAKcor.Fill(p4jet[1].phi());
  }

  //IC unc
  evt.getByLabel(ICJets, Jets);
  jetInd = 0;
  for (i_jet = Jets->begin(); i_jet != Jets->end() && jetInd < 2; ++i_jet) {
    p4jet[jetInd] = i_jet->p4();
    jetInd++;
  }
  if (jetInd == 2 && abs(p4jet[0].eta()) < EtaMax && abs(p4jet[1].eta()) < EtaMax) {
    m2jICunc.Fill((p4jet[0] + p4jet[1]).mass());
    ptICunc.Fill(p4jet[0].Pt());
    ptICunc.Fill(p4jet[1].Pt());
    etaICunc.Fill(p4jet[0].eta());
    etaICunc.Fill(p4jet[1].eta());
    phiICunc.Fill(p4jet[0].phi());
    phiICunc.Fill(p4jet[1].phi());
  }

  //IC corrected
  evt.getByLabel(ICCorJets, Jets);
  jetInd = 0;
  for (i_jet = Jets->begin(); i_jet != Jets->end() && jetInd < 2; ++i_jet) {
    p4jet[jetInd] = i_jet->p4();
    jetInd++;
  }
  if (jetInd == 2 && abs(p4jet[0].eta()) < EtaMax && abs(p4jet[1].eta()) < EtaMax) {
    m2jICcor.Fill((p4jet[0] + p4jet[1]).mass());
    ptICcor.Fill(p4jet[0].Pt());
    ptICcor.Fill(p4jet[1].Pt());
    etaICcor.Fill(p4jet[0].eta());
    etaICcor.Fill(p4jet[1].eta());
    phiICcor.Fill(p4jet[0].phi());
    phiICcor.Fill(p4jet[1].phi());
  }

  //KT unc
  evt.getByLabel(KTJets, Jets);
  jetInd = 0;
  for (i_jet = Jets->begin(); i_jet != Jets->end() && jetInd < 2; ++i_jet) {
    p4jet[jetInd] = i_jet->p4();
    jetInd++;
  }
  if (jetInd == 2 && abs(p4jet[0].eta()) < EtaMax && abs(p4jet[1].eta()) < EtaMax) {
    m2jKTunc.Fill((p4jet[0] + p4jet[1]).mass());
    ptKTunc.Fill(p4jet[0].Pt());
    ptKTunc.Fill(p4jet[1].Pt());
    etaKTunc.Fill(p4jet[0].eta());
    etaKTunc.Fill(p4jet[1].eta());
    phiKTunc.Fill(p4jet[0].phi());
    phiKTunc.Fill(p4jet[1].phi());
  }

  //KT corrected
  evt.getByLabel(KTCorJets, Jets);
  jetInd = 0;
  for (i_jet = Jets->begin(); i_jet != Jets->end() && jetInd < 2; ++i_jet) {
    p4jet[jetInd] = i_jet->p4();
    jetInd++;
  }
  if (jetInd == 2 && abs(p4jet[0].eta()) < EtaMax && abs(p4jet[1].eta()) < EtaMax) {
    m2jKTcor.Fill((p4jet[0] + p4jet[1]).mass());
    ptKTcor.Fill(p4jet[0].Pt());
    ptKTcor.Fill(p4jet[1].Pt());
    etaKTcor.Fill(p4jet[0].eta());
    etaKTcor.Fill(p4jet[1].eta());
    phiKTcor.Fill(p4jet[0].phi());
    phiKTcor.Fill(p4jet[1].phi());
  }

  //SC unc
  evt.getByLabel(SCJets, Jets);
  jetInd = 0;
  for (i_jet = Jets->begin(); i_jet != Jets->end() && jetInd < 2; ++i_jet) {
    p4jet[jetInd] = i_jet->p4();
    jetInd++;
  }
  if (jetInd == 2 && abs(p4jet[0].eta()) < EtaMax && abs(p4jet[1].eta()) < EtaMax) {
    m2jSCunc.Fill((p4jet[0] + p4jet[1]).mass());
    ptSCunc.Fill(p4jet[0].Pt());
    ptSCunc.Fill(p4jet[1].Pt());
    etaSCunc.Fill(p4jet[0].eta());
    etaSCunc.Fill(p4jet[1].eta());
    phiSCunc.Fill(p4jet[0].phi());
    phiSCunc.Fill(p4jet[1].phi());
  }

  //SC corrected
  evt.getByLabel(SCCorJets, Jets);
  jetInd = 0;
  for (i_jet = Jets->begin(); i_jet != Jets->end() && jetInd < 2; ++i_jet) {
    p4jet[jetInd] = i_jet->p4();
    jetInd++;
  }
  if (jetInd == 2 && abs(p4jet[0].eta()) < EtaMax && abs(p4jet[1].eta()) < EtaMax) {
    m2jSCcor.Fill((p4jet[0] + p4jet[1]).mass());
    ptSCcor.Fill(p4jet[0].Pt());
    ptSCcor.Fill(p4jet[1].Pt());
    etaSCcor.Fill(p4jet[0].eta());
    etaSCcor.Fill(p4jet[1].eta());
    phiSCcor.Fill(p4jet[0].phi());
    phiSCcor.Fill(p4jet[1].phi());
  }
}

template <class Jet>
void DijetMass<Jet>::endJob() {
  //Write out the histogram file.
  m_file->cd();

  ptAKunc.Write();
  etaAKunc.Write();
  phiAKunc.Write();
  m2jAKunc.Write();

  ptAKcor.Write();
  etaAKcor.Write();
  phiAKcor.Write();
  m2jAKcor.Write();

  ptICunc.Write();
  etaICunc.Write();
  phiICunc.Write();
  m2jICunc.Write();

  ptICcor.Write();
  etaICcor.Write();
  phiICcor.Write();
  m2jICcor.Write();

  ptKTunc.Write();
  etaKTunc.Write();
  phiKTunc.Write();
  m2jKTunc.Write();

  ptKTcor.Write();
  etaKTcor.Write();
  phiKTcor.Write();
  m2jKTcor.Write();

  ptSCunc.Write();
  etaSCunc.Write();
  phiSCunc.Write();
  m2jSCunc.Write();

  ptSCcor.Write();
  etaSCcor.Write();
  phiSCcor.Write();
  m2jSCcor.Write();

  m_file->Close();
}
#include "FWCore/Framework/interface/MakerMacros.h"
typedef DijetMass<CaloJet> DijetMassCaloJets;
DEFINE_FWK_MODULE(DijetMassCaloJets);
typedef DijetMass<GenJet> DijetMassGenJets;
DEFINE_FWK_MODULE(DijetMassGenJets);
typedef DijetMass<PFJet> DijetMassPFJets;
DEFINE_FWK_MODULE(DijetMassPFJets);