DQMHcalPhiSymAlCaReco.cc

Go to the documentation of this file.

/*
 * \file DQMHcalPhiSymAlCaReco.cc
 *
 * \author Olga Kodolova
 *
 *
 *
 * Description: Monitoring of Phi Symmetry Calibration Stream
 */

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

// DQM include files

#include "DQMServices/Core/interface/MonitorElement.h"

// work on collections

#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"

#include "DataFormats/FEDRawData/interface/FEDHeader.h"
#include "DataFormats/FEDRawData/interface/FEDNumbering.h"
#include "DataFormats/FEDRawData/interface/FEDRawData.h"
#include "EventFilter/HcalRawToDigi/interface/HcalDCCHeader.h"
#include "EventFilter/HcalRawToDigi/interface/HcalHTRData.h"

#include "DQMOffline/CalibCalo/src/DQMHcalPhiSymAlCaReco.h"
#include "DQMServices/Core/interface/DQMStore.h"

using namespace std;
using namespace edm;

// ******************************************
// constructors
// *****************************************

DQMHcalPhiSymAlCaReco::DQMHcalPhiSymAlCaReco(const edm::ParameterSet &ps) : eventCounter_(0) {
  //
  // Input from configurator file
  //
  folderName_ = ps.getUntrackedParameter<string>("FolderName", "ALCAStreamHcalPhiSym");

  hbherecoMB = consumes<HBHERecHitCollection>(ps.getParameter<edm::InputTag>("hbheInputMB"));
  horecoMB = ps.getParameter<edm::InputTag>("hoInputMB");
  hfrecoMB = consumes<HFRecHitCollection>(ps.getParameter<edm::InputTag>("hfInputMB"));

  hbherecoNoise = consumes<HBHERecHitCollection>(ps.getParameter<edm::InputTag>("hbheInputNoise"));
  horecoNoise = ps.getParameter<edm::InputTag>("hoInputNoise");
  hfrecoNoise = consumes<HFRecHitCollection>(ps.getParameter<edm::InputTag>("hfInputNoise"));

  rawInLabel_ = consumes<FEDRawDataCollection>(ps.getParameter<edm::InputTag>("rawInputLabel"));

  period_ = ps.getParameter<unsigned int>("period");

  saveToFile_ = ps.getUntrackedParameter<bool>("SaveToFile", false);
  fileName_ = ps.getUntrackedParameter<string>("FileName", "MonitorAlCaHcalPhiSym.root");

  // histogram parameters

  // Distribution of rechits in iPhi, iEta
  hiDistr_y_nbin_ = ps.getUntrackedParameter<int>("hiDistr_y_nbin", 72);
  hiDistr_y_min_ = ps.getUntrackedParameter<double>("hiDistr_y_min", 0.5);
  hiDistr_y_max_ = ps.getUntrackedParameter<double>("hiDistr_y_max", 72.5);
  hiDistr_x_nbin_ = ps.getUntrackedParameter<int>("hiDistr_x_nbin", 41);
  hiDistr_x_min_ = ps.getUntrackedParameter<double>("hiDistr_x_min", 0.5);
  hiDistr_x_max_ = ps.getUntrackedParameter<double>("hiDistr_x_max", 41.5);
  // Check for NZS
  hiDistr_r_nbin_ = ps.getUntrackedParameter<int>("hiDistr_r_nbin", 100);
  ihbhe_size_ = ps.getUntrackedParameter<double>("ihbhe_size_", 5184.);
  ihf_size_ = ps.getUntrackedParameter<double>("ihf_size_", 1728.);
}

DQMHcalPhiSymAlCaReco::~DQMHcalPhiSymAlCaReco() {}

