DTMonitorModule/src/DTChamberEfficiencyTask.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author G. Mila - INFN Torino
0005  */
0006
0007 #include "DTChamberEfficiencyTask.h"
0008
0009 //Framework
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013
0014 #include "DQMServices/Core/interface/DQMStore.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include <iterator>
0018 #include <iostream>
0019 #include <cmath>
0020 using namespace edm;
0021 using namespace std;
0022
0023 DTChamberEfficiencyTask::DTChamberEfficiencyTask(const ParameterSet& pset)
0024     : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()) {
0025   debug = pset.getUntrackedParameter<bool>("debug", false);
0026
0027   edm::LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiencyTask") << "[DTChamberEfficiencyTask] Constructor called!";
0028
0029   parameters = pset;
0030
0031   // the name of the 4D rec hits collection
0032   recHits4DToken_ =
0033       consumes<DTRecSegment4DCollection>(edm::InputTag(parameters.getUntrackedParameter<string>("recHits4DLabel")));
0034
0035   // parameters to use for the segment quality check
0036   theMinHitsSegment = static_cast<unsigned int>(parameters.getParameter<int>("minHitsSegment"));
0037   theMinChi2NormSegment = parameters.getParameter<double>("minChi2NormSegment");
0038   // parameter to use for the exstrapolated segment check
0039   theMinCloseDist = parameters.getParameter<double>("minCloseDist");
0040
0041   // the running modality
0042   onlineMonitor = parameters.getUntrackedParameter<bool>("onlineMonitor");
0043
0044   // the analysis mode
0045   detailedAnalysis = parameters.getUntrackedParameter<bool>("detailedAnalysis");
0046 }
0047
0048 DTChamberEfficiencyTask::~DTChamberEfficiencyTask() {
0049   edm::LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiencyTask") << "[DTChamberEfficiencyTask] Destructor called!";
0050 }
0051
0052 void DTChamberEfficiencyTask::beginLuminosityBlock(LuminosityBlock const& lumiSeg, EventSetup const& context) {
0053   edm::LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiencyTask")
0054       << "[DTChamberEfficiencyTask]: Begin of LS transition";
0055
0056   if (lumiSeg.id().luminosityBlock() % parameters.getUntrackedParameter<int>("ResetCycle", 3) == 0 && onlineMonitor) {
0057     for (map<DTChamberId, vector<MonitorElement*> >::const_iterator histo = histosPerCh.begin();
0058          histo != histosPerCh.end();
0059          histo++) {
0060       int size = (*histo).second.size();
0061       for (int i = 0; i < size; i++) {
0062         (*histo).second[i]->Reset();
0063       }
0064     }
0065   }
0066 }
0067
0068 void DTChamberEfficiencyTask::dqmBeginRun(const edm::Run& run, const edm::EventSetup& setup) {
0069   // Get the DT Geometry
0070   dtGeom = &setup.getData(muonGeomToken_);
0071 }
0072
0073 void DTChamberEfficiencyTask::bookHistograms(DQMStore::IBooker& ibooker,
0074                                              edm::Run const& iRun,
0075                                              edm::EventSetup const& context) {
0076   ibooker.setCurrentFolder("DT/DTChamberEfficiencyTask");
0077
0078   // Loop over all the chambers
0079   vector<const DTChamber*>::const_iterator ch_it = dtGeom->chambers().begin();
0080   vector<const DTChamber*>::const_iterator ch_end = dtGeom->chambers().end();
0081   for (; ch_it != ch_end; ++ch_it) {
0082     // histo booking
0083     bookHistos(ibooker, (*ch_it)->id());
0084   }
0085 }
0086
0087 // Book a set of histograms for a given Layer
0088 void DTChamberEfficiencyTask::bookHistos(DQMStore::IBooker& ibooker, DTChamberId chId) {
0089   edm::LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiencyTask") << "   Booking histos for CH : " << chId;
0090
0091   // Compose the chamber name
0092   stringstream wheel;
0093   wheel << chId.wheel();
0094   stringstream station;
0095   station << chId.station();
0096   stringstream sector;
0097   sector << chId.sector();
0098
0099   string HistoName = "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str();
0100
0101   ibooker.setCurrentFolder("DT/01-DTChamberEfficiency/Task/Wheel" + wheel.str() + "/Sector" + sector.str() +
0102                            "/Station" + station.str());
0103
0104   // Create the monitor elements
0105   vector<MonitorElement*> histos;
0106
0107   //efficiency selection cuts
0108   // a- number of segments of the top chamber > 0 && number of segments of the bottom chamber > 0
0109   // b- number of segments of the middle chamber > 0
0110   // c- check of the top and bottom segment quality
0111   // d- check if interpolation falls inside the middle chamber
0112   // e- check of the middle segment quality
0113   // f- check if the distance between the reconstructed and the exstrapolated segments is ok
0114
0115   // histo for efficiency with cuts a-/c-/d-
0116   histos.push_back(ibooker.book2D(
0117       "hEffGoodSegVsPosDen" + HistoName, "Eff vs local position (good) ", 25, -250., 250., 25, -250., 250.));
0118   // histo for efficiency with cuts a-/b-/c-/d-/e-/f-
0119   histos.push_back(ibooker.book2D("hEffGoodCloseSegVsPosNum" + HistoName,
0120                                   "Eff vs local position (good and close segs) ",
0121                                   25,
0122                                   -250.,
0123                                   250.,
0124                                   25,
0125                                   -250.,
0126                                   250.));
0127   if (detailedAnalysis) {
0128     histos.push_back(
0129         ibooker.book1D("hDistSegFromExtrap" + HistoName, "Distance segments from extrap position ", 200, 0., 200.));
0130     // histo for efficiency from segment counting
0131     histos.push_back(ibooker.book1D("hNaiveEffSeg" + HistoName, "Naive eff ", 10, 0., 10.));
0132     // histo for efficiency with cuts a-/c-
0133     histos.push_back(ibooker.book2D(
0134         "hEffSegVsPosDen" + HistoName, "Eff vs local position (all) ", 25, -250., 250., 25, -250., 250.));
0135     // histo for efficiency with cuts a-/b-/c-/d-
0136     histos.push_back(
0137         ibooker.book2D("hEffSegVsPosNum" + HistoName, "Eff vs local position ", 25, -250., 250., 25, -250., 250.));
0138     // histo for efficiency with cuts a-/b-/c-/d-/e-
0139     histos.push_back(ibooker.book2D(
0140         "hEffGoodSegVsPosNum" + HistoName, "Eff vs local position (good segs) ", 25, -250., 250., 25, -250., 250.));
0141   }
0142   histosPerCh[chId] = histos;
0143 }
0144
0145 void DTChamberEfficiencyTask::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0146   edm::LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiencyTask")
0147       << "[DTChamberEfficiencyTask] Analyze #Run: " << event.id().run() << " #Event: " << event.id().event();
0148
0149   // Get the 4D rechit collection from the event
0150   event.getByToken(recHits4DToken_, segs);
0151
0152   int bottom = 0, top = 0;
0153
0154   // Loop over all the chambers
0155   vector<const DTChamber*>::const_iterator ch_it = dtGeom->chambers().begin();
0156   vector<const DTChamber*>::const_iterator ch_end = dtGeom->chambers().end();
0157   for (; ch_it != ch_end; ++ch_it) {
0158     DTChamberId ch = (*ch_it)->id();
0159     int wheel = ch.wheel();
0160     int sector = ch.sector();
0161     int station = ch.station();
0162
0163     DTChamberId MidId(wheel, station, sector);
0164
0165     // get efficiency for MB1 using MB2 and MB3
0166     if (station == 1) {
0167       bottom = 2;
0168       top = 3;
0169     }
0170
0171     // get efficiency for MB2 using MB1 and MB3
0172     if (station == 2) {
0173       bottom = 1;
0174       top = 3;
0175     }
0176
0177     // get efficiency for MB2 using MB2 and MB4
0178     if (station == 3) {
0179       bottom = 2;
0180       top = 4;
0181     }
0182
0183     // get efficiency for MB4 using MB2 and MB3
0184     if (station == 4) {
0185       bottom = 2;
0186       top = 3;
0187     }
0188
0189     // Select events with (good) segments in Bot and Top
0190     DTChamberId BotId(wheel, bottom, sector);
0191     DTChamberId TopId(wheel, top, sector);
0192
0193     // Get segments in the bottom chambers (if any)
0194     DTRecSegment4DCollection::range segsBot = segs->get(BotId);
0195     int nSegsBot = segsBot.second - segsBot.first;
0196     // check if any segments is there
0197     if (nSegsBot == 0)
0198       continue;
0199
0200     vector<MonitorElement*> histos = histosPerCh[MidId];
0201
0202     // Get segments in the top chambers (if any)
0203     DTRecSegment4DCollection::range segsTop = segs->get(TopId);
0204     int nSegsTop = segsTop.second - segsTop.first;
0205
0206     // Select one segment for the bottom chamber
0207     const DTRecSegment4D& bestBotSeg = getBestSegment(segsBot);
0208
0209     // Select one segment for the top chamber
0210     DTRecSegment4D* pBestTopSeg = nullptr;
0211     if (nSegsTop > 0)
0212       pBestTopSeg = const_cast<DTRecSegment4D*>(&getBestSegment(segsTop));
0213     //if top chamber is MB4 sector 10, consider also sector 14
0214     if (TopId.station() == 4 && TopId.sector() == 10) {
0215       DTChamberId TopId14(wheel, top, 14);
0216       DTRecSegment4DCollection::range segsTop14 = segs->get(TopId14);
0217       int nSegsTop14 = segsTop14.second - segsTop14.first;
0218       nSegsTop += nSegsTop14;
0219       if (nSegsTop14) {
0220         DTRecSegment4D* pBestTopSeg14 = const_cast<DTRecSegment4D*>(&getBestSegment(segsTop14));
0221         // get best between sector 10 and 14
0222         pBestTopSeg = const_cast<DTRecSegment4D*>(getBestSegment(pBestTopSeg, pBestTopSeg14));
0223       }
0224     }
0225     if (!pBestTopSeg)
0226       continue;
0227     const DTRecSegment4D& bestTopSeg = *pBestTopSeg;
0228
0229     DTRecSegment4DCollection::range segsMid = segs->get(MidId);
0230     int nSegsMid = segsMid.second - segsMid.first;
0231
0232     if (detailedAnalysis) {
0233       // very trivial efficiency, just count segments
0234       histos[3]->Fill(0);
0235       if (nSegsMid > 0)
0236         histos[3]->Fill(1);
0237     }
0238
0239     // get position at Mid by interpolating the position (not direction) of best
0240     // segment in Bot and Top to Mid surface
0241     LocalPoint posAtMid = interpolate(bestBotSeg, bestTopSeg, MidId);
0242
0243     // is best segment good enough?
0244     if (isGoodSegment(bestBotSeg) && isGoodSegment(bestTopSeg)) {
0245       if (detailedAnalysis)
0246         histos[4]->Fill(posAtMid.x(), posAtMid.y());
0247       //check if interpolation fall inside middle chamber
0248       if ((dtGeom->chamber(MidId))->surface().bounds().inside(posAtMid)) {
0249         histos[0]->Fill(posAtMid.x(), posAtMid.y());
0250         if (nSegsMid > 0) {
0251           if (detailedAnalysis) {
0252             histos[3]->Fill(2);
0253             histos[5]->Fill(posAtMid.x(), posAtMid.y());
0254           }
0255
0256           const DTRecSegment4D& bestMidSeg = getBestSegment(segsMid);
0257           // check if middle segments is good enough
0258           if (isGoodSegment(bestMidSeg)) {
0259             if (detailedAnalysis)
0260               histos[6]->Fill(posAtMid.x(), posAtMid.y());
0261             LocalPoint midSegPos = bestMidSeg.localPosition();
0262
0263             // check if middle segments is also close enough
0264             double dist;
0265             if (bestMidSeg.hasPhi()) {
0266               if (bestTopSeg.hasZed() && bestBotSeg.hasZed() && bestMidSeg.hasZed()) {
0267                 dist = (midSegPos - posAtMid).mag();
0268               } else {
0269                 dist = fabs((midSegPos - posAtMid).x());
0270               }
0271             } else {
0272               dist = fabs((midSegPos - posAtMid).y());
0273             }
0274             if (dist < theMinCloseDist) {
0275               histos[1]->Fill(posAtMid.x(), posAtMid.y());
0276             }
0277             if (detailedAnalysis)
0278               histos[2]->Fill(dist);
0279           }
0280         }
0281       }
0282     }
0283   }  // loop over stations
0284 }
0285
0286 // requirements : max number of hits and min chi2
0287 const DTRecSegment4D& DTChamberEfficiencyTask::getBestSegment(const DTRecSegment4DCollection::range& segs) const {
0288   DTRecSegment4DCollection::const_iterator bestIter;
0289   unsigned int nHitBest = 0;
0290   double chi2Best = 99999.;
0291   for (DTRecSegment4DCollection::const_iterator seg = segs.first; seg != segs.second; ++seg) {
0292     unsigned int nHits = ((*seg).hasPhi() ? (*seg).phiSegment()->recHits().size() : 0);
0293     nHits += ((*seg).hasZed() ? (*seg).zSegment()->recHits().size() : 0);
0294
0295     if (nHits == nHitBest) {
0296       if ((*seg).chi2() / (*seg).degreesOfFreedom() < chi2Best) {
0297         chi2Best = (*seg).chi2() / (*seg).degreesOfFreedom();
0298         bestIter = seg;
0299       }
0300     } else if (nHits > nHitBest) {
0301       nHitBest = nHits;
0302       bestIter = seg;
0303     }
0304   }
0305   return *bestIter;
0306 }
0307
0308 const DTRecSegment4D* DTChamberEfficiencyTask::getBestSegment(const DTRecSegment4D* s1,
0309                                                               const DTRecSegment4D* s2) const {
0310   if (!s1)
0311     return s2;
0312   if (!s2)
0313     return s1;
0314   unsigned int nHits1 = (s1->hasPhi() ? s1->phiSegment()->recHits().size() : 0);
0315   nHits1 += (s1->hasZed() ? s1->zSegment()->recHits().size() : 0);
0316
0317   unsigned int nHits2 = (s2->hasPhi() ? s2->phiSegment()->recHits().size() : 0);
0318   nHits2 += (s2->hasZed() ? s2->zSegment()->recHits().size() : 0);
0319
0320   if (nHits1 == nHits2) {
0321     if (s1->chi2() / s1->degreesOfFreedom() < s2->chi2() / s2->degreesOfFreedom())
0322       return s1;
0323     else
0324       return s2;
0325   } else if (nHits1 > nHits2)
0326     return s1;
0327   return s2;
0328 }
0329
0330 LocalPoint DTChamberEfficiencyTask::interpolate(const DTRecSegment4D& seg1,
0331                                                 const DTRecSegment4D& seg3,
0332                                                 const DTChamberId& id2) const {
0333   // Get GlobalPoition of Seg in MB1
0334   GlobalPoint gpos1 = (dtGeom->chamber(seg1.chamberId()))->toGlobal(seg1.localPosition());
0335
0336   // Get GlobalPoition of Seg in MB3
0337   GlobalPoint gpos3 = (dtGeom->chamber(seg3.chamberId()))->toGlobal(seg3.localPosition());
0338
0339   // interpolate
0340   // get all in MB2 frame
0341   LocalPoint pos1 = (dtGeom->chamber(id2))->toLocal(gpos1);
0342   LocalPoint pos3 = (dtGeom->chamber(id2))->toLocal(gpos3);
0343
0344   // case 1: 1 and 3 has both projection. No problem
0345
0346   // case 2: one projection is missing for one of the segments. Keep the other's segment position
0347   if (!seg1.hasZed())
0348     pos1 = LocalPoint(pos1.x(), pos3.y(), pos1.z());
0349   if (!seg3.hasZed())
0350     pos3 = LocalPoint(pos3.x(), pos1.y(), pos3.z());
0351
0352   if (!seg1.hasPhi())
0353     pos1 = LocalPoint(pos3.x(), pos1.y(), pos1.z());
0354   if (!seg3.hasPhi())
0355     pos3 = LocalPoint(pos1.x(), pos3.y(), pos3.z());
0356
0357   // direction
0358   LocalVector dir = (pos3 - pos1).unit();  // z points inward!
0359   LocalPoint pos2 = pos1 + dir * pos1.z() / (-dir.z());
0360
0361   return pos2;
0362 }
0363
0364 bool DTChamberEfficiencyTask::isGoodSegment(const DTRecSegment4D& seg) const {
0365   if (seg.chamberId().station() != 4 && !seg.hasZed())
0366     return false;
0367   unsigned int nHits = (seg.hasPhi() ? seg.phiSegment()->recHits().size() : 0);
0368   nHits += (seg.hasZed() ? seg.zSegment()->recHits().size() : 0);
0369   return (nHits >= theMinHitsSegment && seg.chi2() / seg.degreesOfFreedom() < theMinChi2NormSegment);
0370 }
0371
0372 // Local Variables:
0373 // show-trailing-whitespace: t
0374 // truncate-lines: t
0375 // End: