File indexing completed on 2023-03-17 11:27:05
0001
0002
0003
0004
0005
0006
0007 #include <iostream>
0008 #include <map>
0009
0010 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0011 #include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
0012 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0017 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
0018 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0019 #include "Validation/DTRecHits/interface/DTHitQualityUtils.h"
0020
0021 #include "DTSegment2DQuality.h"
0022 #include "Histograms.h"
0023
0024 using namespace std;
0025 using namespace edm;
0026
0027 namespace dtsegment2d {
0028 struct Histograms {
0029 std::unique_ptr<HRes2DHit> h2DHitRPhi;
0030 std::unique_ptr<HRes2DHit> h2DHitRZ;
0031 std::unique_ptr<HRes2DHit> h2DHitRZ_W0;
0032 std::unique_ptr<HRes2DHit> h2DHitRZ_W1;
0033 std::unique_ptr<HRes2DHit> h2DHitRZ_W2;
0034
0035 std::unique_ptr<HEff2DHit> h2DHitEff_RPhi;
0036 std::unique_ptr<HEff2DHit> h2DHitEff_RZ;
0037 std::unique_ptr<HEff2DHit> h2DHitEff_RZ_W0;
0038 std::unique_ptr<HEff2DHit> h2DHitEff_RZ_W1;
0039 std::unique_ptr<HEff2DHit> h2DHitEff_RZ_W2;
0040 };
0041 }
0042
0043 using namespace dtsegment2d;
0044
0045
0046 DTSegment2DQuality::DTSegment2DQuality(const ParameterSet &pset) : muonGeomToken_(esConsumes()) {
0047
0048 debug_ = pset.getUntrackedParameter<bool>("debug");
0049 DTHitQualityUtils::debug = debug_;
0050
0051 simHitLabel_ = pset.getUntrackedParameter<InputTag>("simHitLabel");
0052 simHitToken_ = consumes<PSimHitContainer>(pset.getUntrackedParameter<InputTag>("simHitLabel"));
0053
0054 segment2DLabel_ = pset.getUntrackedParameter<InputTag>("segment2DLabel");
0055 segment2DToken_ = consumes<DTRecSegment2DCollection>(pset.getUntrackedParameter<InputTag>("segment2DLabel"));
0056
0057
0058 sigmaResPos_ = pset.getParameter<double>("sigmaResPos");
0059
0060 sigmaResAngle_ = pset.getParameter<double>("sigmaResAngle");
0061
0062 if (debug_) {
0063 cout << "[DTSegment2DQuality] Constructor called " << endl;
0064 }
0065 }
0066
0067 void DTSegment2DQuality::bookHistograms(DQMStore::IBooker &booker,
0068 edm::Run const &run,
0069 edm::EventSetup const &setup,
0070 Histograms &histograms) const {
0071 histograms.h2DHitRPhi = std::make_unique<HRes2DHit>("RPhi", booker, true, true);
0072 histograms.h2DHitRZ = std::make_unique<HRes2DHit>("RZ", booker, true, true);
0073 histograms.h2DHitRZ_W0 = std::make_unique<HRes2DHit>("RZ_W0", booker, true, true);
0074 histograms.h2DHitRZ_W1 = std::make_unique<HRes2DHit>("RZ_W1", booker, true, true);
0075 histograms.h2DHitRZ_W2 = std::make_unique<HRes2DHit>("RZ_W2", booker, true, true);
0076
0077 histograms.h2DHitEff_RPhi = std::make_unique<HEff2DHit>("RPhi", booker);
0078 histograms.h2DHitEff_RZ = std::make_unique<HEff2DHit>("RZ", booker);
0079 histograms.h2DHitEff_RZ_W0 = std::make_unique<HEff2DHit>("RZ_W0", booker);
0080 histograms.h2DHitEff_RZ_W1 = std::make_unique<HEff2DHit>("RZ_W1", booker);
0081 histograms.h2DHitEff_RZ_W2 = std::make_unique<HEff2DHit>("RZ_W2", booker);
0082 if (debug_) {
0083 cout << "[DTSegment2DQuality] hitsos created " << endl;
0084 }
0085 }
0086
0087
0088 void DTSegment2DQuality::dqmAnalyze(edm::Event const &event,
0089 edm::EventSetup const &setup,
0090 Histograms const &histograms) const {
0091
0092 const DTGeometry &dtGeom = setup.getData(muonGeomToken_);
0093
0094
0095 edm::Handle<PSimHitContainer> simHits;
0096 event.getByToken(simHitToken_, simHits);
0097
0098
0099 map<DTSuperLayerId, PSimHitContainer> simHitsPerSl;
0100 for (const auto &simHit : *simHits) {
0101
0102 DTSuperLayerId slId = ((DTWireId(simHit.detUnitId())).layerId()).superlayerId();
0103
0104 simHitsPerSl[slId].push_back(simHit);
0105 }
0106
0107
0108 Handle<DTRecSegment2DCollection> segment2Ds;
0109 event.getByToken(segment2DToken_, segment2Ds);
0110
0111 if (not segment2Ds.isValid()) {
0112 if (debug_) {
0113 cout << "[DTSegment2DQuality]**Warning: no 2DSegments with label: " << segment2DLabel_
0114 << " in this event, skipping !" << endl;
0115 }
0116 return;
0117 }
0118
0119
0120 DTRecSegment2DCollection::id_iterator slId;
0121 for (slId = segment2Ds->id_begin(); slId != segment2Ds->id_end(); ++slId) {
0122
0123
0124 PSimHitContainer simHits = simHitsPerSl[(*slId)];
0125
0126
0127 map<DTWireId, PSimHitContainer> simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(simHits);
0128 map<DTWireId, const PSimHit *> muSimHitPerWire = DTHitQualityUtils::mapMuSimHitsPerWire(simHitsPerWire);
0129 int nMuSimHit = muSimHitPerWire.size();
0130 if (nMuSimHit == 0 or nMuSimHit == 1) {
0131 if (debug_ and nMuSimHit == 1) {
0132 cout << "[DTSegment2DQuality] Only " << nMuSimHit << " mu SimHit in this SL, skipping !" << endl;
0133 }
0134 continue;
0135 }
0136 if (debug_) {
0137 cout << "=== SL " << (*slId) << " has " << nMuSimHit << " SimHits" << endl;
0138 }
0139
0140
0141 pair<const PSimHit *, const PSimHit *> inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire);
0142
0143 if (inAndOutSimHit.first == inAndOutSimHit.second) {
0144 cout << "[DTHitQualityUtils]***Warning: outermost and innermost SimHit "
0145 "are the same !"
0146 << endl;
0147 continue;
0148 }
0149
0150
0151 pair<LocalVector, LocalPoint> dirAndPosSimSegm =
0152 DTHitQualityUtils::findMuSimSegmentDirAndPos(inAndOutSimHit, (*slId), dtGeom);
0153
0154 LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
0155 LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
0156 if (debug_) {
0157 cout << " Simulated segment: local direction " << simSegmLocalDir << endl
0158 << " local position " << simSegmLocalPos << endl;
0159 }
0160 const DTSuperLayer *superLayer = dtGeom.superLayer(*slId);
0161 GlobalPoint simSegmGlobalPos = superLayer->toGlobal(simSegmLocalPos);
0162
0163
0164 float angleSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).first;
0165 float posSimSeg = simSegmLocalPos.x();
0166
0167 float etaSimSeg = simSegmGlobalPos.eta();
0168 float phiSimSeg = simSegmGlobalPos.phi();
0169
0170
0171
0172 bool recHitFound = false;
0173 DTRecSegment2DCollection::range range = segment2Ds->get(*slId);
0174 int nsegm = distance(range.first, range.second);
0175 if (debug_) {
0176 cout << " Sl: " << *slId << " has " << nsegm << " 2D segments" << endl;
0177 }
0178
0179 if (nsegm != 0) {
0180
0181
0182
0183
0184
0185 const DTRecSegment2D *bestRecHit = nullptr;
0186 bool bestRecHitFound = false;
0187 double deltaAlpha = 99999;
0188
0189
0190 for (DTRecSegment2DCollection::const_iterator segment2D = range.first; segment2D != range.second; ++segment2D) {
0191
0192 if ((*segment2D).dimension() != 2) {
0193 if (debug_) {
0194 cout << "[DTSegment2DQuality]***Error: This is not 2D segment !!!" << endl;
0195 }
0196 abort();
0197 }
0198
0199 LocalVector recSegDirection = (*segment2D).localDirection();
0200 LocalPoint recSegPosition = (*segment2D).localPosition();
0201
0202 float recSegAlpha = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).first;
0203 if (debug_) {
0204 cout << " RecSegment direction: " << recSegDirection << endl
0205 << " position : " << recSegPosition << endl
0206 << " alpha : " << recSegAlpha << endl;
0207 }
0208
0209 if (fabs(recSegAlpha - angleSimSeg) < deltaAlpha) {
0210 deltaAlpha = fabs(recSegAlpha - angleSimSeg);
0211 bestRecHit = &(*segment2D);
0212 bestRecHitFound = true;
0213 }
0214 }
0215
0216 if (bestRecHitFound) {
0217
0218 LocalPoint bestRecHitLocalPos = bestRecHit->localPosition();
0219 LocalVector bestRecHitLocalDir = bestRecHit->localDirection();
0220
0221 LocalError bestRecHitLocalPosErr = bestRecHit->localPositionError();
0222 LocalError bestRecHitLocalDirErr = bestRecHit->localDirectionError();
0223
0224 float angleBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).first;
0225
0226 if (fabs(angleBestRHit - angleSimSeg) < 5 * sigmaResAngle_ and
0227 fabs(bestRecHitLocalPos.x() - posSimSeg) < 5 * sigmaResPos_) {
0228 recHitFound = true;
0229 }
0230
0231
0232 HRes2DHit *hRes = nullptr;
0233 if ((*slId).superlayer() == 1 or (*slId).superlayer() == 3) {
0234 hRes = histograms.h2DHitRPhi.get();
0235 } else if ((*slId).superlayer() == 2) {
0236 histograms.h2DHitRZ->fill(angleSimSeg,
0237 angleBestRHit,
0238 posSimSeg,
0239 bestRecHitLocalPos.x(),
0240 etaSimSeg,
0241 phiSimSeg,
0242 sqrt(bestRecHitLocalPosErr.xx()),
0243 sqrt(bestRecHitLocalDirErr.xx()));
0244 if (abs((*slId).wheel()) == 0) {
0245 hRes = histograms.h2DHitRZ_W0.get();
0246 } else if (abs((*slId).wheel()) == 1) {
0247 hRes = histograms.h2DHitRZ_W1.get();
0248 } else if (abs((*slId).wheel()) == 2) {
0249 hRes = histograms.h2DHitRZ_W2.get();
0250 }
0251 }
0252 hRes->fill(angleSimSeg,
0253 angleBestRHit,
0254 posSimSeg,
0255 bestRecHitLocalPos.x(),
0256 etaSimSeg,
0257 phiSimSeg,
0258 sqrt(bestRecHitLocalPosErr.xx()),
0259 sqrt(bestRecHitLocalDirErr.xx()));
0260 }
0261 }
0262
0263
0264 HEff2DHit *hEff = nullptr;
0265 if ((*slId).superlayer() == 1 or (*slId).superlayer() == 3) {
0266 hEff = histograms.h2DHitEff_RPhi.get();
0267 } else if ((*slId).superlayer() == 2) {
0268 histograms.h2DHitEff_RZ->fill(etaSimSeg, phiSimSeg, posSimSeg, angleSimSeg, recHitFound);
0269 if (abs((*slId).wheel()) == 0) {
0270 hEff = histograms.h2DHitEff_RZ_W0.get();
0271 } else if (abs((*slId).wheel()) == 1) {
0272 hEff = histograms.h2DHitEff_RZ_W1.get();
0273 } else if (abs((*slId).wheel()) == 2) {
0274 hEff = histograms.h2DHitEff_RZ_W2.get();
0275 }
0276 }
0277 hEff->fill(etaSimSeg, phiSimSeg, posSimSeg, angleSimSeg, recHitFound);
0278 }
0279 }
0280
0281
0282 #include "FWCore/Framework/interface/MakerMacros.h"
0283 DEFINE_FWK_MODULE(DTSegment2DQuality);