DTChamberEfficiencyClient.cc

Go to the documentation of this file.

/*
 *  See header file for a description of this class.
 *
 *  \author M. Pelliccioni - INFN Torino
 *
 *  threadsafe version (//-) oct/nov 2014 - WATWanAbdullah ncpp-um-my
 *
 */

#include "DQM/DTMonitorClient/src/DTChamberEfficiencyClient.h"
#include "DQMServices/Core/interface/DQMStore.h"

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

#include "Geometry/DTGeometry/interface/DTGeometry.h"

#include <cstdio>
#include <sstream>
#include <cmath>

using namespace edm;
using namespace std;

//two words about conventions: "All" histograms are those made for all segments
//while "Qual" histograms are those for segments with at least 12 hits

DTChamberEfficiencyClient::DTChamberEfficiencyClient(const ParameterSet& pSet)
    : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()) {
  LogVerbatim("DTDQM|DTMonitorClient|DTChamberEfficiencyClient") << "DTChamberEfficiencyClient: Constructor called";

  prescaleFactor = pSet.getUntrackedParameter<int>("diagnosticPrescale", 1);
}

DTChamberEfficiencyClient::~DTChamberEfficiencyClient() {
  LogVerbatim("DTDQM|DTMonitorClient|DTChamberEfficiencyClient") << "DTChamberEfficiencyClient: Destructor called";
}

void DTChamberEfficiencyClient::beginRun(const edm::Run& run, const edm::EventSetup& setup) {
  // Get the DT Geometry
  muonGeom = &setup.getData(muonGeomToken_);
}

void DTChamberEfficiencyClient::dqmEndLuminosityBlock(DQMStore::IBooker& ibooker,
                                                      DQMStore::IGetter& igetter,
                                                      edm::LuminosityBlock const& lumiSeg,
                                                      edm::EventSetup const& setup) {
  LogVerbatim("DTDQM|DTMonitorClient|DTChamberEfficiencyClient") << "DTChamberEfficiencyClient: endluminosityBlock";
}

void DTChamberEfficiencyClient::dqmEndJob(DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter) {
  LogVerbatim("DTDQM|DTMonitorClient|DTChamberEfficiencyClient") << "DTChamberEfficiencyClient: endRun";

  bookHistos(ibooker);

  // reset the global summary
  globalEffSummary->Reset();

  //Loop over the wheels
  for (int wheel = -2; wheel <= 2; wheel++) {
    stringstream wheel_str;
    wheel_str << wheel;

    // Get the ME produced by EfficiencyTask Source
    // All means no selection on segments, Qual means segments with at least 12 hits

    MonitorElement* MECountAll = igetter.get("DT/05-ChamberEff/Task/hCountSectVsChamb_All_W" + wheel_str.str());
    MonitorElement* MECountQual = igetter.get("DT/05-ChamberEff/Task/hCountSectVsChamb_Qual_W" + wheel_str.str());
    MonitorElement* MEExtrap = igetter.get("DT/05-ChamberEff/Task/hExtrapSectVsChamb_W" + wheel_str.str());

    //get the TH2F
    if (!MECountAll || !(MECountAll->getTH2F())) {
      edm::LogWarning("DTChamberEfficiencyClient") << "ME not available" << std::endl;
      return;
    }

    TH2F* hCountAll = MECountAll->getTH2F();
    TH2F* hCountQual = MECountQual->getTH2F();
    TH2F* hExtrap = MEExtrap->getTH2F();

    const int nBinX = summaryHistos[wheel + 2][0]->getNbinsX();
    const int nBinY = summaryHistos[wheel + 2][0]->getNbinsY();

    for (int j = 1; j <= nBinX; j++) {
      for (int k = 1; k <= nBinY; k++) {
        summaryHistos[wheel + 2][0]->setBinContent(j, k, 0.);
        summaryHistos[wheel + 2][1]->setBinContent(j, k, 0.);

        const float numerAll = hCountAll->GetBinContent(j, k);
        const float numerQual = hCountQual->GetBinContent(j, k);
        const float denom = hExtrap->GetBinContent(j, k);

        if (denom != 0.) {
          const float effAll = numerAll / denom;
          const float eff_error_All = sqrt((effAll + effAll * effAll) / denom);

          const float effQual = numerQual / denom;
          const float eff_error_Qual = sqrt((effQual + effQual * effQual) / denom);

          summaryHistos[wheel + 2][0]->setBinContent(j, k, effAll);
          summaryHistos[wheel + 2][0]->setBinError(j, k, eff_error_All);

          summaryHistos[wheel + 2][1]->setBinContent(j, k, effQual);
          summaryHistos[wheel + 2][1]->setBinError(j, k, eff_error_Qual);

          // Fill 1D eff distributions
          globalEffDistr->Fill(effAll);
          EffDistrPerWh[wheel + 2]->Fill(effAll);
        }
      }
    }
  }

  // fill the global eff. summary
  // problems at a granularity smaller than the chamber are ignored
  for (int wheel = -2; wheel <= 2; wheel++) {  // loop over wheels
    // retrieve the chamber efficiency summary
    MonitorElement* segmentWheelSummary = summaryHistos[wheel + 2][0];
    if (segmentWheelSummary != nullptr) {
      for (int sector = 1; sector <= 12; sector++) {  // loop over sectors
        float nFailingChambers = 0.;

        double meaneff = 0.;
        double errorsum = 0.;

        for (int station = 1; station != 5; ++station) {  // loop over stations

          const double tmpefficiency = segmentWheelSummary->getBinContent(sector, station);
          const double tmpvariance = pow(segmentWheelSummary->getBinError(sector, station), 2);

          if (tmpefficiency < 0.2 || tmpvariance == 0) {
            nFailingChambers++;
            continue;
          }

          meaneff += tmpefficiency / tmpvariance;
          errorsum += 1. / tmpvariance;

          LogTrace("DTDQM|DTMonitorClient|DTChamberEfficiencyClient")
              << "Wheel: " << wheel << " Stat: " << station << " Sect: " << sector << " status: " << meaneff / errorsum
              << endl;
        }

        if (sector == 4 || sector == 10) {
          int whichSector = (sector == 4) ? 13 : 14;

          const double tmpefficiency = segmentWheelSummary->getBinContent(whichSector, 4);
          const double tmpvariance = pow(segmentWheelSummary->getBinError(whichSector, 4), 2);

          if (tmpefficiency > 0.2 && tmpvariance != 0) {
            meaneff += tmpefficiency / tmpvariance;
            errorsum += 1. / tmpvariance;
          } else
            nFailingChambers++;
        }

        double eff_result = 0;
        if (errorsum != 0)
          eff_result = meaneff / errorsum;

        if (nFailingChambers != 0) {
          if (sector != 4 && sector != 10)
            eff_result = eff_result * (4. - nFailingChambers) / 4.;
          else
            eff_result = eff_result * (5. - nFailingChambers) / 5.;
        }

        if (eff_result > 0.7)
          globalEffSummary->Fill(sector, wheel, 1.);
        else if (eff_result < 0.7 && eff_result > 0.5)
          globalEffSummary->Fill(sector, wheel, 0.6);
        else if (eff_result < 0.5 && eff_result > 0.3)
          globalEffSummary->Fill(sector, wheel, 0.4);
        else if (eff_result < 0.3 && eff_result > 0.)
          globalEffSummary->Fill(sector, wheel, 0.15);
        else
          globalEffSummary->Fill(sector, wheel, 0.01);
      }
    }
  }
  return;
}

