File indexing completed on 2024-04-06 12:26:10
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include "RecoLocalMuon/DTSegment/test/DTEffAnalyzer.h"
0015
0016
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Framework/interface/ESHandle.h"
0022 #include "FWCore/Utilities/interface/Exception.h"
0023
0024 using namespace edm;
0025
0026 #include "TFile.h"
0027 #include "TH1F.h"
0028 #include "TH2F.h"
0029
0030 #include "DataFormats/LTCDigi/interface/LTCDigi.h"
0031
0032 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0033 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0034 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0035 #include "DataFormats/DTRecHit/interface/DTRangeMapAccessor.h"
0036
0037 #include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
0038 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0039
0040
0041 #include <iostream>
0042 #include <cmath>
0043 using namespace std;
0044
0045
0046
0047
0048 DTEffAnalyzer::DTEffAnalyzer(const ParameterSet& pset) : theDtGeomToken(esConsumes()) {
0049
0050 debug = pset.getUntrackedParameter<bool>("debug");
0051 theRootFileName = pset.getUntrackedParameter<string>("rootFileName");
0052
0053
0054 theRecHits1DLabel = pset.getParameter<string>("recHits1DLabel");
0055
0056
0057 theRecHits2DLabel = pset.getParameter<string>("recHits2DLabel");
0058
0059
0060 theRecHits4DLabel = pset.getParameter<string>("recHits4DLabel");
0061
0062 theMinHitsSegment = static_cast<unsigned int>(pset.getParameter<int>("minHitsSegment"));
0063 theMinChi2NormSegment = pset.getParameter<double>("minChi2NormSegment");
0064 theMinCloseDist = pset.getParameter<double>("minCloseDist");
0065
0066 consumes<DTRecSegment4DCollection>(theRecHits4DLabel);
0067 }
0068
0069 void DTEffAnalyzer::beginJob() {
0070 if (debug)
0071 cout << "beginOfJob" << endl;
0072
0073
0074 theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0075 bool dirStat = TH1::AddDirectoryStatus();
0076 TH1::AddDirectory(kTRUE);
0077
0078
0079 new TH1F("hTrigBits", "All trigger bits", 10, 0., 10.);
0080
0081 for (int w = -2; w <= 2; ++w) {
0082 stringstream nameWheel;
0083 nameWheel << "_Wh" << w;
0084
0085 for (int sec = 1; sec <= 14; ++sec) {
0086 stringstream nameSector;
0087 nameSector << nameWheel.str() << "_Sec" << sec;
0088
0089 for (int st = 1; st <= 4; ++st) {
0090
0091 stringstream nameChamber;
0092 nameChamber << nameSector.str() << "_St" << st;
0093
0094
0095 createTH1F("hDistSegFromExtrap", "Distance segments from extrap position ", nameChamber.str(), 200, 0., 200.);
0096 createTH1F("hNaiveEffSeg", "Naive eff ", nameChamber.str(), 10, 0., 10.);
0097 createTH2F(
0098 "hEffSegVsPosDen", "Eff vs local position (all) ", nameChamber.str(), 25, -250., 250., 25, -250., 250.);
0099 createTH2F(
0100 "hEffGoodSegVsPosDen", "Eff vs local position (good) ", nameChamber.str(), 25, -250., 250., 25, -250., 250.);
0101 createTH2F("hEffSegVsPosNum", "Eff vs local position ", nameChamber.str(), 25, -250., 250., 25, -250., 250.);
0102 createTH2F("hEffGoodSegVsPosNum",
0103 "Eff vs local position (good segs) ",
0104 nameChamber.str(),
0105 25,
0106 -250.,
0107 250.,
0108 25,
0109 -250.,
0110 250.);
0111 createTH2F("hEffGoodCloseSegVsPosNum",
0112 "Eff vs local position (good aand close segs) ",
0113 nameChamber.str(),
0114 25,
0115 -250.,
0116 250.,
0117 25,
0118 -250.,
0119 250.);
0120 }
0121 }
0122 }
0123
0124
0125 TH1::AddDirectory(dirStat);
0126 }
0127
0128
0129 DTEffAnalyzer::~DTEffAnalyzer() {
0130 theFile->cd();
0131 theFile->Write();
0132 theFile->Close();
0133 }
0134
0135
0136 void DTEffAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0137 dtGeom = eventSetup.getHandle(theDtGeomToken);
0138
0139 if (debug)
0140 cout << endl
0141 << "--- [DTEffAnalyzer] Event analysed #Run: " << event.id().run() << " #Event: " << event.id().event()
0142 << endl;
0143
0144 effSegments(event, eventSetup);
0145 }
0146
0147 void DTEffAnalyzer::effSegments(const Event& event, const EventSetup& eventSetup) {
0148
0149
0150 event.getByLabel(theRecHits4DLabel, segs);
0151 if (debug) {
0152 cout << "4d " << segs->size() << endl;
0153 for (DTRecSegment4DCollection::const_iterator seg = segs->begin(); seg != segs->end(); ++seg)
0154 cout << *seg << endl;
0155 }
0156
0157
0158
0159
0160
0161
0162
0163 for (int wheel = -2; wheel <= 2; ++wheel) {
0164 for (int sector = 1; sector <= 12; ++sector) {
0165 evaluateEff(DTChamberId(wheel, 1, sector), 2, 3);
0166 evaluateEff(DTChamberId(wheel, 2, sector), 1, 3);
0167 evaluateEff(DTChamberId(wheel, 3, sector), 2, 4);
0168 evaluateEff(DTChamberId(wheel, 4, sector), 2, 3);
0169 }
0170 }
0171 }
0172
0173 void DTEffAnalyzer::evaluateEff(const DTChamberId& MidId, int bottom, int top) const {
0174 if (debug)
0175 cout << "evaluateEff " << MidId << " bott/top " << bottom << "/" << top << endl;
0176
0177 DTChamberId BotId(MidId.wheel(), bottom, MidId.sector());
0178 DTChamberId TopId(MidId.wheel(), top, MidId.sector());
0179
0180
0181 DTRecSegment4DCollection::range segsBot = segs->get(BotId);
0182 int nSegsBot = segsBot.second - segsBot.first;
0183
0184 if (nSegsBot == 0)
0185 return;
0186
0187
0188 DTRecSegment4DCollection::range segsTop = segs->get(TopId);
0189 int nSegsTop = segsTop.second - segsTop.first;
0190
0191
0192 const DTRecSegment4D& bestBotSeg = getBestSegment(segsBot);
0193
0194
0195 DTRecSegment4D* pBestTopSeg = 0;
0196 if (nSegsTop > 0)
0197 pBestTopSeg = const_cast<DTRecSegment4D*>(&getBestSegment(segsTop));
0198
0199 if (TopId.station() == 4 && TopId.sector() == 10) {
0200
0201
0202
0203 DTChamberId TopId14(MidId.wheel(), top, 14);
0204 DTRecSegment4DCollection::range segsTop14 = segs->get(TopId14);
0205 int nSegsTop14 = segsTop14.second - segsTop14.first;
0206
0207 nSegsTop += nSegsTop;
0208 if (nSegsTop14) {
0209 DTRecSegment4D* pBestTopSeg14 = const_cast<DTRecSegment4D*>(&getBestSegment(segsTop14));
0210
0211
0212 pBestTopSeg = const_cast<DTRecSegment4D*>(getBestSegment(pBestTopSeg, pBestTopSeg14));
0213
0214 }
0215 }
0216 if (!pBestTopSeg)
0217 return;
0218 const DTRecSegment4D& bestTopSeg = *pBestTopSeg;
0219
0220
0221 DTRecSegment4DCollection::range segsMid = segs->get(MidId);
0222 int nSegsMid = segsMid.second - segsMid.first;
0223
0224
0225
0226
0227
0228
0229 histo(hName("hNaiveEffSeg", MidId))->Fill(0);
0230 if (nSegsMid > 0)
0231 histo(hName("hNaiveEffSeg", MidId))->Fill(1);
0232
0233
0234
0235 LocalPoint posAtMid = interpolate(bestBotSeg, bestTopSeg, MidId);
0236
0237
0238
0239
0240 if (isGoodSegment(bestBotSeg) && isGoodSegment(bestTopSeg)) {
0241 histo2d(hName("hEffSegVsPosDen", MidId))->Fill(posAtMid.x(), posAtMid.y());
0242
0243 if ((dtGeom->chamber(MidId))->surface().bounds().inside(posAtMid)) {
0244
0245
0246
0247 histo2d(hName("hEffGoodSegVsPosDen", MidId))->Fill(posAtMid.x(), posAtMid.y());
0248
0249 if (nSegsMid > 0) {
0250 histo(hName("hNaiveEffSeg", MidId))->Fill(2);
0251 histo2d(hName("hEffSegVsPosNum", MidId))->Fill(posAtMid.x(), posAtMid.y());
0252 const DTRecSegment4D& bestMidSeg = getBestSegment(segsMid);
0253
0254 if (isGoodSegment(bestMidSeg)) {
0255 histo2d(hName("hEffGoodSegVsPosNum", MidId))->Fill(posAtMid.x(), posAtMid.y());
0256 LocalPoint midSegPos = bestMidSeg.localPosition();
0257
0258 double dist;
0259
0260
0261
0262
0263
0264
0265
0266
0267 if (bestMidSeg.hasPhi()) {
0268 if (bestTopSeg.hasZed() && bestBotSeg.hasZed() && bestMidSeg.hasZed()) {
0269 dist = (midSegPos - posAtMid).mag();
0270 } else {
0271 dist = fabs((midSegPos - posAtMid).x());
0272 }
0273 } else {
0274 dist = fabs((midSegPos - posAtMid).y());
0275 }
0276
0277 if (dist < theMinCloseDist) {
0278 histo2d(hName("hEffGoodCloseSegVsPosNum", MidId))->Fill(posAtMid.x(), posAtMid.y());
0279 }
0280 histo(hName("hDistSegFromExtrap", MidId))->Fill(dist);
0281 }
0282 }
0283 }
0284 }
0285
0286 }
0287
0288
0289 const DTRecSegment4D& DTEffAnalyzer::getBestSegment(const DTRecSegment4DCollection::range& segs) const {
0290 DTRecSegment4DCollection::const_iterator bestIter;
0291 unsigned int nHitBest = 0;
0292 double chi2Best = 99999.;
0293 for (DTRecSegment4DCollection::const_iterator seg = segs.first; seg != segs.second; ++seg) {
0294 unsigned int nHits = ((*seg).hasPhi() ? (*seg).phiSegment()->recHits().size() : 0);
0295 nHits += ((*seg).hasZed() ? (*seg).zSegment()->recHits().size() : 0);
0296
0297 if (nHits == nHitBest) {
0298 if ((*seg).chi2() / (*seg).degreesOfFreedom() < chi2Best) {
0299 chi2Best = (*seg).chi2() / (*seg).degreesOfFreedom();
0300 bestIter = seg;
0301 }
0302 } else if (nHits > nHitBest) {
0303 nHitBest = nHits;
0304 bestIter = seg;
0305 }
0306 }
0307 return *bestIter;
0308 }
0309
0310 const DTRecSegment4D* DTEffAnalyzer::getBestSegment(const DTRecSegment4D* s1, const DTRecSegment4D* s2) const {
0311 if (!s1)
0312 return s2;
0313 if (!s2)
0314 return s1;
0315 unsigned int nHits1 = (s1->hasPhi() ? s1->phiSegment()->recHits().size() : 0);
0316 nHits1 += (s1->hasZed() ? s1->zSegment()->recHits().size() : 0);
0317
0318 unsigned int nHits2 = (s2->hasPhi() ? s2->phiSegment()->recHits().size() : 0);
0319 nHits2 += (s2->hasZed() ? s2->zSegment()->recHits().size() : 0);
0320
0321 if (nHits1 == nHits2) {
0322 if (s1->chi2() / s1->degreesOfFreedom() < s2->chi2() / s2->degreesOfFreedom())
0323 return s1;
0324 else
0325 return s2;
0326 } else if (nHits1 > nHits2)
0327 return s1;
0328 return s2;
0329 }
0330
0331 bool DTEffAnalyzer::isGoodSegment(const DTRecSegment4D& seg) const {
0332 if (seg.chamberId().station() != 4 && !seg.hasZed())
0333 return false;
0334 unsigned int nHits = (seg.hasPhi() ? seg.phiSegment()->recHits().size() : 0);
0335 nHits += (seg.hasZed() ? seg.zSegment()->recHits().size() : 0);
0336 return (nHits >= theMinHitsSegment && seg.chi2() / seg.degreesOfFreedom() < theMinChi2NormSegment);
0337 }
0338
0339 LocalPoint DTEffAnalyzer::interpolate(const DTRecSegment4D& seg1,
0340 const DTRecSegment4D& seg3,
0341 const DTChamberId& id2) const {
0342
0343 GlobalPoint gpos1 = (dtGeom->chamber(seg1.chamberId()))->toGlobal(seg1.localPosition());
0344
0345
0346 GlobalPoint gpos3 = (dtGeom->chamber(seg3.chamberId()))->toGlobal(seg3.localPosition());
0347
0348
0349
0350 LocalPoint pos1 = (dtGeom->chamber(id2))->toLocal(gpos1);
0351 LocalPoint pos3 = (dtGeom->chamber(id2))->toLocal(gpos3);
0352
0353
0354
0355
0356
0357
0358
0359 if (!seg1.hasZed())
0360 pos1 = LocalPoint(pos1.x(), pos3.y(), pos1.z());
0361 if (!seg3.hasZed())
0362 pos3 = LocalPoint(pos3.x(), pos1.y(), pos3.z());
0363
0364 if (!seg1.hasPhi())
0365 pos1 = LocalPoint(pos3.x(), pos1.y(), pos1.z());
0366 if (!seg3.hasPhi())
0367 pos3 = LocalPoint(pos1.x(), pos3.y(), pos3.z());
0368
0369
0370
0371
0372 LocalVector dir = (pos3 - pos1).unit();
0373
0374 LocalPoint pos2 = pos1 + dir * pos1.z() / (-dir.z());
0375
0376
0377 return pos2;
0378 }
0379
0380 TH1F* DTEffAnalyzer::histo(const string& name) const {
0381 if (!theFile->Get(name.c_str()))
0382 throw cms::Exception("DTEffAnalyzer") << " TH1F not existing " << name;
0383 if (TH1F* h = dynamic_cast<TH1F*>(theFile->Get(name.c_str())))
0384 return h;
0385 else
0386 throw cms::Exception("DTEffAnalyzer") << " Not a TH1F " << name;
0387 }
0388
0389 TH2F* DTEffAnalyzer::histo2d(const string& name) const {
0390 if (!theFile->Get(name.c_str()))
0391 throw cms::Exception("DTEffAnalyzer") << " TH1F not existing " << name;
0392 if (TH2F* h = dynamic_cast<TH2F*>(theFile->Get(name.c_str())))
0393 return h;
0394 else
0395 throw cms::Exception("DTEffAnalyzer") << " Not a TH2F " << name;
0396 }
0397
0398 string DTEffAnalyzer::toString(const DTChamberId& id) const {
0399 stringstream result;
0400 result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station();
0401 return result.str();
0402 }
0403
0404 template <class T>
0405 string DTEffAnalyzer::hName(const string& s, const T& id) const {
0406 string name(toString(id));
0407 stringstream hName;
0408 hName << s << name;
0409 return hName.str();
0410 }
0411
0412 void DTEffAnalyzer::createTH1F(const std::string& name,
0413 const std::string& title,
0414 const std::string& suffix,
0415 int nbin,
0416 const double& binMin,
0417 const double& binMax) const {
0418 stringstream hName;
0419 stringstream hTitle;
0420 hName << name << suffix;
0421 hTitle << title << suffix;
0422 new TH1F(hName.str().c_str(), hTitle.str().c_str(), nbin, binMin, binMax);
0423 }
0424
0425 void DTEffAnalyzer::createTH2F(const std::string& name,
0426 const std::string& title,
0427 const std::string& suffix,
0428 int nBinX,
0429 const double& binXMin,
0430 const double& binXMax,
0431 int nBinY,
0432 const double& binYMin,
0433 const double& binYMax) const {
0434 stringstream hName;
0435 stringstream hTitle;
0436 hName << name << suffix;
0437 hTitle << title << suffix;
0438 new TH2F(hName.str().c_str(), hTitle.str().c_str(), nBinX, binXMin, binXMax, nBinY, binYMin, binYMax);
0439 }
0440
0441 DEFINE_FWK_MODULE(DTEffAnalyzer);