L1GctValidation.cc

Go to the documentation of this file.

#include "L1Trigger/GlobalCaloTrigger/plugins/L1GctValidation.h"

#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"

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

#include "CondFormats/L1TObjects/interface/L1GctJetFinderParams.h"
#include "CondFormats/DataRecord/interface/L1GctJetFinderParamsRcd.h"

#include "CondFormats/L1TObjects/interface/L1CaloEtScale.h"
#include "CondFormats/DataRecord/interface/L1HtMissScaleRcd.h"
#include "CondFormats/DataRecord/interface/L1HfRingEtScaleRcd.h"

#include "DataFormats/L1CaloTrigger/interface/L1CaloCollections.h"

#include "DataFormats/Math/interface/deltaPhi.h"

#include <cmath>

L1GctValidation::L1GctValidation(const edm::ParameterSet& iConfig)
    : m_gctinp_tag(iConfig.getUntrackedParameter<edm::InputTag>("rctInputTag", edm::InputTag("rctDigis"))),
      m_energy_tag(iConfig.getUntrackedParameter<edm::InputTag>("gctInputTag", edm::InputTag("gctDigis"))) {}

L1GctValidation::~L1GctValidation() {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}

//
// member functions
//

// ------------ method called to for each event  ------------
void L1GctValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  using namespace edm;

  // Get the scales from the event setup
  ESHandle<L1GctJetFinderParams> jfPars;
  iSetup.get<L1GctJetFinderParamsRcd>().get(jfPars);  // which record?

  double lsbForEt = jfPars.product()->getRgnEtLsbGeV();
  double lsbForHt = jfPars.product()->getHtLsbGeV();

  unsigned httJetThreshold = static_cast<int>(jfPars.product()->getHtJetEtThresholdGeV() / lsbForHt);
  unsigned htmJetThreshold = static_cast<int>(jfPars.product()->getMHtJetEtThresholdGeV() / lsbForHt);

  ESHandle<L1CaloEtScale> htMissScale;
  iSetup.get<L1HtMissScaleRcd>().get(htMissScale);  // which record?
  ESHandle<L1CaloEtScale> hfRingEtScale;
  iSetup.get<L1HfRingEtScaleRcd>().get(hfRingEtScale);  // which record?

  // Get the Gct energy sums from the event
  Handle<L1GctEtTotalCollection> sumEtColl;
  iEvent.getByLabel(m_energy_tag, sumEtColl);
  Handle<L1GctEtHadCollection> sumHtColl;
  iEvent.getByLabel(m_energy_tag, sumHtColl);
  Handle<L1GctEtMissCollection> missEtColl;
  iEvent.getByLabel(m_energy_tag, missEtColl);
  Handle<L1GctHtMissCollection> missHtColl;
  iEvent.getByLabel(m_energy_tag, missHtColl);

  // Get the input calo regions from the event (for checking MEt)
  Handle<L1CaloRegionCollection> inputColl;
  iEvent.getByLabel(m_gctinp_tag, inputColl);

  // Get the internal jet data from the event (for checking Ht)
  Handle<L1GctInternJetDataCollection> internalJetsColl;
  iEvent.getByLabel(m_energy_tag, internalJetsColl);

  double etTot = 0.0;
  for (L1GctEtTotalCollection::const_iterator jbx = sumEtColl->begin(); jbx != sumEtColl->end(); jbx++) {
    if (jbx->bx() == 0) {
      etTot = static_cast<double>(jbx->et());
    }
  }

  double etHad = 0.0;
  for (L1GctEtHadCollection::const_iterator jbx = sumHtColl->begin(); jbx != sumHtColl->end(); jbx++) {
    if (jbx->bx() == 0) {
      etHad = static_cast<double>(jbx->et());
    }
  }

  double etMiss = 0.0;
  double etMAng = 0.0;
  for (L1GctEtMissCollection::const_iterator jbx = missEtColl->begin(); jbx != missEtColl->end(); jbx++) {
    if (jbx->bx() == 0) {
      etMiss = static_cast<double>(jbx->et());
      int phibin = jbx->phi();
      if (phibin >= 36)
        phibin -= 72;
      double etMPhi = static_cast<double>(phibin);

      etMAng = (etMPhi + 0.5) * M_PI / 36.;
    }
  }

  double etTotFromRegions = 0.0;
  double exTotFromRegions = 0.0;
  double eyTotFromRegions = 0.0;
  for (L1CaloRegionCollection::const_iterator jrg = inputColl->begin(); jrg != inputColl->end(); jrg++) {
    if (jrg->bx() == 0) {
      double rgEt = static_cast<double>(jrg->et()) * lsbForEt;
      double rgPhibin = static_cast<double>(jrg->id().iphi());
      double rgPh = (rgPhibin + 0.5) * M_PI / 9.;

      etTotFromRegions += rgEt;
      exTotFromRegions += rgEt * cos(rgPh);
      eyTotFromRegions += rgEt * sin(rgPh);
    }
  }

  double htMissGct = 0.0;
  double htMissAng = 0.0;
  double htMissGeV = 0.0;
  for (L1GctHtMissCollection::const_iterator jbx = missHtColl->begin(); jbx != missHtColl->end(); jbx++) {
    if (jbx->bx() == 0) {
      htMissGct = static_cast<double>(jbx->et());
      htMissGeV = htMissScale->et(jbx->et());
      int phibin = jbx->phi();
      if (phibin >= 9)
        phibin -= 18;
      double htMPhi = static_cast<double>(phibin);
      htMissAng = (htMPhi + 0.5) * M_PI / 9.;
    }
  }

  double htFromJets = 0.0;
  double hxFromJets = 0.0;
  double hyFromJets = 0.0;
  for (L1GctInternJetDataCollection::const_iterator jet = internalJetsColl->begin(); jet != internalJetsColl->end();
       jet++) {
    if (jet->bx() == 0 && !jet->empty()) {
      unsigned jetEtGct = jet->et();
      double jetEt = static_cast<double>(jetEtGct);
      int phibin = jet->regionId().iphi();
      if (phibin >= 9)
        phibin -= 18;
      // The phi bin centres are at 0, 20, 40, ... degrees
      double jetAng = (static_cast<double>(phibin)) * M_PI / 9.;
      if (jetEtGct > httJetThreshold) {
        htFromJets += jetEt;
      }
      if (jetEtGct > htmJetThreshold) {
        hxFromJets += jetEt * cos(jetAng);
        hyFromJets += jetEt * sin(jetAng);
      }
    }
  }

  double dPhiMetMht = deltaPhi(etMAng, htMissAng);

  theSumEtInLsb->Fill(etTot);
  theSumHtInLsb->Fill(etHad);
  theMissEtInLsb->Fill(etMiss);
  theMissHtInLsb->Fill(htMissGct);
  theSumEtInGeV->Fill(etTot * lsbForEt);
  theSumHtInGeV->Fill(etHad * lsbForHt);
  theMissEtInGeV->Fill(etMiss * lsbForEt);
  theMissEtAngle->Fill(etMAng);
  theMissEtVector->Fill(etMiss * lsbForEt * cos(etMAng), etMiss * lsbForEt * sin(etMAng));
  if (htMissGct < 126.5) {
    theMissHtInGeV->Fill(htMissGeV);
    theMissHtAngle->Fill(htMissAng);
    theMissHtVector->Fill(htMissGeV * cos(htMissAng), htMissGeV * sin(htMissAng));
  }

  theSumEtVsInputRegions->Fill(etTot * lsbForEt, etTotFromRegions);
  theMissEtMagVsInputRegions->Fill(etMiss * lsbForEt,
                                   sqrt(exTotFromRegions * exTotFromRegions + eyTotFromRegions * eyTotFromRegions));
  theMissEtAngleVsInputRegions->Fill(etMAng, atan2(-eyTotFromRegions, -exTotFromRegions));
  theMissHtMagVsInputRegions->Fill(htMissGeV,
                                   sqrt(exTotFromRegions * exTotFromRegions + eyTotFromRegions * eyTotFromRegions));

  theMissEtVsMissHt->Fill(etMiss * lsbForEt, htMissGeV);
  theMissEtVsMissHtAngle->Fill(etMAng, htMissAng);
  theDPhiVsMissEt->Fill(dPhiMetMht, etMiss * lsbForEt);
  theDPhiVsMissHt->Fill(dPhiMetMht, htMissGeV);

  theHtVsInternalJetsSum->Fill(etHad * lsbForHt, htFromJets * lsbForHt);
  if (htMissGct < 126.5) {
    theMissHtVsInternalJetsSum->Fill(htMissGeV, sqrt(hxFromJets * hxFromJets + hyFromJets * hyFromJets) * lsbForHt);
    theMissHtPhiVsInternalJetsSum->Fill(htMissAng, atan2(-hyFromJets, -hxFromJets));
    theMissHxVsInternalJetsSum->Fill(htMissGeV * cos(htMissAng), hxFromJets * lsbForHt);
    theMissHyVsInternalJetsSum->Fill(htMissGeV * sin(htMissAng), hyFromJets * lsbForHt);
  }

  // Get minbias trigger quantities from HF
  Handle<L1GctHFRingEtSumsCollection> HFEtSumsColl;
  Handle<L1GctHFBitCountsCollection> HFCountsColl;
  iEvent.getByLabel(m_energy_tag, HFEtSumsColl);
  iEvent.getByLabel(m_energy_tag, HFCountsColl);

  for (L1GctHFRingEtSumsCollection::const_iterator es = HFEtSumsColl->begin(); es != HFEtSumsColl->end(); es++) {
    if (es->bx() == 0) {
      theHfRing0EtSumPositiveEta->Fill(hfRingEtScale->et(es->etSum(0)));
      theHfRing0EtSumNegativeEta->Fill(hfRingEtScale->et(es->etSum(1)));
      theHfRing1EtSumPositiveEta->Fill(hfRingEtScale->et(es->etSum(2)));
      theHfRing1EtSumNegativeEta->Fill(hfRingEtScale->et(es->etSum(3)));
    }
  }

  for (L1GctHFBitCountsCollection::const_iterator bc = HFCountsColl->begin(); bc != HFCountsColl->end(); bc++) {
    if (bc->bx() == 0) {
      theHfRing0CountPositiveEta->Fill(bc->bitCount(0));
      theHfRing0CountNegativeEta->Fill(bc->bitCount(1));
      theHfRing1CountPositiveEta->Fill(bc->bitCount(2));
      theHfRing1CountNegativeEta->Fill(bc->bitCount(3));
    }
  }
}

