DTTriggerEfficiencyTask.cc

Go to the documentation of this file.

/*
 * \file DTTriggerEfficiencyTask.cc
 *
 * \author C.Battilana - CIEMAT
 *
 */

#include "DQM/DTMonitorModule/src/DTTriggerEfficiencyTask.h"

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

// DT trigger
#include "DQM/DTMonitorModule/interface/DTTrigGeomUtils.h"

// Geometry
#include "DataFormats/GeometryVector/interface/Pi.h"
#include "Geometry/Records/interface/MuonGeometryRecord.h"
#include "Geometry/DTGeometry/interface/DTGeometry.h"
#include "Geometry/DTGeometry/interface/DTLayer.h"
#include "Geometry/DTGeometry/interface/DTSuperLayer.h"
#include "Geometry/DTGeometry/interface/DTTopology.h"

// DT Digi
#include <DataFormats/DTDigi/interface/DTDigi.h>
#include <DataFormats/DTDigi/interface/DTDigiCollection.h>

// Muon tracks
#include <DataFormats/MuonReco/interface/Muon.h>

//Root
#include "TH1.h"
#include "TAxis.h"

#include <sstream>
#include <iostream>
#include <fstream>

using namespace edm;
using namespace std;

DTTriggerEfficiencyTask::DTTriggerEfficiencyTask(const edm::ParameterSet& ps) : trigGeomUtils(nullptr) {
  LogTrace("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask") << "[DTTriggerEfficiencyTask]: Constructor" << endl;

  parameters = ps;

  muons_Token_ = consumes<reco::MuonCollection>(parameters.getUntrackedParameter<edm::InputTag>("inputTagMuons"));
  tm_Token_ = consumes<L1MuDTChambPhContainer>(parameters.getUntrackedParameter<edm::InputTag>("inputTagTM"));
  inputTagSEG = parameters.getUntrackedParameter<edm::InputTag>("inputTagSEG");
  gmt_Token_ = consumes<L1MuGMTReadoutCollection>(parameters.getUntrackedParameter<edm::InputTag>("inputTagGMT"));

  SegmArbitration = parameters.getUntrackedParameter<std::string>("SegmArbitration");

  detailedPlots = parameters.getUntrackedParameter<bool>("detailedAnalysis");
  processTM = parameters.getUntrackedParameter<bool>("processTM");

  checkRPCtriggers = parameters.getUntrackedParameter<bool>("checkRPCtriggers");
  nMinHitsPhi = parameters.getUntrackedParameter<int>("nMinHitsPhi");
  phiAccRange = parameters.getUntrackedParameter<double>("phiAccRange");

  if (processTM)
    processTags.push_back("TM");
  if (!processTM)
    LogError("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask")
        << "[DTTriggerEfficiencyTask]: Error, no trigger source (Twinmux) has been selected!!" << endl;
}

DTTriggerEfficiencyTask::~DTTriggerEfficiencyTask() {
  LogTrace("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask")
      << "[DTTriggerEfficiencyTask]: analyzed " << nevents << " events" << endl;
}

void DTTriggerEfficiencyTask::dqmBeginRun(const edm::Run& run, const edm::EventSetup& context) {
  // Get the geometry
  context.get<MuonGeometryRecord>().get(muonGeom);
  trigGeomUtils = new DTTrigGeomUtils(muonGeom);
}

void DTTriggerEfficiencyTask::bookHistograms(DQMStore::IBooker& ibooker,
                                             edm::Run const& run,
                                             edm::EventSetup const& context) {
  LogTrace("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask") << "[DTTriggerEfficiencyTask]: bookHistograms" << endl;

  nevents = 0;
  for (int wh = -2; wh <= 2; ++wh) {
    vector<string>::const_iterator tagIt = processTags.begin();
    vector<string>::const_iterator tagEnd = processTags.end();
    for (; tagIt != tagEnd; ++tagIt) {
      bookWheelHistos(ibooker, wh, (*tagIt), "Task");
      if (detailedPlots) {
        for (int stat = 1; stat <= 4; ++stat) {
          for (int sect = 1; sect <= 12; ++sect) {
            bookChamberHistos(ibooker, DTChamberId(wh, stat, sect), (*tagIt), "Segment");
          }
        }
      }
    }
  }
}

void DTTriggerEfficiencyTask::analyze(const edm::Event& e, const edm::EventSetup& c) {
  nevents++;

  if (checkRPCtriggers) {  //For pre-2016 Era compatibility
    if (!hasRPCTriggers(e)) {
      return;
    }
  }

  map<DTChamberId, const L1MuDTChambPhDigi*> phBestTM;
  // Getting best TM Stuff
  edm::Handle<L1MuDTChambPhContainer> l1DTTPGPh;
  e.getByToken(tm_Token_, l1DTTPGPh);
  vector<L1MuDTChambPhDigi> const* phTrigs = l1DTTPGPh->getContainer();
  //empty from dttfDigis, needs emulator working?
  vector<L1MuDTChambPhDigi>::const_iterator iph = phTrigs->begin();
  vector<L1MuDTChambPhDigi>::const_iterator iphe = phTrigs->end();
  for (; iph != iphe; ++iph) {
    int phwheel = iph->whNum();
    int phsec = iph->scNum() + 1;  // DTTF numbering [0:11] -> DT numbering [1:12]
    int phst = iph->stNum();
    int phcode = iph->code();

    DTChamberId chId(phwheel, phst, phsec);

    if (phcode < 7 && (phBestTM.find(chId) == phBestTM.end() || phcode > phBestTM[chId]->code()))
      phBestTM[chId] = &(*iph);
  }

  //Getting Best Segments
  vector<const DTRecSegment4D*> best4DSegments;

  Handle<reco::MuonCollection> muons;
  e.getByToken(muons_Token_, muons);
  reco::MuonCollection::const_iterator mu;

  for (mu = muons->begin(); mu != muons->end(); ++mu) {
    // Make sure that is standalone muon
    if (!((*mu).isStandAloneMuon())) {
      continue;
    }

    // Get the chambers compatible with the muon
    const vector<reco::MuonChamberMatch> matchedChambers = (*mu).matches();
    vector<reco::MuonChamberMatch>::const_iterator chamber;

    for (chamber = matchedChambers.begin(); chamber != matchedChambers.end(); ++chamber) {
      // look only in DTs
      if (chamber->detector() != MuonSubdetId::DT) {
        continue;
      }

      // Get the matched segments in the chamber
      const vector<reco::MuonSegmentMatch> matchedSegments = (*chamber).segmentMatches;
      vector<reco::MuonSegmentMatch>::const_iterator segment;

      for (segment = matchedSegments.begin(); segment != matchedSegments.end(); ++segment) {
        edm::Ref<DTRecSegment4DCollection> dtSegment = segment->dtSegmentRef;

        // Segment Arbitration
        if (SegmArbitration == "SegmentArbitration" &&
            !((*segment).isMask(reco::MuonSegmentMatch::BestInChamberByDR))) {
          continue;
        }

        if (SegmArbitration == "SegmentAndTrackArbitration" &&
            (!((*segment).isMask(reco::MuonSegmentMatch::BestInChamberByDR)) ||
             !((*segment).isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)))) {
          continue;
        }

        if (SegmArbitration == "SegmentAndTrackArbitrationCleaned" &&
            (!((*segment).isMask(reco::MuonSegmentMatch::BestInChamberByDR)) ||
             !((*segment).isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) ||
             !((*segment).isMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning)))) {
          continue;
        }

        if ((*dtSegment).hasPhi()) {
          best4DSegments.push_back(&(*dtSegment));
        }

      }  // end loop on matched segments
    }    // end loop on compatible chambers
  }      // end loop on muons

  // Plot filling
  vector<const DTRecSegment4D*>::const_iterator btrack;
  for (btrack = best4DSegments.begin(); btrack != best4DSegments.end(); ++btrack) {
    int wheel = (*btrack)->chamberId().wheel();
    int station = (*btrack)->chamberId().station();
    int scsector = 0;
    float x, xdir, y, ydir;
    trigGeomUtils->computeSCCoordinates((*btrack), scsector, x, xdir, y, ydir);
    int nHitsPhi = (*btrack)->phiSegment()->degreesOfFreedom() + 2;
    DTChamberId dtChId(wheel, station, scsector);
    uint32_t indexCh = dtChId.rawId();
    map<string, MonitorElement*>& innerChME = chamberHistos[indexCh];
    map<string, MonitorElement*>& innerWhME = wheelHistos[wheel];

    if (fabs(xdir) < phiAccRange && nHitsPhi >= nMinHitsPhi) {
      vector<string>::const_iterator tagIt = processTags.begin();
      vector<string>::const_iterator tagEnd = processTags.end();
      for (; tagIt != tagEnd; ++tagIt) {
        int qual = phBestTM.find(dtChId) != phBestTM.end() ? phBestTM[dtChId]->code() : -1;
        innerWhME.find((*tagIt) + "_TrigEffDenum")->second->Fill(scsector, station);
        if (qual >= 0 && qual < 7) {
          innerWhME.find((*tagIt) + "_TrigEffNum")->second->Fill(scsector, station);
          if (qual >= 4) {
            innerWhME.find((*tagIt) + "_TrigEffCorrNum")->second->Fill(scsector, station);
          }
        }
        if (detailedPlots) {
          innerChME.find((*tagIt) + "_TrackPosvsAngle")->second->Fill(xdir, x);
          if (qual >= 0 && qual < 7) {
            innerChME.find((*tagIt) + "_TrackPosvsAngleAnyQual")->second->Fill(xdir, x);
            if (qual >= 4) {
              innerChME.find((*tagIt) + "_TrackPosvsAngleCorr")->second->Fill(xdir, x);
            }
          }
        }
      }
    }
  }
}