//--------------------------------------------------------
void DQMHcalPhiSymAlCaReco::bookHistograms(DQMStore::IBooker &ibooker,
                                           edm::Run const &irun,
                                           edm::EventSetup const &isetup) {
  // create and cd into new folder
  ibooker.setCurrentFolder(folderName_);

  eventCounter_ = 0;

  hFEDsize = ibooker.book1D("hFEDsize", "HCAL FED size (kB)", 200, -0.5, 20.5);
  hFEDsize->setAxisTitle("kB", 1);

  hHcalIsZS = ibooker.book1D("hHcalIsZS", "Hcal Is ZS", 4, -1.5, 2.5);
  hHcalIsZS->setBinLabel(2, "NZS");
  hHcalIsZS->setBinLabel(3, "ZS");

  char hname[50];
  sprintf(hname, "L1 Event Number %% %i", period_);
  hL1Id = ibooker.book1D("hL1Id", hname, 4200, -99.5, 4099.5);
  hL1Id->setAxisTitle(hname);

  // book some histograms 1D
  double xmin = 0.1;
  double xmax = 1.1;
  hiDistrHBHEsize1D_ = ibooker.book1D("DistrHBHEsize", "Size of HBHE Collection", hiDistr_r_nbin_, xmin, xmax);
  hiDistrHFsize1D_ = ibooker.book1D("DistrHFsize", "Size of HF Collection", hiDistr_r_nbin_, xmin, xmax);

  // First moment
  hiDistrMBPl2D_ = ibooker.book2D("MBdepthPl1",
                                  "iphi- +ieta signal distribution at depth1",
                                  hiDistr_x_nbin_,
                                  hiDistr_x_min_,
                                  hiDistr_x_max_,
                                  hiDistr_y_nbin_,
                                  hiDistr_y_min_,
                                  hiDistr_y_max_);

  hiDistrMBPl2D_->setAxisTitle("i#phi ", 2);
  hiDistrMBPl2D_->setAxisTitle("i#eta ", 1);

  hiDistrNoisePl2D_ = ibooker.book2D("NoisedepthPl1",
                                     "iphi-ieta noise distribution at depth1",
                                     hiDistr_x_nbin_ + 1,
                                     hiDistr_x_min_ - 1.,
                                     hiDistr_x_max_,
                                     hiDistr_y_nbin_ + 1,
                                     hiDistr_y_min_ - 1.,
                                     hiDistr_y_max_);

  hiDistrNoisePl2D_->setAxisTitle("i#phi ", 2);
  hiDistrNoisePl2D_->setAxisTitle("i#eta ", 1);
  // Second moment
  hiDistrMB2Pl2D_ = ibooker.book2D("MB2depthPl1",
                                   "iphi- +ieta signal distribution at depth1",
                                   hiDistr_x_nbin_,
                                   hiDistr_x_min_,
                                   hiDistr_x_max_,
                                   hiDistr_y_nbin_,
                                   hiDistr_y_min_,
                                   hiDistr_y_max_);

  hiDistrMB2Pl2D_->setAxisTitle("i#phi ", 2);
  hiDistrMB2Pl2D_->setAxisTitle("i#eta ", 1);

  hiDistrNoise2Pl2D_ = ibooker.book2D("Noise2depthPl1",
                                      "iphi-ieta noise distribution at depth1",
                                      hiDistr_x_nbin_,
                                      hiDistr_x_min_,
                                      hiDistr_x_max_,
                                      hiDistr_y_nbin_,
                                      hiDistr_y_min_,
                                      hiDistr_y_max_);

  hiDistrNoise2Pl2D_->setAxisTitle("i#phi ", 2);
  hiDistrNoise2Pl2D_->setAxisTitle("i#eta ", 1);

  // Variance
  hiDistrVarMBPl2D_ = ibooker.book2D("VarMBdepthPl1",
                                     "iphi- +ieta signal distribution at depth1",
                                     hiDistr_x_nbin_,
                                     hiDistr_x_min_,
                                     hiDistr_x_max_,
                                     hiDistr_y_nbin_,
                                     hiDistr_y_min_,
                                     hiDistr_y_max_);

  hiDistrVarMBPl2D_->setAxisTitle("i#phi ", 2);
  hiDistrVarMBPl2D_->setAxisTitle("i#eta ", 1);

  hiDistrVarNoisePl2D_ = ibooker.book2D("VarNoisedepthPl1",
                                        "iphi-ieta noise distribution at depth1",
                                        hiDistr_x_nbin_,
                                        hiDistr_x_min_,
                                        hiDistr_x_max_,
                                        hiDistr_y_nbin_,
                                        hiDistr_y_min_,
                                        hiDistr_y_max_);

  hiDistrVarNoisePl2D_->setAxisTitle("i#phi ", 2);
  hiDistrVarNoisePl2D_->setAxisTitle("i#eta ", 1);

  //==================================================================================
  // First moment
  hiDistrMBMin2D_ = ibooker.book2D("MBdepthMin1",
                                   "iphi- +ieta signal distribution at depth1",
                                   hiDistr_x_nbin_,
                                   hiDistr_x_min_,
                                   hiDistr_x_max_,
                                   hiDistr_y_nbin_,
                                   hiDistr_y_min_,
                                   hiDistr_y_max_);

  hiDistrMBMin2D_->setAxisTitle("i#phi ", 2);
  hiDistrMBMin2D_->setAxisTitle("i#eta ", 1);

  hiDistrNoiseMin2D_ = ibooker.book2D("NoisedepthMin1",
                                      "iphi-ieta noise distribution at depth1",
                                      hiDistr_x_nbin_,
                                      hiDistr_x_min_,
                                      hiDistr_x_max_,
                                      hiDistr_y_nbin_,
                                      hiDistr_y_min_,
                                      hiDistr_y_max_);

  hiDistrNoiseMin2D_->setAxisTitle("i#phi ", 2);
  hiDistrNoiseMin2D_->setAxisTitle("i#eta ", 1);
  // Second moment
  hiDistrMB2Min2D_ = ibooker.book2D("MB2depthMin1",
                                    "iphi- +ieta signal distribution at depth1",
                                    hiDistr_x_nbin_,
                                    hiDistr_x_min_,
                                    hiDistr_x_max_,
                                    hiDistr_y_nbin_,
                                    hiDistr_y_min_,
                                    hiDistr_y_max_);

  hiDistrMB2Min2D_->setAxisTitle("i#phi ", 2);
  hiDistrMB2Min2D_->setAxisTitle("i#eta ", 1);

  hiDistrNoise2Min2D_ = ibooker.book2D("Noise2depthMin1",
                                       "iphi-ieta noise distribution at depth1",
                                       hiDistr_x_nbin_,
                                       hiDistr_x_min_,
                                       hiDistr_x_max_,
                                       hiDistr_y_nbin_,
                                       hiDistr_y_min_,
                                       hiDistr_y_max_);

  hiDistrNoise2Min2D_->setAxisTitle("i#phi ", 2);
  hiDistrNoise2Min2D_->setAxisTitle("i#eta ", 1);

  // Variance
  hiDistrVarMBMin2D_ = ibooker.book2D("VarMBdepthMin1",
                                      "iphi- +ieta signal distribution at depth1",
                                      hiDistr_x_nbin_,
                                      hiDistr_x_min_,
                                      hiDistr_x_max_,
                                      hiDistr_y_nbin_,
                                      hiDistr_y_min_,
                                      hiDistr_y_max_);

  hiDistrVarMBMin2D_->setAxisTitle("i#phi ", 2);
  hiDistrVarMBMin2D_->setAxisTitle("i#eta ", 1);

  hiDistrVarNoiseMin2D_ = ibooker.book2D("VarNoisedepthMin1",
                                         "iphi-ieta noise distribution at depth1",
                                         hiDistr_x_nbin_,
                                         hiDistr_x_min_,
                                         hiDistr_x_max_,
                                         hiDistr_y_nbin_,
                                         hiDistr_y_min_,
                                         hiDistr_y_max_);

  hiDistrVarNoiseMin2D_->setAxisTitle("i#phi ", 2);
  hiDistrVarNoiseMin2D_->setAxisTitle("i#eta ", 1);
}