// ------------ method called once each job just before starting event loop  ------------
void L1GctValidation::beginJob() {
  edm::Service<TFileService> fs;

  TFileDirectory dir0 = fs->mkdir("L1GctEtSums");

  theSumEtInLsb = dir0.make<TH1F>("SumEtInLsb", "Total Et (GCT units)", 128, 0., 2048.);
  theSumHtInLsb = dir0.make<TH1F>("SumHtInLsb", "Total Ht (GCT units)", 128, 0., 2048.);
  theMissEtInLsb = dir0.make<TH1F>("MissEtInLsb", "Missing Et magnitude (GCT units)", 128, 0., 1024.);
  theMissHtInLsb = dir0.make<TH1F>("MissHtInLsb", "Missing Ht magnitude (GCT units)", 128, 0., 127.);
  theSumEtInGeV = dir0.make<TH1F>("SumEtInGeV", "Total Et (in GeV)", 100, 0., 1000.);
  theSumHtInGeV = dir0.make<TH1F>("SumHtInGeV", "Total Ht (in GeV)", 100, 0., 1000.);
  theMissEtInGeV = dir0.make<TH1F>("MissEtInGeV", "Missing Et magnitude (in GeV)", 100, 0., 500.);
  theMissEtAngle = dir0.make<TH1F>("MissEtAngle", "Missing Et angle", 72, -M_PI, M_PI);
  theMissEtVector = dir0.make<TH2F>("MissEtVector", "Missing Ex vs Missing Ey", 100, -100., 100., 100, -100., 100.);
  theMissHtInGeV = dir0.make<TH1F>("MissHtInGeV", "Missing Ht magnitude (in GeV)", 100, 0., 500.);
  theMissHtAngle = dir0.make<TH1F>("MissHtAngle", "Missing Ht angle", 72, -M_PI, M_PI);
  theMissHtVector = dir0.make<TH2F>("MissHtVector", "Missing Hx vs Missing Hy", 100, -100., 100., 100, -100., 100.);
  theSumEtVsInputRegions =
      dir0.make<TH2F>("SumEtVsInputRegions", "Total Et vs sum of input regions", 100, 0., 1000., 100, 0., 1000.);
  theMissEtMagVsInputRegions = dir0.make<TH2F>(
      "MissEtMagVsInputRegions", "Missing Et magnitude vs sum of input regions", 100, 0., 500., 100, 0., 500.);
  theMissEtAngleVsInputRegions = dir0.make<TH2F>(
      "MissEtAngleVsInputRegions", "Missing Et angle vs sum of input regions", 72, -M_PI, M_PI, 72, -M_PI, M_PI);
  theMissHtMagVsInputRegions = dir0.make<TH2F>(
      "MissHtMagVsInputRegions", "Missing Ht magnitude vs sum of input regions", 100, 0., 500., 100, 0., 500.);
  theMissEtVsMissHt = dir0.make<TH2F>("MissEtVsMissHt", "Missing Et vs Missing Ht", 100, 0., 500., 100, 0., 500.);
  theMissEtVsMissHtAngle = dir0.make<TH2F>(
      "MissEtVsMissHtAngle", "Angle correlation Missing Et vs Missing Ht", 72, -M_PI, M_PI, 72, -M_PI, M_PI);
  theDPhiVsMissEt =
      dir0.make<TH2F>("theDPhiVsMissEt", "Angle difference MET-MHT vs MET magnitude", 72, -M_PI, M_PI, 100, 0., 500.);
  theDPhiVsMissHt =
      dir0.make<TH2F>("theDPhiVsMissHt", "Angle difference MET-MHT vs MHT magnitude", 72, -M_PI, M_PI, 100, 0., 500.);

  theHtVsInternalJetsSum = dir0.make<TH2F>(
      "HtVsInternalJetsSum", "Ht vs scalar sum of jet Et values (in GeV)", 128, 0., 2048., 128, 0., 2048.);
  theMissHtVsInternalJetsSum = dir0.make<TH2F>(
      "MissHtVsInternalJetsSum", "Missing Ht vs vector sum of jet Et values (in GeV)", 128, 0., 512., 128, 0., 512.);
  theMissHtPhiVsInternalJetsSum = dir0.make<TH2F>("MissHtPhiVsInternalJetsSum",
                                                  "Angle correlation Missing Ht vs vector sum of jet Et values",
                                                  72,
                                                  -M_PI,
                                                  M_PI,
                                                  72,
                                                  -M_PI,
                                                  M_PI);
  theMissHxVsInternalJetsSum = dir0.make<TH2F>("MissHxVsInternalJetsSum",
                                               "Missing Ht x component vs sum of jet Et values (in GeV)",
                                               128,
                                               -256.,
                                               256.,
                                               128,
                                               -256.,
                                               256.);
  theMissHyVsInternalJetsSum = dir0.make<TH2F>("MissHyVsInternalJetsSum",
                                               "Missing Ht y component vs sum of jet Et values (in GeV)",
                                               128,
                                               -256.,
                                               256.,
                                               128,
                                               -256.,
                                               256.);

  TFileDirectory dir1 = fs->mkdir("L1GctHfSumsAndJetCounts");

  // Minimum bias triggers from Hf inner rings
  theHfRing0EtSumPositiveEta = dir1.make<TH1F>("HfRing0EtSumPositiveEta", "Hf Inner Ring0 Et eta+", 60, 0., 30.);
  theHfRing0EtSumNegativeEta = dir1.make<TH1F>("HfRing0EtSumNegativeEta", "Hf Inner Ring0 Et eta-", 60, 0., 30.);
  theHfRing1EtSumPositiveEta = dir1.make<TH1F>("HfRing1EtSumPositiveEta", "Hf Inner Ring1 Et eta+", 60, 0., 30.);
  theHfRing1EtSumNegativeEta = dir1.make<TH1F>("HfRing1EtSumNegativeEta", "Hf Inner Ring1 Et eta-", 60, 0., 30.);
  theHfRing0CountPositiveEta = dir1.make<TH1F>("HfRing0CountPositiveEta", "Hf Threshold bits Ring0 eta+", 20, 0., 20.);
  theHfRing0CountNegativeEta = dir1.make<TH1F>("HfRing0CountNegativeEta", "Hf Threshold bits Ring0 eta-", 20, 0., 20.);
  theHfRing1CountPositiveEta = dir1.make<TH1F>("HfRing1CountPositiveEta", "Hf Threshold bits Ring1 eta+", 20, 0., 20.);
  theHfRing1CountNegativeEta = dir1.make<TH1F>("HfRing1CountNegativeEta", "Hf Threshold bits Ring1 eta-", 20, 0., 20.);
}

// ------------ method called once each job just after ending the event loop  ------------
void L1GctValidation::endJob() {}

DEFINE_FWK_MODULE(L1GctValidation);