void DTChamberEfficiencyClient::bookHistos(DQMStore::IBooker& ibooker) {
  ibooker.setCurrentFolder("DT/05-ChamberEff");

  globalEffSummary = ibooker.book2D("EfficiencyGlbSummary", "Efficiency Summary", 12, 1, 13, 5, -2, 3);
  globalEffSummary->setAxisTitle("sector", 1);
  globalEffSummary->setAxisTitle("wheel", 2);

  globalEffDistr = ibooker.book1D("TotalEfficiency", "Total efficiency", 51, 0., 1.02);
  globalEffDistr->setAxisTitle("Eff", 1);

  for (int wh = -2; wh <= 2; wh++) {
    stringstream wheel;
    wheel << wh;
    string histoNameAll = "EfficiencyMap_All_W" + wheel.str();
    string histoTitleAll = "Efficiency map for all segments for wheel " + wheel.str();

    string histoNameQual = "EfficiencyMap_Qual_W" + wheel.str();
    string histoTitleQual = "Efficiency map for quality segments for wheel " + wheel.str();

    string histoNameEff = "Efficiency_W" + wheel.str();
    string histoTitleEff = "Segment efficiency, wheel " + wheel.str();

    ibooker.setCurrentFolder("DT/05-ChamberEff");

    summaryHistos[wh + 2][0] = ibooker.book2D(histoNameAll.c_str(), histoTitleAll.c_str(), 14, 1., 15., 4, 1., 5.);
    summaryHistos[wh + 2][0]->setAxisTitle("Sector", 1);
    summaryHistos[wh + 2][0]->setBinLabel(1, "MB1", 2);
    summaryHistos[wh + 2][0]->setBinLabel(2, "MB2", 2);
    summaryHistos[wh + 2][0]->setBinLabel(3, "MB3", 2);
    summaryHistos[wh + 2][0]->setBinLabel(4, "MB4", 2);

    EffDistrPerWh[wh + 2] = ibooker.book1D(histoNameEff.c_str(), histoTitleEff.c_str(), 51, 0., 1.02);
    EffDistrPerWh[wh + 2]->setAxisTitle("Eff", 1);

    ibooker.setCurrentFolder("DT/05-ChamberEff/HighQual");

    summaryHistos[wh + 2][1] = ibooker.book2D(histoNameQual.c_str(), histoTitleQual.c_str(), 14, 1., 15., 4, 1., 5.);
    summaryHistos[wh + 2][1]->setAxisTitle("Sector", 1);
    summaryHistos[wh + 2][1]->setBinLabel(1, "MB1", 2);
    summaryHistos[wh + 2][1]->setBinLabel(2, "MB2", 2);
    summaryHistos[wh + 2][1]->setBinLabel(3, "MB3", 2);
    summaryHistos[wh + 2][1]->setBinLabel(4, "MB4", 2);
  }

  return;
}