//--------------------------------------------------------
// void DQMHcalPhiSymAlCaReco::beginRun(const edm::Run& r, const EventSetup&
// context) {
//}

//--------------------------------------------------------

//-------------------------------------------------------------

void DQMHcalPhiSymAlCaReco::analyze(const Event &iEvent, const EventSetup &iSetup) {
  eventCounter_++;

  edm::Handle<FEDRawDataCollection> rawIn;
  iEvent.getByToken(rawInLabel_, rawIn);

  if (!rawIn.isValid()) {
    LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
    return;
  }

  // get HCAL FEDs:
  std::vector<int> selFEDs;
  for (int i = FEDNumbering::MINHCALFEDID; i <= FEDNumbering::MAXHCALFEDID; i++) {
    selFEDs.push_back(i);
  }

  //  std::cout<<" Size of FED "<<selFEDs.size()<<std::endl;

  const FEDRawDataCollection *rdc = rawIn.product();

  bool hcalIsZS = false;
  int lvl1ID = 0;
  bool lvl1IDFound = false;
  for (unsigned int k = 0; k < selFEDs.size(); k++) {
    const FEDRawData &fedData = rdc->FEDData(selFEDs[k]);
    // std::cout<<fedData.size()*std::pow(1024.,-1)<<std::endl;
    hFEDsize->Fill(fedData.size() * std::pow(1024., -1), 1);

    // get HCAL DCC Header for each FEDRawData
    const HcalDCCHeader *dccHeader = (const HcalDCCHeader *)(fedData.data());
    if (dccHeader) {
      // walk through the HTR data...
      HcalHTRData htr;

      int nspigot = 0;
      for (int spigot = 0; spigot < HcalDCCHeader::SPIGOT_COUNT; spigot++) {
        nspigot++;

        if (!dccHeader->getSpigotPresent(spigot))
          continue;

        // Load the given decoder with the pointer and length from this spigot.
        dccHeader->getSpigotData(spigot, htr, fedData.size());

        if (k != 20 && nspigot != 14) {
          if (!htr.isUnsuppressed()) {
            hcalIsZS = true;
          }
        }
      }
    }

    // try to get the lvl1ID from the HCAL fed
    if (!lvl1IDFound && (fedData.size() > 0)) {
      // get FED Header for FEDRawData
      FEDHeader fedHeader(fedData.data());
      lvl1ID = fedHeader.lvl1ID();
      lvl1IDFound = true;
    }
  }  // loop over HcalFEDs

  hHcalIsZS->Fill(hcalIsZS);
  hL1Id->Fill(lvl1ID % period_);

  edm::Handle<HBHERecHitCollection> hbheNS;
  iEvent.getByToken(hbherecoNoise, hbheNS);

  if (!hbheNS.isValid()) {
    LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
    return;
  }

  edm::Handle<HBHERecHitCollection> hbheMB;
  iEvent.getByToken(hbherecoMB, hbheMB);

  if (!hbheMB.isValid()) {
    LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
    return;
  }

  edm::Handle<HFRecHitCollection> hfNS;
  iEvent.getByToken(hfrecoNoise, hfNS);

  if (!hfNS.isValid()) {
    LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
    return;
  }

  edm::Handle<HFRecHitCollection> hfMB;
  iEvent.getByToken(hfrecoMB, hfMB);

  if (!hfMB.isValid()) {
    LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
    return;
  }

  const HBHERecHitCollection HithbheNS = *(hbheNS.product());

  hiDistrHBHEsize1D_->Fill(HithbheNS.size() / ihbhe_size_);

  for (HBHERecHitCollection::const_iterator hbheItr = HithbheNS.begin(); hbheItr != HithbheNS.end(); hbheItr++) {
    DetId id = (*hbheItr).detid();
    HcalDetId hid = HcalDetId(id);

    if (hid.depth() == 1) {
      if (hid.ieta() > 0) {
        hiDistrNoisePl2D_->Fill(hid.ieta(), hid.iphi(), hbheItr->energy());
        hiDistrNoise2Pl2D_->Fill(hid.ieta(), hid.iphi(), hbheItr->energy() * hbheItr->energy());
      } else {
        hiDistrNoiseMin2D_->Fill(fabs(hid.ieta()), hid.iphi(), hbheItr->energy());
        hiDistrNoise2Min2D_->Fill(fabs(hid.ieta()), hid.iphi(), hbheItr->energy() * hbheItr->energy());
      }
    }
  }

  const HBHERecHitCollection HithbheMB = *(hbheMB.product());

  for (HBHERecHitCollection::const_iterator hbheItr = HithbheMB.begin(); hbheItr != HithbheMB.end(); hbheItr++) {
    DetId id = (*hbheItr).detid();
    HcalDetId hid = HcalDetId(id);

    if (hid.depth() == 1) {
      if (hid.ieta() > 0) {
        hiDistrMBPl2D_->Fill(hid.ieta(), hid.iphi(), hbheItr->energy());
        hiDistrMB2Pl2D_->Fill(hid.ieta(), hid.iphi(), hbheItr->energy() * hbheItr->energy());
      } else {
        hiDistrMBMin2D_->Fill(fabs(hid.ieta()), hid.iphi(), hbheItr->energy());
        hiDistrMB2Min2D_->Fill(fabs(hid.ieta()), hid.iphi(), hbheItr->energy() * hbheItr->energy());
      }
    }
  }

  const HFRecHitCollection HithfNS = *(hfNS.product());

  hiDistrHFsize1D_->Fill(HithfNS.size() / ihf_size_);

  for (HFRecHitCollection::const_iterator hbheItr = HithfNS.begin(); hbheItr != HithfNS.end(); hbheItr++) {
    DetId id = (*hbheItr).detid();
    HcalDetId hid = HcalDetId(id);

    if (hid.depth() == 1) {
      if (hid.ieta() > 0) {
        hiDistrNoisePl2D_->Fill(hid.ieta(), hid.iphi(), hbheItr->energy());
        hiDistrNoise2Pl2D_->Fill(hid.ieta(), hid.iphi(), hbheItr->energy() * hbheItr->energy());
      } else {
        hiDistrNoiseMin2D_->Fill(fabs(hid.ieta()), hid.iphi(), hbheItr->energy());
        hiDistrNoise2Min2D_->Fill(fabs(hid.ieta()), hid.iphi(), hbheItr->energy() * hbheItr->energy());
      }
    }
  }

  const HFRecHitCollection HithfMB = *(hfMB.product());

  for (HFRecHitCollection::const_iterator hbheItr = HithfMB.begin(); hbheItr != HithfMB.end(); hbheItr++) {
    DetId id = (*hbheItr).detid();
    HcalDetId hid = HcalDetId(id);

    if (hid.depth() == 1) {
      if (hid.ieta() > 0) {
        hiDistrMBPl2D_->Fill(hid.ieta(), hid.iphi(), hbheItr->energy());
        hiDistrMB2Pl2D_->Fill(hid.ieta(), hid.iphi(), hbheItr->energy() * hbheItr->energy());
      } else {
        hiDistrMBMin2D_->Fill(fabs(hid.ieta()), hid.iphi(), hbheItr->energy());
        hiDistrMB2Min2D_->Fill(fabs(hid.ieta()), hid.iphi(), hbheItr->energy() * hbheItr->energy());
      }
    }
  }

}  // analyze

