File indexing completed on 2024-04-06 12:32:00
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/DTRecSegment4DCollection.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/DTChamber.h"
0017 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0018 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0019 #include "Validation/DTRecHits/interface/DTHitQualityUtils.h"
0020
0021 #include "DTSegment2DSLPhiQuality.h"
0022 #include "Histograms.h"
0023
0024 using namespace std;
0025 using namespace edm;
0026
0027 namespace dtsegment2dsl {
0028 struct Histograms {
0029 std::unique_ptr<HRes2DHit> h2DHitSuperPhi;
0030 std::unique_ptr<HEff2DHit> h2DHitEff_SuperPhi;
0031 };
0032 }
0033
0034 using namespace dtsegment2dsl;
0035
0036
0037 DTSegment2DSLPhiQuality::DTSegment2DSLPhiQuality(const ParameterSet &pset) : muonGeomToken_(esConsumes()) {
0038
0039 debug_ = pset.getUntrackedParameter<bool>("debug");
0040 DTHitQualityUtils::debug = debug_;
0041
0042
0043 simHitLabel_ = pset.getUntrackedParameter<InputTag>("simHitLabel");
0044 simHitToken_ = consumes<PSimHitContainer>(pset.getUntrackedParameter<InputTag>("simHitLabel"));
0045
0046 segment4DLabel_ = pset.getUntrackedParameter<InputTag>("segment4DLabel");
0047 segment4DToken_ = consumes<DTRecSegment4DCollection>(pset.getUntrackedParameter<InputTag>("segment4DLabel"));
0048
0049
0050 sigmaResPos_ = pset.getParameter<double>("sigmaResPos");
0051
0052 sigmaResAngle_ = pset.getParameter<double>("sigmaResAngle");
0053 doall_ = pset.getUntrackedParameter<bool>("doall", false);
0054 local_ = pset.getUntrackedParameter<bool>("local", false);
0055 }
0056
0057 void DTSegment2DSLPhiQuality::bookHistograms(DQMStore::IBooker &booker,
0058 edm::Run const &run,
0059 edm::EventSetup const &setup,
0060 Histograms &histograms) const {
0061
0062 histograms.h2DHitSuperPhi = std::make_unique<HRes2DHit>("SuperPhi", booker, doall_, local_);
0063 if (doall_) {
0064 histograms.h2DHitEff_SuperPhi = std::make_unique<HEff2DHit>("SuperPhi", booker);
0065 }
0066 }
0067
0068
0069 void DTSegment2DSLPhiQuality::dqmAnalyze(edm::Event const &event,
0070 edm::EventSetup const &setup,
0071 Histograms const &histograms) const {
0072
0073 const DTGeometry &dtGeom = setup.getData(muonGeomToken_);
0074
0075
0076 edm::Handle<PSimHitContainer> simHits;
0077 event.getByToken(simHitToken_, simHits);
0078
0079
0080 map<DTChamberId, PSimHitContainer> simHitsPerCh;
0081 for (const auto &simHit : *simHits) {
0082
0083 DTChamberId chamberId = (((DTWireId(simHit.detUnitId())).layerId()).superlayerId()).chamberId();
0084
0085 simHitsPerCh[chamberId].push_back(simHit);
0086 }
0087
0088
0089 Handle<DTRecSegment4DCollection> segment4Ds;
0090 event.getByToken(segment4DToken_, segment4Ds);
0091
0092 if (!segment4Ds.isValid()) {
0093 if (debug_) {
0094 cout << "[DTSegment2DSLPhiQuality]**Warning: no 4D Segments with label: " << segment4DLabel_
0095 << " in this event, skipping!" << endl;
0096 }
0097 return;
0098 }
0099
0100
0101 DTRecSegment4DCollection::id_iterator chamberId;
0102 for (chamberId = segment4Ds->id_begin(); chamberId != segment4Ds->id_end(); ++chamberId) {
0103
0104
0105 PSimHitContainer simHits = simHitsPerCh[(*chamberId)];
0106
0107
0108 map<DTWireId, PSimHitContainer> simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(simHits);
0109 map<DTWireId, const PSimHit *> muSimHitPerWire = DTHitQualityUtils::mapMuSimHitsPerWire(simHitsPerWire);
0110 int nMuSimHit = muSimHitPerWire.size();
0111 if (nMuSimHit == 0 || nMuSimHit == 1) {
0112 if (debug_ && nMuSimHit == 1) {
0113 cout << "[DTSegment2DSLPhiQuality] Only " << nMuSimHit << " mu SimHit in this chamber, skipping!" << endl;
0114 }
0115 continue;
0116 }
0117 if (debug_) {
0118 cout << "=== Chamber " << (*chamberId) << " has " << nMuSimHit << " SimHits" << endl;
0119 }
0120
0121
0122 pair<const PSimHit *, const PSimHit *> inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire);
0123
0124
0125 pair<LocalVector, LocalPoint> dirAndPosSimSegm =
0126 DTHitQualityUtils::findMuSimSegmentDirAndPos(inAndOutSimHit, (*chamberId), dtGeom);
0127
0128 LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
0129 LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
0130 const DTChamber *chamber = dtGeom.chamber(*chamberId);
0131 GlobalPoint simSegmGlobalPos = chamber->toGlobal(simSegmLocalPos);
0132
0133
0134 float angleSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).first;
0135 float posSimSeg = simSegmLocalPos.x();
0136
0137 float etaSimSeg = simSegmGlobalPos.eta();
0138 float phiSimSeg = simSegmGlobalPos.phi();
0139
0140 if (debug_) {
0141 cout << " Simulated segment: local direction " << simSegmLocalDir << endl
0142 << " local position " << simSegmLocalPos << endl
0143 << " angle " << angleSimSeg << endl;
0144 }
0145
0146
0147
0148 bool recHitFound = false;
0149 DTRecSegment4DCollection::range range = segment4Ds->get(*chamberId);
0150 int nsegm = distance(range.first, range.second);
0151 if (debug_) {
0152 cout << " Chamber: " << *chamberId << " has " << nsegm << " 4D segments" << endl;
0153 }
0154
0155 if (nsegm != 0) {
0156
0157
0158
0159
0160
0161 const DTRecSegment2D *bestRecHit = nullptr;
0162 bool bestRecHitFound = false;
0163 double deltaAlpha = 99999;
0164
0165
0166 for (DTRecSegment4DCollection::const_iterator segment4D = range.first; segment4D != range.second; ++segment4D) {
0167
0168 if ((*segment4D).dimension() != 4) {
0169 if (debug_) {
0170 cout << "[DTSegment2DSLPhiQuality]***Error: This is not 4D "
0171 "segment!!!"
0172 << endl;
0173 }
0174 continue;
0175 }
0176
0177
0178 const DTChamberRecSegment2D *phiSegment2D = (*segment4D).phiSegment();
0179 if ((*phiSegment2D).dimension() != 2) {
0180 if (debug_) {
0181 cout << "[DTSegment2DQuality]***Error: This is not 2D segment!!!" << endl;
0182 }
0183 abort();
0184 }
0185
0186
0187 LocalVector recSegDirection = (*phiSegment2D).localDirection();
0188
0189 float recSegAlpha = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).first;
0190 if (debug_) {
0191 cout << " RecSegment direction: " << recSegDirection << endl
0192 << " position : " << (*phiSegment2D).localPosition() << endl
0193 << " alpha : " << recSegAlpha << endl;
0194 }
0195
0196 if (fabs(recSegAlpha - angleSimSeg) < deltaAlpha) {
0197 deltaAlpha = fabs(recSegAlpha - angleSimSeg);
0198 bestRecHit = &(*phiSegment2D);
0199 bestRecHitFound = true;
0200 }
0201 }
0202
0203 if (bestRecHitFound) {
0204
0205 LocalPoint bestRecHitLocalPos = bestRecHit->localPosition();
0206 LocalVector bestRecHitLocalDir = bestRecHit->localDirection();
0207
0208 LocalError bestRecHitLocalPosErr = bestRecHit->localPositionError();
0209 LocalError bestRecHitLocalDirErr = bestRecHit->localDirectionError();
0210
0211 float angleBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).first;
0212 if (fabs(angleBestRHit - angleSimSeg) < 5 * sigmaResAngle_ &&
0213 fabs(bestRecHitLocalPos.x() - posSimSeg) < 5 * sigmaResPos_) {
0214 recHitFound = true;
0215 }
0216
0217
0218 histograms.h2DHitSuperPhi->fill(angleSimSeg,
0219 angleBestRHit,
0220 posSimSeg,
0221 bestRecHitLocalPos.x(),
0222 etaSimSeg,
0223 phiSimSeg,
0224 sqrt(bestRecHitLocalPosErr.xx()),
0225 sqrt(bestRecHitLocalDirErr.xx()));
0226 }
0227 }
0228
0229
0230 if (doall_) {
0231 histograms.h2DHitEff_SuperPhi->fill(etaSimSeg, phiSimSeg, posSimSeg, angleSimSeg, recHitFound);
0232 }
0233 }
0234 }
0235
0236
0237 #include "FWCore/Framework/interface/MakerMacros.h"
0238 DEFINE_FWK_MODULE(DTSegment2DSLPhiQuality);