File indexing completed on 2024-04-06 12:07:08
0001
0002
0003
0004
0005
0006
0007 #include "DTChamberEfficiencyTask.h"
0008
0009
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
0032 recHits4DToken_ =
0033 consumes<DTRecSegment4DCollection>(edm::InputTag(parameters.getUntrackedParameter<string>("recHits4DLabel")));
0034
0035
0036 theMinHitsSegment = static_cast<unsigned int>(parameters.getParameter<int>("minHitsSegment"));
0037 theMinChi2NormSegment = parameters.getParameter<double>("minChi2NormSegment");
0038
0039 theMinCloseDist = parameters.getParameter<double>("minCloseDist");
0040
0041
0042 onlineMonitor = parameters.getUntrackedParameter<bool>("onlineMonitor");
0043
0044
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
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
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
0083 bookHistos(ibooker, (*ch_it)->id());
0084 }
0085 }
0086
0087
0088 void DTChamberEfficiencyTask::bookHistos(DQMStore::IBooker& ibooker, DTChamberId chId) {
0089 edm::LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiencyTask") << " Booking histos for CH : " << chId;
0090
0091
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
0105 vector<MonitorElement*> histos;
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116 histos.push_back(ibooker.book2D(
0117 "hEffGoodSegVsPosDen" + HistoName, "Eff vs local position (good) ", 25, -250., 250., 25, -250., 250.));
0118
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
0131 histos.push_back(ibooker.book1D("hNaiveEffSeg" + HistoName, "Naive eff ", 10, 0., 10.));
0132
0133 histos.push_back(ibooker.book2D(
0134 "hEffSegVsPosDen" + HistoName, "Eff vs local position (all) ", 25, -250., 250., 25, -250., 250.));
0135
0136 histos.push_back(
0137 ibooker.book2D("hEffSegVsPosNum" + HistoName, "Eff vs local position ", 25, -250., 250., 25, -250., 250.));
0138
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
0150 event.getByToken(recHits4DToken_, segs);
0151
0152 int bottom = 0, top = 0;
0153
0154
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
0166 if (station == 1) {
0167 bottom = 2;
0168 top = 3;
0169 }
0170
0171
0172 if (station == 2) {
0173 bottom = 1;
0174 top = 3;
0175 }
0176
0177
0178 if (station == 3) {
0179 bottom = 2;
0180 top = 4;
0181 }
0182
0183
0184 if (station == 4) {
0185 bottom = 2;
0186 top = 3;
0187 }
0188
0189
0190 DTChamberId BotId(wheel, bottom, sector);
0191 DTChamberId TopId(wheel, top, sector);
0192
0193
0194 DTRecSegment4DCollection::range segsBot = segs->get(BotId);
0195 int nSegsBot = segsBot.second - segsBot.first;
0196
0197 if (nSegsBot == 0)
0198 continue;
0199
0200 vector<MonitorElement*> histos = histosPerCh[MidId];
0201
0202
0203 DTRecSegment4DCollection::range segsTop = segs->get(TopId);
0204 int nSegsTop = segsTop.second - segsTop.first;
0205
0206
0207 const DTRecSegment4D& bestBotSeg = getBestSegment(segsBot);
0208
0209
0210 DTRecSegment4D* pBestTopSeg = nullptr;
0211 if (nSegsTop > 0)
0212 pBestTopSeg = const_cast<DTRecSegment4D*>(&getBestSegment(segsTop));
0213
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
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
0234 histos[3]->Fill(0);
0235 if (nSegsMid > 0)
0236 histos[3]->Fill(1);
0237 }
0238
0239
0240
0241 LocalPoint posAtMid = interpolate(bestBotSeg, bestTopSeg, MidId);
0242
0243
0244 if (isGoodSegment(bestBotSeg) && isGoodSegment(bestTopSeg)) {
0245 if (detailedAnalysis)
0246 histos[4]->Fill(posAtMid.x(), posAtMid.y());
0247
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
0258 if (isGoodSegment(bestMidSeg)) {
0259 if (detailedAnalysis)
0260 histos[6]->Fill(posAtMid.x(), posAtMid.y());
0261 LocalPoint midSegPos = bestMidSeg.localPosition();
0262
0263
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 }
0284 }
0285
0286
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
0334 GlobalPoint gpos1 = (dtGeom->chamber(seg1.chamberId()))->toGlobal(seg1.localPosition());
0335
0336
0337 GlobalPoint gpos3 = (dtGeom->chamber(seg3.chamberId()))->toGlobal(seg3.localPosition());
0338
0339
0340
0341 LocalPoint pos1 = (dtGeom->chamber(id2))->toLocal(gpos1);
0342 LocalPoint pos3 = (dtGeom->chamber(id2))->toLocal(gpos3);
0343
0344
0345
0346
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
0358 LocalVector dir = (pos3 - pos1).unit();
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
0373
0374
0375