bool DTTriggerEfficiencyTask::hasRPCTriggers(const edm::Event& e) {
  edm::Handle<L1MuGMTReadoutCollection> gmtrc;
  e.getByToken(gmt_Token_, gmtrc);

  std::vector<L1MuGMTReadoutRecord> gmt_records = gmtrc->getRecords();
  std::vector<L1MuGMTReadoutRecord>::const_iterator igmtrr = gmt_records.begin();
  std::vector<L1MuGMTReadoutRecord>::const_iterator egmtrr = gmt_records.end();
  for (; igmtrr != egmtrr; igmtrr++) {
    std::vector<L1MuGMTExtendedCand> candsGMT = igmtrr->getGMTCands();
    std::vector<L1MuGMTExtendedCand>::const_iterator candGMTIt = candsGMT.begin();
    std::vector<L1MuGMTExtendedCand>::const_iterator candGMTEnd = candsGMT.end();

    for (; candGMTIt != candGMTEnd; ++candGMTIt) {
      if (!candGMTIt->empty()) {
        int quality = candGMTIt->quality();
        if (candGMTIt->bx() == 0 && (quality == 5 || quality == 7)) {
          return true;
        }
      }
    }
  }

  return false;
}

void DTTriggerEfficiencyTask::bookChamberHistos(DQMStore::IBooker& ibooker,
                                                const DTChamberId& dtCh,
                                                string histoType,
                                                string folder) {
  int wh = dtCh.wheel();
  int sc = dtCh.sector();
  int st = dtCh.station();
  stringstream wheel;
  wheel << wh;
  stringstream station;
  station << st;
  stringstream sector;
  sector << sc;

  string hwFolder = topFolder();
  string bookingFolder =
      hwFolder + "Wheel" + wheel.str() + "/Sector" + sector.str() + "/Station" + station.str() + "/" + folder;
  string histoTag = "_W" + wheel.str() + "_Sec" + sector.str() + "_St" + station.str();

  ibooker.setCurrentFolder(bookingFolder);

  LogTrace("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask")
      << "[DTTriggerEfficiencyTask]: booking histos in " << bookingFolder << endl;

  float min, max;
  int nbins;
  trigGeomUtils->phiRange(dtCh, min, max, nbins, 20);

  string histoName = histoType + "_TrackPosvsAngle" + histoTag;
  string histoLabel = "Position vs Angle (phi)";

  (chamberHistos[dtCh.rawId()])[histoType + "_TrackPosvsAngle"] =
      ibooker.book2D(histoName, histoLabel, 12, -30., 30., nbins, min, max);

  histoName = histoType + "_TrackPosvsAngleAnyQual" + histoTag;
  histoLabel = "Position vs Angle (phi) for any qual triggers";

  (chamberHistos[dtCh.rawId()])[histoType + "_TrackPosvsAngleAnyQual"] =
      ibooker.book2D(histoName, histoLabel, 12, -30., 30., nbins, min, max);

  histoName = histoType + "_TrackPosvsAngleCorr" + histoTag;
  histoLabel = "Position vs Angle (phi) for correlated triggers";

  (chamberHistos[dtCh.rawId()])[histoType + "_TrackPosvsAngleCorr"] =
      ibooker.book2D(histoName, histoLabel, 12, -30., 30., nbins, min, max);
}

void DTTriggerEfficiencyTask::bookWheelHistos(DQMStore::IBooker& ibooker, int wheel, string hTag, string folder) {
  stringstream wh;
  wh << wheel;
  string basedir;
  if (hTag.find("Summary") != string::npos) {
    basedir = topFolder();  //Book summary histo outside folder directory
  } else {
    basedir = topFolder() + folder + "/";
  }

  ibooker.setCurrentFolder(basedir);

  string hTagName = "_W" + wh.str();

  LogTrace("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask")
      << "[DTTriggerEfficiencyTask]: booking histos in " << basedir << endl;

  string hName = hTag + "_TrigEffDenum" + hTagName;

  MonitorElement* me = ibooker.book2D(hName.c_str(), hName.c_str(), 12, 1, 13, 4, 1, 5);

  me->setBinLabel(1, "MB1", 2);
  me->setBinLabel(2, "MB2", 2);
  me->setBinLabel(3, "MB3", 2);
  me->setBinLabel(4, "MB4", 2);
  me->setAxisTitle("Sector", 1);

  wheelHistos[wheel][hTag + "_TrigEffDenum"] = me;

  hName = hTag + "_TrigEffNum" + hTagName;
  me = ibooker.book2D(hName.c_str(), hName.c_str(), 12, 1, 13, 4, 1, 5);

  me->setBinLabel(1, "MB1", 2);
  me->setBinLabel(2, "MB2", 2);
  me->setBinLabel(3, "MB3", 2);
  me->setBinLabel(4, "MB4", 2);
  me->setAxisTitle("Sector", 1);

  wheelHistos[wheel][hTag + "_TrigEffNum"] = me;

  hName = hTag + "_TrigEffCorrNum" + hTagName;
  me = ibooker.book2D(hName.c_str(), hName.c_str(), 12, 1, 13, 4, 1, 5);

  me->setBinLabel(1, "MB1", 2);
  me->setBinLabel(2, "MB2", 2);
  me->setBinLabel(3, "MB3", 2);
  me->setBinLabel(4, "MB4", 2);
  me->setAxisTitle("Sector", 1);

  wheelHistos[wheel][hTag + "_TrigEffCorrNum"] = me;

  return;
}