//--------------------------------------------------------
//--------------------------------------------------------
void DQMHcalPhiSymAlCaReco::endRun(const Run &r, const EventSetup &context) {
  // Keep Variances
  if (eventCounter_ > 0) {
    for (int k = 0; k <= hiDistr_x_nbin_; k++) {
      for (int j = 0; j <= hiDistr_y_nbin_; j++) {
        // First moment
        float cc1 = hiDistrMBPl2D_->getBinContent(k, j);
        cc1 = cc1 * 1. / eventCounter_;
        float cc2 = hiDistrNoisePl2D_->getBinContent(k, j);
        cc2 = cc2 * 1. / eventCounter_;
        float cc3 = hiDistrMBMin2D_->getBinContent(k, j);
        cc3 = cc3 * 1. / eventCounter_;
        float cc4 = hiDistrNoiseMin2D_->getBinContent(k, j);
        cc4 = cc4 * 1. / eventCounter_;
        // Second moment
        float cc11 = hiDistrMB2Pl2D_->getBinContent(k, j);
        cc11 = cc11 * 1. / eventCounter_;
        hiDistrVarMBPl2D_->setBinContent(k, j, cc11 - cc1 * cc1);
        float cc22 = hiDistrNoise2Pl2D_->getBinContent(k, j);
        cc22 = cc22 * 1. / eventCounter_;
        hiDistrVarNoisePl2D_->setBinContent(k, j, cc22 - cc2 * cc2);
        float cc33 = hiDistrMB2Min2D_->getBinContent(k, j);
        cc33 = cc33 * 1. / eventCounter_;
        hiDistrVarMBMin2D_->setBinContent(k, j, cc33 - cc3 * cc3);
        float cc44 = hiDistrNoise2Min2D_->getBinContent(k, j);
        cc44 = cc44 * 1. / eventCounter_;
        hiDistrVarNoiseMin2D_->setBinContent(k, j, cc44 - cc4 * cc4);
      }
    }
  }
}