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/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 "DTSegment4DQuality.h"
0022 #include "Histograms.h"
0023
0024 using namespace std;
0025 using namespace edm;
0026
0027 namespace dtsegment4d {
0028 struct Histograms {
0029 std::unique_ptr<HRes4DHit> h4DHit;
0030 std::unique_ptr<HRes4DHit> h4DHit_W0;
0031 std::unique_ptr<HRes4DHit> h4DHit_W1;
0032 std::unique_ptr<HRes4DHit> h4DHit_W2;
0033 std::unique_ptr<HRes4DHit> h4DHitWS[3][4];
0034
0035 std::unique_ptr<HEff4DHit> hEff_All;
0036 std::unique_ptr<HEff4DHit> hEff_W0;
0037 std::unique_ptr<HEff4DHit> hEff_W1;
0038 std::unique_ptr<HEff4DHit> hEff_W2;
0039 std::unique_ptr<HEff4DHit> hEffWS[3][4];
0040 };
0041 }
0042 using namespace dtsegment4d;
0043
0044
0045
0046
0047
0048
0049
0050 namespace {
0051 constexpr bool mirrorMinusWheels = true;
0052 }
0053
0054
0055 DTSegment4DQuality::DTSegment4DQuality(const ParameterSet &pset) : muonGeomToken_(esConsumes()) {
0056
0057 debug_ = pset.getUntrackedParameter<bool>("debug");
0058 DTHitQualityUtils::debug = debug_;
0059
0060
0061 simHitLabel_ = pset.getUntrackedParameter<InputTag>("simHitLabel");
0062 simHitToken_ = consumes<PSimHitContainer>(pset.getUntrackedParameter<InputTag>("simHitLabel"));
0063
0064 segment4DLabel_ = pset.getUntrackedParameter<InputTag>("segment4DLabel");
0065 segment4DToken_ = consumes<DTRecSegment4DCollection>(pset.getUntrackedParameter<InputTag>("segment4DLabel"));
0066
0067
0068 sigmaResX_ = pset.getParameter<double>("sigmaResX");
0069 sigmaResY_ = pset.getParameter<double>("sigmaResY");
0070
0071 sigmaResAlpha_ = pset.getParameter<double>("sigmaResAlpha");
0072 sigmaResBeta_ = pset.getParameter<double>("sigmaResBeta");
0073 doall_ = pset.getUntrackedParameter<bool>("doall", false);
0074 local_ = pset.getUntrackedParameter<bool>("local", false);
0075 }
0076
0077 void DTSegment4DQuality::bookHistograms(DQMStore::IBooker &booker,
0078 edm::Run const &run,
0079 edm::EventSetup const &setup,
0080 Histograms &histograms) const {
0081 histograms.h4DHit = std::make_unique<HRes4DHit>("All", booker, doall_, local_);
0082 histograms.h4DHit_W0 = std::make_unique<HRes4DHit>("W0", booker, doall_, local_);
0083 histograms.h4DHit_W1 = std::make_unique<HRes4DHit>("W1", booker, doall_, local_);
0084 histograms.h4DHit_W2 = std::make_unique<HRes4DHit>("W2", booker, doall_, local_);
0085
0086 if (doall_) {
0087 histograms.hEff_All = std::make_unique<HEff4DHit>("All", booker);
0088 histograms.hEff_W0 = std::make_unique<HEff4DHit>("W0", booker);
0089 histograms.hEff_W1 = std::make_unique<HEff4DHit>("W1", booker);
0090 histograms.hEff_W2 = std::make_unique<HEff4DHit>("W2", booker);
0091 }
0092
0093 if (local_) {
0094
0095 TString name = "W";
0096 for (long w = 0; w <= 2; ++w) {
0097 for (long s = 1; s <= 4; ++s) {
0098
0099 TString nameWS = (name + w + "_St" + s);
0100 histograms.h4DHitWS[w][s - 1] = std::make_unique<HRes4DHit>(nameWS.Data(), booker, doall_, local_);
0101 histograms.hEffWS[w][s - 1] = std::make_unique<HEff4DHit>(nameWS.Data(), booker);
0102 }
0103 }
0104 }
0105 };
0106
0107
0108 void DTSegment4DQuality::dqmAnalyze(edm::Event const &event,
0109 edm::EventSetup const &setup,
0110 Histograms const &histograms) const {
0111 const float epsilon = 5e-5;
0112
0113
0114 const DTGeometry &dtGeom = setup.getData(muonGeomToken_);
0115
0116
0117 edm::Handle<PSimHitContainer> simHits;
0118 event.getByToken(simHitToken_, simHits);
0119
0120
0121 map<DTChamberId, PSimHitContainer> simHitsPerCh;
0122 for (const auto &simHit : *simHits) {
0123
0124
0125 if (abs(simHit.particleType()) == 13) {
0126
0127 DTChamberId chamberId = (((DTWireId(simHit.detUnitId())).layerId()).superlayerId()).chamberId();
0128
0129 simHitsPerCh[chamberId].push_back(simHit);
0130 }
0131 }
0132
0133
0134 Handle<DTRecSegment4DCollection> segment4Ds;
0135 event.getByToken(segment4DToken_, segment4Ds);
0136
0137 if (!segment4Ds.isValid()) {
0138 if (debug_) {
0139 cout << "[DTSegment4DQuality]**Warning: no 4D Segments with label: " << segment4DLabel_
0140 << " in this event, skipping!" << endl;
0141 }
0142 return;
0143 }
0144
0145
0146 for (auto &simHitsInChamber : simHitsPerCh) {
0147 DTChamberId chamberId = simHitsInChamber.first;
0148 int station = chamberId.station();
0149 if (station == 4 && !(local_)) {
0150 continue;
0151 }
0152 int wheel = chamberId.wheel();
0153
0154
0155
0156 const PSimHitContainer &simHits = simHitsInChamber.second;
0157
0158
0159 auto const &simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(simHits);
0160 auto const &muSimHitPerWire = DTHitQualityUtils::mapMuSimHitsPerWire(simHitsPerWire);
0161 int nMuSimHit = muSimHitPerWire.size();
0162 if (nMuSimHit < 2) {
0163 continue;
0164 }
0165 if (debug_) {
0166 cout << "=== Chamber " << chamberId << " has " << nMuSimHit << " SimHits" << endl;
0167 }
0168
0169
0170 pair<const PSimHit *, const PSimHit *> inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire);
0171
0172
0173 if ((DTWireId(inAndOutSimHit.first->detUnitId())).superlayer() ==
0174 (DTWireId(inAndOutSimHit.second->detUnitId())).superLayer()) {
0175 continue;
0176 }
0177
0178
0179 pair<LocalVector, LocalPoint> dirAndPosSimSegm =
0180 DTHitQualityUtils::findMuSimSegmentDirAndPos(inAndOutSimHit, chamberId, dtGeom);
0181
0182 LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
0183 LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
0184 const DTChamber *chamber = dtGeom.chamber(chamberId);
0185 GlobalPoint simSegmGlobalPos = chamber->toGlobal(simSegmLocalPos);
0186 GlobalVector simSegmGlobalDir = chamber->toGlobal(simSegmLocalDir);
0187
0188
0189 float alphaSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).first;
0190 float betaSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).second;
0191
0192 float xSimSeg = simSegmLocalPos.x();
0193 float ySimSeg = simSegmLocalPos.y();
0194
0195 float etaSimSeg = simSegmGlobalPos.eta();
0196 float phiSimSeg = simSegmGlobalPos.phi();
0197
0198 double count_seg = 0;
0199
0200 if (debug_) {
0201 cout << " Simulated segment: local direction " << simSegmLocalDir << endl
0202 << " local position " << simSegmLocalPos << endl
0203 << " alpha " << alphaSimSeg << endl
0204 << " beta " << betaSimSeg << endl;
0205 }
0206
0207
0208
0209 bool recHitFound = false;
0210 DTRecSegment4DCollection::range range = segment4Ds->get(chamberId);
0211 int nsegm = distance(range.first, range.second);
0212 if (debug_) {
0213 cout << " Chamber: " << chamberId << " has " << nsegm << " 4D segments" << endl;
0214 }
0215
0216 if (nsegm != 0) {
0217
0218
0219
0220
0221
0222 const DTRecSegment4D *bestRecHit = nullptr;
0223 double deltaAlpha = 99999;
0224 double deltaBeta = 99999;
0225
0226
0227 for (DTRecSegment4DCollection::const_iterator segment4D = range.first; segment4D != range.second; ++segment4D) {
0228
0229 if (station != 4 && (*segment4D).dimension() != 4) {
0230 continue;
0231 }
0232
0233 LocalVector recSegDirection = (*segment4D).localDirection();
0234 LocalPoint recSegPosition = (*segment4D).localPosition();
0235
0236 pair<double, double> ab = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection);
0237 float recSegAlpha = ab.first;
0238 float recSegBeta = ab.second;
0239
0240 if (debug_) {
0241 cout << &(*segment4D) << " RecSegment direction: " << recSegDirection << endl
0242 << " position : " << (*segment4D).localPosition() << endl
0243 << " alpha : " << recSegAlpha << endl
0244 << " beta : " << recSegBeta << endl
0245 << " nhits : " << (*segment4D).phiSegment()->recHits().size() << " "
0246 << (((*segment4D).zSegment() != nullptr) ? (*segment4D).zSegment()->recHits().size() : 0) << endl;
0247 }
0248
0249 float dAlphaRecSim = fabs(recSegAlpha - alphaSimSeg);
0250 float dBetaRecSim = fabs(recSegBeta - betaSimSeg);
0251
0252 if ((fabs(recSegPosition.x() - simSegmLocalPos.x()) <
0253 4)
0254 && ((fabs(recSegPosition.y() - simSegmLocalPos.y()) < 4) ||
0255 (*segment4D).dimension() < 4)) {
0256 ++count_seg;
0257
0258 if (fabs(dAlphaRecSim - deltaAlpha) < epsilon) {
0259 if (dBetaRecSim < deltaBeta) {
0260 deltaAlpha = dAlphaRecSim;
0261 deltaBeta = dBetaRecSim;
0262 bestRecHit = &(*segment4D);
0263 }
0264
0265 } else if (dAlphaRecSim < deltaAlpha) {
0266 deltaAlpha = dAlphaRecSim;
0267 deltaBeta = dBetaRecSim;
0268 bestRecHit = &(*segment4D);
0269 }
0270 }
0271
0272 }
0273
0274 if (bestRecHit) {
0275 if (debug_) {
0276 cout << endl << "Chosen: " << bestRecHit << endl;
0277 }
0278
0279 LocalPoint bestRecHitLocalPos = bestRecHit->localPosition();
0280 LocalVector bestRecHitLocalDir = bestRecHit->localDirection();
0281
0282 LocalError bestRecHitLocalPosErr = bestRecHit->localPositionError();
0283 LocalError bestRecHitLocalDirErr = bestRecHit->localDirectionError();
0284
0285 pair<double, double> ab = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir);
0286 float alphaBestRHit = ab.first;
0287 float betaBestRHit = ab.second;
0288
0289
0290
0291
0292
0293
0294 const DTSLRecSegment2D *zedRecSeg = bestRecHit->zSegment();
0295 LocalPoint bestRecHitLocalPosRZ;
0296 LocalVector bestRecHitLocalDirRZ;
0297 LocalError bestRecHitLocalPosErrRZ;
0298 LocalError bestRecHitLocalDirErrRZ;
0299 LocalPoint simSegLocalPosRZ;
0300
0301 float alphaBestRHitRZ = 0;
0302 float alphaSimSegRZ = betaSimSeg;
0303 if (zedRecSeg) {
0304 bestRecHitLocalPosRZ = zedRecSeg->localPosition();
0305 bestRecHitLocalDirRZ = zedRecSeg->localDirection();
0306
0307 bestRecHitLocalPosErrRZ = zedRecSeg->localPositionError();
0308 bestRecHitLocalDirErrRZ = zedRecSeg->localDirectionError();
0309
0310
0311 alphaBestRHitRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDirRZ).first;
0312
0313
0314 const DTSuperLayer *sl = dtGeom.superLayer(zedRecSeg->superLayerId());
0315 LocalPoint simSegLocalPosRZTmp = sl->toLocal(simSegmGlobalPos);
0316 LocalVector simSegLocalDirRZ = sl->toLocal(simSegmGlobalDir);
0317 simSegLocalPosRZ =
0318 simSegLocalPosRZTmp + simSegLocalDirRZ * (-simSegLocalPosRZTmp.z() / (cos(simSegLocalDirRZ.theta())));
0319 alphaSimSegRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegLocalDirRZ).first;
0320
0321 if (debug_) {
0322 cout << "RZ SL: recPos " << bestRecHitLocalPosRZ << "recDir " << bestRecHitLocalDirRZ << "recAlpha "
0323 << alphaBestRHitRZ << endl
0324 << "RZ SL: simPos " << simSegLocalPosRZ << "simDir " << simSegLocalDirRZ << "simAlpha "
0325 << alphaSimSegRZ << endl;
0326 }
0327 }
0328
0329
0330 const DTChamberRecSegment2D *phiSeg = bestRecHit->phiSegment();
0331
0332 float t0phi = -999;
0333 float t0theta = -999;
0334 int nHitPhi = 0;
0335 int nHitTheta = 0;
0336
0337 if (phiSeg) {
0338 t0phi = phiSeg->t0();
0339 nHitPhi = phiSeg->recHits().size();
0340 }
0341
0342 if (zedRecSeg) {
0343 t0theta = zedRecSeg->t0();
0344 nHitTheta = zedRecSeg->recHits().size();
0345 }
0346
0347 recHitFound = true;
0348
0349
0350
0351 if (mirrorMinusWheels && wheel < 0) {
0352 alphaSimSeg *= -1.;
0353 alphaBestRHit *= -1.;
0354
0355
0356
0357 }
0358
0359
0360 HRes4DHit *histo = nullptr;
0361
0362 if (wheel == 0) {
0363 histo = histograms.h4DHit_W0.get();
0364 } else if (abs(wheel) == 1) {
0365 histo = histograms.h4DHit_W1.get();
0366 } else if (abs(wheel) == 2) {
0367 histo = histograms.h4DHit_W2.get();
0368 }
0369
0370 float sigmaAlphaBestRhit = sqrt(DTHitQualityUtils::sigmaAngle(alphaBestRHit, bestRecHitLocalDirErr.xx()));
0371 float sigmaBetaBestRhit =
0372 sqrt(DTHitQualityUtils::sigmaAngle(betaBestRHit,
0373 bestRecHitLocalDirErr.yy()));
0374
0375 float sigmaAlphaBestRhitRZ = sqrt(DTHitQualityUtils::sigmaAngle(alphaBestRHitRZ, bestRecHitLocalDirErrRZ.xx()));
0376
0377 histo->fill(alphaSimSeg,
0378 alphaBestRHit,
0379 betaSimSeg,
0380 betaBestRHit,
0381 xSimSeg,
0382 bestRecHitLocalPos.x(),
0383 ySimSeg,
0384 bestRecHitLocalPos.y(),
0385 etaSimSeg,
0386 phiSimSeg,
0387 bestRecHitLocalPosRZ.x(),
0388 simSegLocalPosRZ.x(),
0389 alphaBestRHitRZ,
0390 alphaSimSegRZ,
0391 sigmaAlphaBestRhit,
0392 sigmaBetaBestRhit,
0393 sqrt(bestRecHitLocalPosErr.xx()),
0394 sqrt(bestRecHitLocalPosErr.yy()),
0395 sigmaAlphaBestRhitRZ,
0396 sqrt(bestRecHitLocalPosErrRZ.xx()),
0397 nHitPhi,
0398 nHitTheta,
0399 t0phi,
0400 t0theta);
0401
0402 histograms.h4DHit->fill(alphaSimSeg,
0403 alphaBestRHit,
0404 betaSimSeg,
0405 betaBestRHit,
0406 xSimSeg,
0407 bestRecHitLocalPos.x(),
0408 ySimSeg,
0409 bestRecHitLocalPos.y(),
0410 etaSimSeg,
0411 phiSimSeg,
0412 bestRecHitLocalPosRZ.x(),
0413 simSegLocalPosRZ.x(),
0414 alphaBestRHitRZ,
0415 alphaSimSegRZ,
0416 sigmaAlphaBestRhit,
0417 sigmaBetaBestRhit,
0418 sqrt(bestRecHitLocalPosErr.xx()),
0419 sqrt(bestRecHitLocalPosErr.yy()),
0420 sigmaAlphaBestRhitRZ,
0421 sqrt(bestRecHitLocalPosErrRZ.xx()),
0422 nHitPhi,
0423 nHitTheta,
0424 t0phi,
0425 t0theta);
0426
0427 if (local_) {
0428 histograms.h4DHitWS[abs(wheel)][station - 1]->fill(alphaSimSeg,
0429 alphaBestRHit,
0430 betaSimSeg,
0431 betaBestRHit,
0432 xSimSeg,
0433 bestRecHitLocalPos.x(),
0434 ySimSeg,
0435 bestRecHitLocalPos.y(),
0436 etaSimSeg,
0437 phiSimSeg,
0438 bestRecHitLocalPosRZ.x(),
0439 simSegLocalPosRZ.x(),
0440 alphaBestRHitRZ,
0441 alphaSimSegRZ,
0442 sigmaAlphaBestRhit,
0443 sigmaBetaBestRhit,
0444 sqrt(bestRecHitLocalPosErr.xx()),
0445 sqrt(bestRecHitLocalPosErr.yy()),
0446 sigmaAlphaBestRhitRZ,
0447 sqrt(bestRecHitLocalPosErrRZ.xx()),
0448 nHitPhi,
0449 nHitTheta,
0450 t0phi,
0451 t0theta);
0452 }
0453
0454 }
0455
0456 }
0457
0458
0459 if (doall_) {
0460 HEff4DHit *heff = nullptr;
0461
0462 if (wheel == 0) {
0463 heff = histograms.hEff_W0.get();
0464 } else if (abs(wheel) == 1) {
0465 heff = histograms.hEff_W1.get();
0466 } else if (abs(wheel) == 2) {
0467 heff = histograms.hEff_W2.get();
0468 }
0469 heff->fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound, count_seg);
0470 histograms.hEff_All->fill(
0471 etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound, count_seg);
0472 if (local_) {
0473 histograms.hEffWS[abs(wheel)][station - 1]->fill(
0474 etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound, count_seg);
0475 }
0476 }
0477 }
0478 }
0479
0480
0481 #include "FWCore/Framework/interface/MakerMacros.h"
0482 DEFINE_FWK_MODULE(DTSegment4DQuality);