DTMonitorModule/src/DTTriggerEfficiencyTask.cc

0001 /*
0002  * \file DTTriggerEfficiencyTask.cc
0003  *
0004  * \author C.Battilana - CIEMAT
0005  *
0006  */
0007
0008 #include "DQM/DTMonitorModule/src/DTTriggerEfficiencyTask.h"
0009
0010 // Framework
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012
0013 // DT trigger
0014 #include "DQM/DTMonitorModule/interface/DTTrigGeomUtils.h"
0015
0016 // Geometry
0017 #include "DataFormats/GeometryVector/interface/Pi.h"
0018 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0019 #include "Geometry/DTGeometry/interface/DTLayer.h"
0020 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
0021 #include "Geometry/DTGeometry/interface/DTTopology.h"
0022
0023 // DT Digi
0024 #include "DataFormats/DTDigi/interface/DTDigi.h"
0025 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
0026
0027 // Muon tracks
0028 #include "DataFormats/MuonReco/interface/Muon.h"
0029
0030 //Root
0031 #include "TH1.h"
0032 #include "TAxis.h"
0033
0034 #include <sstream>
0035 #include <iostream>
0036 #include <fstream>
0037
0038 using namespace edm;
0039 using namespace std;
0040
0041 DTTriggerEfficiencyTask::DTTriggerEfficiencyTask(const edm::ParameterSet& ps)
0042     : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()), trigGeomUtils(nullptr) {
0043   LogTrace("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask") << "[DTTriggerEfficiencyTask]: Constructor" << endl;
0044
0045   parameters = ps;
0046
0047   muons_Token_ = consumes<reco::MuonCollection>(parameters.getUntrackedParameter<edm::InputTag>("inputTagMuons"));
0048   tm_Token_ = consumes<L1MuDTChambPhContainer>(parameters.getUntrackedParameter<edm::InputTag>("inputTagTM"));
0049   inputTagSEG = parameters.getUntrackedParameter<edm::InputTag>("inputTagSEG");
0050   gmt_Token_ = consumes<L1MuGMTReadoutCollection>(parameters.getUntrackedParameter<edm::InputTag>("inputTagGMT"));
0051
0052   SegmArbitration = parameters.getUntrackedParameter<std::string>("SegmArbitration");
0053
0054   detailedPlots = parameters.getUntrackedParameter<bool>("detailedAnalysis");
0055   processTM = parameters.getUntrackedParameter<bool>("processTM");
0056
0057   checkRPCtriggers = parameters.getUntrackedParameter<bool>("checkRPCtriggers");
0058   nMinHitsPhi = parameters.getUntrackedParameter<int>("nMinHitsPhi");
0059   phiAccRange = parameters.getUntrackedParameter<double>("phiAccRange");
0060
0061   if (processTM)
0062     processTags.push_back("TM");
0063   if (!processTM)
0064     LogError("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask")
0065         << "[DTTriggerEfficiencyTask]: Error, no trigger source (Twinmux) has been selected!!" << endl;
0066 }
0067
0068 DTTriggerEfficiencyTask::~DTTriggerEfficiencyTask() {
0069   LogTrace("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask")
0070       << "[DTTriggerEfficiencyTask]: analyzed " << nevents << " events" << endl;
0071 }
0072
0073 void DTTriggerEfficiencyTask::dqmBeginRun(const edm::Run& run, const edm::EventSetup& context) {
0074   // Get the geometry
0075   muonGeom = &context.getData(muonGeomToken_);
0076   trigGeomUtils = new DTTrigGeomUtils(muonGeom);
0077 }
0078
0079 void DTTriggerEfficiencyTask::bookHistograms(DQMStore::IBooker& ibooker,
0080                                              edm::Run const& run,
0081                                              edm::EventSetup const& context) {
0082   LogTrace("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask") << "[DTTriggerEfficiencyTask]: bookHistograms" << endl;
0083
0084   nevents = 0;
0085   for (int wh = -2; wh <= 2; ++wh) {
0086     vector<string>::const_iterator tagIt = processTags.begin();
0087     vector<string>::const_iterator tagEnd = processTags.end();
0088     for (; tagIt != tagEnd; ++tagIt) {
0089       bookWheelHistos(ibooker, wh, (*tagIt), "Task");
0090       if (detailedPlots) {
0091         for (int stat = 1; stat <= 4; ++stat) {
0092           for (int sect = 1; sect <= 12; ++sect) {
0093             bookChamberHistos(ibooker, DTChamberId(wh, stat, sect), (*tagIt), "Segment");
0094           }
0095         }
0096       }
0097     }
0098   }
0099 }
0100
0101 void DTTriggerEfficiencyTask::analyze(const edm::Event& e, const edm::EventSetup& c) {
0102   nevents++;
0103
0104   if (checkRPCtriggers) {  //For pre-2016 Era compatibility
0105     if (!hasRPCTriggers(e)) {
0106       return;
0107     }
0108   }
0109
0110   map<DTChamberId, const L1MuDTChambPhDigi*> phBestTM;
0111   // Getting best TM Stuff
0112   edm::Handle<L1MuDTChambPhContainer> l1DTTPGPh;
0113   e.getByToken(tm_Token_, l1DTTPGPh);
0114   vector<L1MuDTChambPhDigi> const* phTrigs = l1DTTPGPh->getContainer();
0115   //empty from dttfDigis, needs emulator working?
0116   vector<L1MuDTChambPhDigi>::const_iterator iph = phTrigs->begin();
0117   vector<L1MuDTChambPhDigi>::const_iterator iphe = phTrigs->end();
0118   for (; iph != iphe; ++iph) {
0119     int phwheel = iph->whNum();
0120     int phsec = iph->scNum() + 1;  // DTTF numbering [0:11] -> DT numbering [1:12]
0121     int phst = iph->stNum();
0122     int phcode = iph->code();
0123
0124     DTChamberId chId(phwheel, phst, phsec);
0125
0126     if (phcode < 7 && (phBestTM.find(chId) == phBestTM.end() || phcode > phBestTM[chId]->code()))
0127       phBestTM[chId] = &(*iph);
0128   }
0129
0130   //Getting Best Segments
0131   vector<const DTRecSegment4D*> best4DSegments;
0132
0133   Handle<reco::MuonCollection> muons;
0134   e.getByToken(muons_Token_, muons);
0135   reco::MuonCollection::const_iterator mu;
0136
0137   for (mu = muons->begin(); mu != muons->end(); ++mu) {
0138     // Make sure that is standalone muon
0139     if (!((*mu).isStandAloneMuon())) {
0140       continue;
0141     }
0142
0143     // Get the chambers compatible with the muon
0144     const vector<reco::MuonChamberMatch> matchedChambers = (*mu).matches();
0145     vector<reco::MuonChamberMatch>::const_iterator chamber;
0146
0147     for (chamber = matchedChambers.begin(); chamber != matchedChambers.end(); ++chamber) {
0148       // look only in DTs
0149       if (chamber->detector() != MuonSubdetId::DT) {
0150         continue;
0151       }
0152
0153       // Get the matched segments in the chamber
0154       const vector<reco::MuonSegmentMatch> matchedSegments = (*chamber).segmentMatches;
0155       vector<reco::MuonSegmentMatch>::const_iterator segment;
0156
0157       for (segment = matchedSegments.begin(); segment != matchedSegments.end(); ++segment) {
0158         edm::Ref<DTRecSegment4DCollection> dtSegment = segment->dtSegmentRef;
0159
0160         // Segment Arbitration
0161         if (SegmArbitration == "SegmentArbitration" &&
0162             !((*segment).isMask(reco::MuonSegmentMatch::BestInChamberByDR))) {
0163           continue;
0164         }
0165
0166         if (SegmArbitration == "SegmentAndTrackArbitration" &&
0167             (!((*segment).isMask(reco::MuonSegmentMatch::BestInChamberByDR)) ||
0168              !((*segment).isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)))) {
0169           continue;
0170         }
0171
0172         if (SegmArbitration == "SegmentAndTrackArbitrationCleaned" &&
0173             (!((*segment).isMask(reco::MuonSegmentMatch::BestInChamberByDR)) ||
0174              !((*segment).isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) ||
0175              !((*segment).isMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning)))) {
0176           continue;
0177         }
0178
0179         if ((*dtSegment).hasPhi()) {
0180           best4DSegments.push_back(&(*dtSegment));
0181         }
0182
0183       }  // end loop on matched segments
0184     }  // end loop on compatible chambers
0185   }  // end loop on muons
0186
0187   // Plot filling
0188   vector<const DTRecSegment4D*>::const_iterator btrack;
0189   for (btrack = best4DSegments.begin(); btrack != best4DSegments.end(); ++btrack) {
0190     int wheel = (*btrack)->chamberId().wheel();
0191     int station = (*btrack)->chamberId().station();
0192     int scsector = 0;
0193     float x, xdir, y, ydir;
0194     trigGeomUtils->computeSCCoordinates((*btrack), scsector, x, xdir, y, ydir);
0195     int nHitsPhi = (*btrack)->phiSegment()->degreesOfFreedom() + 2;
0196     DTChamberId dtChId(wheel, station, scsector);
0197     uint32_t indexCh = dtChId.rawId();
0198     map<string, MonitorElement*>& innerChME = chamberHistos[indexCh];
0199     map<string, MonitorElement*>& innerWhME = wheelHistos[wheel];
0200
0201     if (fabs(xdir) < phiAccRange && nHitsPhi >= nMinHitsPhi) {
0202       vector<string>::const_iterator tagIt = processTags.begin();
0203       vector<string>::const_iterator tagEnd = processTags.end();
0204       for (; tagIt != tagEnd; ++tagIt) {
0205         int qual = phBestTM.find(dtChId) != phBestTM.end() ? phBestTM[dtChId]->code() : -1;
0206         innerWhME.find((*tagIt) + "_TrigEffDenum")->second->Fill(scsector, station);
0207         if (qual >= 0 && qual < 7) {
0208           innerWhME.find((*tagIt) + "_TrigEffNum")->second->Fill(scsector, station);
0209           if (qual >= 4) {
0210             innerWhME.find((*tagIt) + "_TrigEffCorrNum")->second->Fill(scsector, station);
0211           }
0212         }
0213         if (detailedPlots) {
0214           innerChME.find((*tagIt) + "_TrackPosvsAngle")->second->Fill(xdir, x);
0215           if (qual >= 0 && qual < 7) {
0216             innerChME.find((*tagIt) + "_TrackPosvsAngleAnyQual")->second->Fill(xdir, x);
0217             if (qual >= 4) {
0218               innerChME.find((*tagIt) + "_TrackPosvsAngleCorr")->second->Fill(xdir, x);
0219             }
0220           }
0221         }
0222       }
0223     }
0224   }
0225 }
0226
0227 bool DTTriggerEfficiencyTask::hasRPCTriggers(const edm::Event& e) {
0228   edm::Handle<L1MuGMTReadoutCollection> gmtrc;
0229   e.getByToken(gmt_Token_, gmtrc);
0230
0231   std::vector<L1MuGMTReadoutRecord> gmt_records = gmtrc->getRecords();
0232   std::vector<L1MuGMTReadoutRecord>::const_iterator igmtrr = gmt_records.begin();
0233   std::vector<L1MuGMTReadoutRecord>::const_iterator egmtrr = gmt_records.end();
0234   for (; igmtrr != egmtrr; igmtrr++) {
0235     std::vector<L1MuGMTExtendedCand> candsGMT = igmtrr->getGMTCands();
0236     std::vector<L1MuGMTExtendedCand>::const_iterator candGMTIt = candsGMT.begin();
0237     std::vector<L1MuGMTExtendedCand>::const_iterator candGMTEnd = candsGMT.end();
0238
0239     for (; candGMTIt != candGMTEnd; ++candGMTIt) {
0240       if (!candGMTIt->empty()) {
0241         int quality = candGMTIt->quality();
0242         if (candGMTIt->bx() == 0 && (quality == 5 || quality == 7)) {
0243           return true;
0244         }
0245       }
0246     }
0247   }
0248
0249   return false;
0250 }
0251
0252 void DTTriggerEfficiencyTask::bookChamberHistos(DQMStore::IBooker& ibooker,
0253                                                 const DTChamberId& dtCh,
0254                                                 string histoType,
0255                                                 string folder) {
0256   int wh = dtCh.wheel();
0257   int sc = dtCh.sector();
0258   int st = dtCh.station();
0259   stringstream wheel;
0260   wheel << wh;
0261   stringstream station;
0262   station << st;
0263   stringstream sector;
0264   sector << sc;
0265
0266   string hwFolder = topFolder();
0267   string bookingFolder =
0268       hwFolder + "Wheel" + wheel.str() + "/Sector" + sector.str() + "/Station" + station.str() + "/" + folder;
0269   string histoTag = "_W" + wheel.str() + "_Sec" + sector.str() + "_St" + station.str();
0270
0271   ibooker.setCurrentFolder(bookingFolder);
0272
0273   LogTrace("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask")
0274       << "[DTTriggerEfficiencyTask]: booking histos in " << bookingFolder << endl;
0275
0276   float min, max;
0277   int nbins;
0278   trigGeomUtils->phiRange(dtCh, min, max, nbins, 20);
0279
0280   string histoName = histoType + "_TrackPosvsAngle" + histoTag;
0281   string histoLabel = "Position vs Angle (phi)";
0282
0283   (chamberHistos[dtCh.rawId()])[histoType + "_TrackPosvsAngle"] =
0284       ibooker.book2D(histoName, histoLabel, 12, -30., 30., nbins, min, max);
0285
0286   histoName = histoType + "_TrackPosvsAngleAnyQual" + histoTag;
0287   histoLabel = "Position vs Angle (phi) for any qual triggers";
0288
0289   (chamberHistos[dtCh.rawId()])[histoType + "_TrackPosvsAngleAnyQual"] =
0290       ibooker.book2D(histoName, histoLabel, 12, -30., 30., nbins, min, max);
0291
0292   histoName = histoType + "_TrackPosvsAngleCorr" + histoTag;
0293   histoLabel = "Position vs Angle (phi) for correlated triggers";
0294
0295   (chamberHistos[dtCh.rawId()])[histoType + "_TrackPosvsAngleCorr"] =
0296       ibooker.book2D(histoName, histoLabel, 12, -30., 30., nbins, min, max);
0297 }
0298
0299 void DTTriggerEfficiencyTask::bookWheelHistos(DQMStore::IBooker& ibooker, int wheel, string hTag, string folder) {
0300   stringstream wh;
0301   wh << wheel;
0302   string basedir;
0303   if (hTag.find("Summary") != string::npos) {
0304     basedir = topFolder();  //Book summary histo outside folder directory
0305   } else {
0306     basedir = topFolder() + folder + "/";
0307   }
0308
0309   ibooker.setCurrentFolder(basedir);
0310
0311   string hTagName = "_W" + wh.str();
0312
0313   LogTrace("DTDQM|DTMonitorModule|DTTriggerEfficiencyTask")
0314       << "[DTTriggerEfficiencyTask]: booking histos in " << basedir << endl;
0315
0316   string hName = hTag + "_TrigEffDenum" + hTagName;
0317
0318   MonitorElement* me = ibooker.book2D(hName.c_str(), hName.c_str(), 12, 1, 13, 4, 1, 5);
0319
0320   me->setBinLabel(1, "MB1", 2);
0321   me->setBinLabel(2, "MB2", 2);
0322   me->setBinLabel(3, "MB3", 2);
0323   me->setBinLabel(4, "MB4", 2);
0324   me->setAxisTitle("Sector", 1);
0325
0326   wheelHistos[wheel][hTag + "_TrigEffDenum"] = me;
0327
0328   hName = hTag + "_TrigEffNum" + hTagName;
0329   me = ibooker.book2D(hName.c_str(), hName.c_str(), 12, 1, 13, 4, 1, 5);
0330
0331   me->setBinLabel(1, "MB1", 2);
0332   me->setBinLabel(2, "MB2", 2);
0333   me->setBinLabel(3, "MB3", 2);
0334   me->setBinLabel(4, "MB4", 2);
0335   me->setAxisTitle("Sector", 1);
0336
0337   wheelHistos[wheel][hTag + "_TrigEffNum"] = me;
0338
0339   hName = hTag + "_TrigEffCorrNum" + hTagName;
0340   me = ibooker.book2D(hName.c_str(), hName.c_str(), 12, 1, 13, 4, 1, 5);
0341
0342   me->setBinLabel(1, "MB1", 2);
0343   me->setBinLabel(2, "MB2", 2);
0344   me->setBinLabel(3, "MB3", 2);
0345   me->setBinLabel(4, "MB4", 2);
0346   me->setAxisTitle("Sector", 1);
0347
0348   wheelHistos[wheel][hTag + "_TrigEffCorrNum"] = me;
0349
0350   return;
0351 }