File indexing completed on 2024-04-06 12:27:06
0001 #include "RecoMuon/MuonSeedGenerator/plugins/CosmicMuonSeedGenerator.h"
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0011 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0012 #include "DataFormats/Common/interface/Handle.h"
0013
0014 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0015
0016 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0017
0018 #include "FWCore/Framework/interface/EventSetup.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/MessageLogger/interface/MessageLogger.h"
0023
0024 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0025 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0026 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0027
0028 #include "DataFormats/Math/interface/deltaR.h"
0029
0030 #include <vector>
0031
0032 typedef MuonTransientTrackingRecHit::MuonRecHitPointer MuonRecHitPointer;
0033 typedef MuonTransientTrackingRecHit::MuonRecHitContainer MuonRecHitContainer;
0034
0035 using namespace edm;
0036 using namespace std;
0037
0038
0039 CosmicMuonSeedGenerator::CosmicMuonSeedGenerator(const edm::ParameterSet& pset) {
0040 produces<TrajectorySeedCollection>();
0041
0042
0043 theEnableDTFlag = pset.getParameter<bool>("EnableDTMeasurement");
0044
0045 theEnableCSCFlag = pset.getParameter<bool>("EnableCSCMeasurement");
0046
0047 theDTRecSegmentLabel = pset.getParameter<InputTag>("DTRecSegmentLabel");
0048
0049 theCSCRecSegmentLabel = pset.getParameter<InputTag>("CSCRecSegmentLabel");
0050
0051
0052 theMaxSeeds = pset.getParameter<int>("MaxSeeds");
0053
0054 theMaxDTChi2 = pset.getParameter<double>("MaxDTChi2");
0055 theMaxCSCChi2 = pset.getParameter<double>("MaxCSCChi2");
0056
0057 theForcePointDownFlag = pset.existsAs<bool>("ForcePointDown") ? pset.getParameter<bool>("ForcePointDown") : true;
0058
0059
0060 theParameters["topmb41"] = 0.87;
0061 theParameters["bottommb41"] = 1.2;
0062 theParameters["topmb42"] = 0.67;
0063 theParameters["bottommb42"] = 0.98;
0064 theParameters["topmb43"] = 0.34;
0065 theParameters["bottommb43"] = 0.58;
0066 theParameters["topmb31"] = 0.54;
0067 theParameters["bottommb31"] = 0.77;
0068 theParameters["topmb32"] = 0.35;
0069 theParameters["bottommb32"] = 0.55;
0070 theParameters["topmb21"] = 0.21;
0071 theParameters["bottommb21"] = 0.31;
0072
0073 edm::ConsumesCollector iC = consumesCollector();
0074 muonMeasurements = new MuonDetLayerMeasurements(theDTRecSegmentLabel,
0075 theCSCRecSegmentLabel,
0076
0077 InputTag(),
0078 InputTag(),
0079 InputTag(),
0080 iC,
0081 theEnableDTFlag,
0082 theEnableCSCFlag,
0083 false,
0084 false,
0085 false);
0086 muonLayersToken = esConsumes<MuonDetLayerGeometry, MuonRecoGeometryRecord>();
0087 magFieldToken = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0088 }
0089
0090
0091 CosmicMuonSeedGenerator::~CosmicMuonSeedGenerator() {
0092 if (muonMeasurements)
0093 delete muonMeasurements;
0094 }
0095
0096
0097 void CosmicMuonSeedGenerator::produce(edm::Event& event, const edm::EventSetup& eSetup) {
0098 theField = eSetup.getHandle(magFieldToken);
0099
0100 auto output = std::make_unique<TrajectorySeedCollection>();
0101
0102 TrajectorySeedCollection seeds;
0103
0104 std::string category = "Muon|RecoMuon|CosmicMuonSeedGenerator";
0105
0106
0107 theMuonLayers = eSetup.getHandle(muonLayersToken);
0108
0109
0110 vector<const DetLayer*> dtLayers = theMuonLayers->allDTLayers();
0111
0112
0113 vector<const DetLayer*> cscForwardLayers = theMuonLayers->forwardCSCLayers();
0114 vector<const DetLayer*> cscBackwardLayers = theMuonLayers->backwardCSCLayers();
0115
0116 muonMeasurements->setEvent(event);
0117
0118 MuonRecHitContainer allHits;
0119
0120 vector<MuonRecHitContainer> RHMBs;
0121 vector<MuonRecHitContainer> RHMEFs;
0122 vector<MuonRecHitContainer> RHMEBs;
0123
0124 stable_sort(allHits.begin(), allHits.end(), DecreasingGlobalY());
0125
0126 for (vector<const DetLayer*>::reverse_iterator icsclayer = cscForwardLayers.rbegin();
0127 icsclayer != cscForwardLayers.rend() - 1;
0128 ++icsclayer) {
0129 MuonRecHitContainer RHMF = muonMeasurements->recHits(*icsclayer);
0130 allHits.insert(allHits.end(), RHMF.begin(), RHMF.end());
0131 }
0132
0133 for (vector<const DetLayer*>::reverse_iterator icsclayer = cscBackwardLayers.rbegin();
0134 icsclayer != cscBackwardLayers.rend() - 1;
0135 ++icsclayer) {
0136 MuonRecHitContainer RHMF = muonMeasurements->recHits(*icsclayer);
0137 allHits.insert(allHits.end(), RHMF.begin(), RHMF.end());
0138 }
0139
0140 for (vector<const DetLayer*>::reverse_iterator idtlayer = dtLayers.rbegin(); idtlayer != dtLayers.rend();
0141 ++idtlayer) {
0142 MuonRecHitContainer RHMB = muonMeasurements->recHits(*idtlayer);
0143 RHMBs.push_back(RHMB);
0144
0145 if (idtlayer != dtLayers.rbegin())
0146 allHits.insert(allHits.end(), RHMB.begin(), RHMB.end());
0147 }
0148
0149
0150
0151 LogTrace(category) << "all RecHits: " << allHits.size();
0152
0153
0154
0155
0156
0157
0158
0159 CosmicMuonSeedGenerator::MuonRecHitPairVector mb42 = makeSegPairs(RHMBs[0], RHMBs[2], "mb42");
0160 createSeeds(seeds, mb42, eSetup);
0161
0162
0163
0164
0165 CosmicMuonSeedGenerator::MuonRecHitPairVector mb31 = makeSegPairs(RHMBs[1], RHMBs[3], "mb31");
0166 createSeeds(seeds, mb31, eSetup);
0167
0168
0169
0170
0171 if (!allHits.empty()) {
0172 MuonRecHitContainer goodhits = selectSegments(allHits);
0173 LogTrace(category) << "good RecHits: " << goodhits.size();
0174
0175 if (goodhits.empty()) {
0176 LogTrace(category) << "No qualified Segments in Event! ";
0177 LogTrace(category) << "Use 2D RecHit";
0178
0179 createSeeds(seeds, allHits, eSetup);
0180
0181 } else {
0182 createSeeds(seeds, goodhits, eSetup);
0183 }
0184 }
0185
0186 LogTrace(category) << "Seeds built: " << seeds.size();
0187
0188 for (std::vector<TrajectorySeed>::iterator seed = seeds.begin(); seed != seeds.end(); ++seed) {
0189 output->push_back(*seed);
0190 }
0191
0192 event.put(std::move(output));
0193 seeds.clear();
0194 }
0195
0196 bool CosmicMuonSeedGenerator::checkQuality(const MuonRecHitPointer& hit) const {
0197 const std::string category = "Muon|RecoMuon|CosmicMuonSeedGenerator";
0198
0199
0200 if (!hit->isValid())
0201 return false;
0202
0203 if (hit->dimension() < 4) {
0204 LogTrace(category) << "dim < 4";
0205 return false;
0206 }
0207
0208 if (hit->isDT() && (hit->chi2() > theMaxDTChi2)) {
0209 LogTrace(category) << "DT chi2 too large";
0210 return false;
0211 } else if (hit->isCSC() && (hit->chi2() > theMaxCSCChi2)) {
0212 LogTrace(category) << "CSC chi2 too large";
0213 return false;
0214 }
0215 return true;
0216 }
0217
0218 MuonRecHitContainer CosmicMuonSeedGenerator::selectSegments(const MuonRecHitContainer& hits) const {
0219 MuonRecHitContainer result;
0220 const std::string category = "Muon|RecoMuon|CosmicMuonSeedGenerator";
0221
0222
0223 for (MuonRecHitContainer::const_iterator hit = hits.begin(); hit != hits.end(); hit++) {
0224 if (checkQuality(*hit))
0225 result.push_back(*hit);
0226 }
0227
0228 if (result.size() < 2)
0229 return result;
0230
0231 MuonRecHitContainer result2;
0232
0233
0234 for (MuonRecHitContainer::iterator hit = result.begin(); hit != result.end(); hit++) {
0235 if (*hit == nullptr)
0236 continue;
0237 if (!(*hit)->isValid())
0238 continue;
0239 bool good = true;
0240
0241
0242 for (MuonRecHitContainer::iterator hit2 = hit + 1; hit2 != result.end(); hit2++) {
0243 if (*hit2 == nullptr)
0244 continue;
0245 if (!(*hit2)->isValid())
0246 continue;
0247
0248
0249
0250
0251 if (!areCorrelated((*hit), (*hit2)))
0252 continue;
0253
0254 if (!leftIsBetter((*hit), (*hit2))) {
0255 good = false;
0256 } else
0257 (*hit2) = nullptr;
0258 }
0259
0260 if (good)
0261 result2.push_back(*hit);
0262 }
0263
0264 result.clear();
0265
0266 return result2;
0267 }
0268
0269 void CosmicMuonSeedGenerator::createSeeds(TrajectorySeedCollection& results,
0270 const MuonRecHitContainer& hits,
0271 const edm::EventSetup& eSetup) const {
0272 const std::string category = "Muon|RecoMuon|CosmicMuonSeedGenerator";
0273
0274 if (hits.empty() || results.size() >= theMaxSeeds)
0275 return;
0276 for (MuonRecHitContainer::const_iterator ihit = hits.begin(); ihit != hits.end(); ihit++) {
0277 const std::vector<TrajectorySeed>& sds = createSeed((*ihit), eSetup);
0278 LogTrace(category) << "created seeds from rechit " << sds.size();
0279 results.insert(results.end(), sds.begin(), sds.end());
0280 if (results.size() >= theMaxSeeds)
0281 break;
0282 }
0283 return;
0284 }
0285
0286 void CosmicMuonSeedGenerator::createSeeds(TrajectorySeedCollection& results,
0287 const CosmicMuonSeedGenerator::MuonRecHitPairVector& hitpairs,
0288 const edm::EventSetup& eSetup) const {
0289 const std::string category = "Muon|RecoMuon|CosmicMuonSeedGenerator";
0290
0291 if (hitpairs.empty() || results.size() >= theMaxSeeds)
0292 return;
0293 for (CosmicMuonSeedGenerator::MuonRecHitPairVector::const_iterator ihitpair = hitpairs.begin();
0294 ihitpair != hitpairs.end();
0295 ihitpair++) {
0296 const std::vector<TrajectorySeed>& sds = createSeed((*ihitpair), eSetup);
0297 LogTrace(category) << "created seeds from rechit " << sds.size();
0298 results.insert(results.end(), sds.begin(), sds.end());
0299 if (results.size() >= theMaxSeeds)
0300 break;
0301 }
0302 return;
0303 }
0304
0305 std::vector<TrajectorySeed> CosmicMuonSeedGenerator::createSeed(const MuonRecHitPointer& hit,
0306 const edm::EventSetup& eSetup) const {
0307 std::vector<TrajectorySeed> result;
0308
0309 const std::string category = "Muon|RecoMuon|CosmicMuonSeedGenerator";
0310
0311 MuonPatternRecoDumper dumper;
0312
0313
0314 double pt = 10.0;
0315
0316
0317 AlgebraicSymMatrix mat(5, 0);
0318
0319
0320 LocalPoint segPos = hit->localPosition();
0321
0322 GlobalVector polar(GlobalVector::Spherical(hit->globalDirection().theta(), hit->globalDirection().phi(), 1.));
0323
0324 if (theForcePointDownFlag) {
0325 if (hit->geographicalId().subdetId() == MuonSubdetId::DT && fabs(hit->globalDirection().eta()) < 4.0 &&
0326 hit->globalDirection().phi() > 0)
0327 polar = -polar;
0328
0329 if (hit->geographicalId().subdetId() == MuonSubdetId::CSC && fabs(hit->globalDirection().eta()) > 2.3)
0330 polar = -polar;
0331 }
0332
0333 polar *= fabs(pt) / polar.perp();
0334
0335 LocalVector segDir = hit->det()->toLocal(polar);
0336
0337 int charge = 1;
0338 LocalTrajectoryParameters param(segPos, segDir, charge);
0339
0340 charge = -1;
0341 LocalTrajectoryParameters param2(segPos, segDir, charge);
0342
0343 mat = hit->parametersError().similarityT(hit->projectionMatrix());
0344
0345 float p_err = 0.2;
0346 mat[0][0] = p_err;
0347
0348 LocalTrajectoryError error(asSMatrix<5>(mat));
0349
0350
0351 TrajectoryStateOnSurface tsos(param, error, hit->det()->surface(), &*theField);
0352 TrajectoryStateOnSurface tsos2(param2, error, hit->det()->surface(), &*theField);
0353
0354 LogTrace(category) << "Trajectory State on Surface of Seed";
0355 LogTrace(category) << "mom: " << tsos.globalMomentum() << " phi: " << tsos.globalMomentum().phi();
0356 LogTrace(category) << "pos: " << tsos.globalPosition();
0357 LogTrace(category) << "The RecSegment relies on: ";
0358 LogTrace(category) << dumper.dumpMuonId(hit->geographicalId());
0359
0360 edm::OwnVector<TrackingRecHit> container;
0361 container.push_back(hit->hit()->clone());
0362
0363 result.push_back(tsosToSeed(tsos, hit->geographicalId().rawId(), container));
0364 result.push_back(tsosToSeed(tsos2, hit->geographicalId().rawId(), container));
0365
0366 return result;
0367 }
0368
0369 bool CosmicMuonSeedGenerator::areCorrelated(const MuonRecHitPointer& lhs, const MuonRecHitPointer& rhs) const {
0370 bool result = false;
0371
0372 GlobalVector dir1 = lhs->globalDirection();
0373 GlobalPoint pos1 = lhs->globalPosition();
0374 GlobalVector dir2 = rhs->globalDirection();
0375 GlobalPoint pos2 = rhs->globalPosition();
0376
0377 GlobalVector dis = pos2 - pos1;
0378
0379 if ((deltaR<double>(dir1.eta(), dir1.phi(), dir2.eta(), dir2.phi()) < 0.1 ||
0380 deltaR<double>(dir1.eta(), dir1.phi(), -dir2.eta(), -dir2.phi()) < 0.1) &&
0381 dis.mag() < 5.0)
0382 result = true;
0383
0384 if ((deltaR<double>(dir1.eta(), dir1.phi(), dir2.eta(), dir2.phi()) < 0.1 ||
0385 deltaR<double>(dir1.eta(), dir1.phi(), -dir2.eta(), -dir2.phi()) < 0.1) &&
0386 (deltaR<double>(dir1.eta(), dir1.phi(), dis.eta(), dis.phi()) < 0.1 ||
0387 deltaR<double>(dir2.eta(), dir2.phi(), dis.eta(), dis.phi()) < 0.1))
0388 result = true;
0389
0390 if (fabs(dir1.eta()) > 4.0 || fabs(dir2.eta()) > 4.0) {
0391 if ((fabs(dir1.theta() - dir2.theta()) < 0.07 || fabs(dir1.theta() + dir2.theta()) > 3.07) &&
0392 (fabs(dir1.theta() - dis.theta()) < 0.07 || fabs(dir1.theta() - dis.theta()) < 0.07 ||
0393 fabs(dir1.theta() + dis.theta()) > 3.07 || fabs(dir1.theta() + dis.theta()) > 3.07))
0394
0395 result = true;
0396 }
0397
0398 return result;
0399 }
0400
0401 bool CosmicMuonSeedGenerator::leftIsBetter(const MuonTransientTrackingRecHit::MuonRecHitPointer& lhs,
0402 const MuonTransientTrackingRecHit::MuonRecHitPointer& rhs) const {
0403 if ((lhs->degreesOfFreedom() > rhs->degreesOfFreedom()) ||
0404 ((lhs->degreesOfFreedom() == rhs->degreesOfFreedom()) && (lhs)->chi2() < (rhs)->chi2()))
0405 return true;
0406 else
0407 return false;
0408 }
0409
0410 CosmicMuonSeedGenerator::MuonRecHitPairVector CosmicMuonSeedGenerator::makeSegPairs(
0411 const MuonTransientTrackingRecHit::MuonRecHitContainer& hits1,
0412 const MuonTransientTrackingRecHit::MuonRecHitContainer& hits2,
0413 std::string tag) const {
0414 MuonRecHitPairVector result;
0415 const std::string category = "Muon|RecoMuon|CosmicMuonSeedGenerator";
0416
0417 if (hits1.empty() || hits2.empty())
0418 return result;
0419
0420 for (MuonRecHitContainer::const_iterator ihit1 = hits1.begin(); ihit1 != hits1.end(); ihit1++) {
0421 if (!checkQuality(*ihit1))
0422 continue;
0423
0424 for (MuonRecHitContainer::const_iterator ihit2 = hits2.begin(); ihit2 != hits2.end(); ihit2++) {
0425 if (!checkQuality(*ihit2))
0426 continue;
0427
0428 float dphi = deltaPhi((*ihit1)->globalPosition().barePhi(), (*ihit2)->globalPosition().barePhi());
0429 if (dphi < 0.5) {
0430 if ((*ihit1)->globalPosition().y() > 0.0 && ((*ihit1)->globalPosition().y() > (*ihit2)->globalPosition().y())) {
0431 std::string tag2 = "top" + tag;
0432
0433 result.push_back(MuonRecHitPair(*ihit1, *ihit2, tag2));
0434 } else if ((*ihit1)->globalPosition().y() < 0.0 &&
0435 ((*ihit1)->globalPosition().y() < (*ihit2)->globalPosition().y())) {
0436 std::string tag2 = "bottom" + tag;
0437 result.push_back(MuonRecHitPair(*ihit2, *ihit1, tag2));
0438 }
0439 }
0440 }
0441 }
0442
0443 return result;
0444 }
0445
0446 std::vector<TrajectorySeed> CosmicMuonSeedGenerator::createSeed(const CosmicMuonSeedGenerator::MuonRecHitPair& hitpair,
0447 const edm::EventSetup& eSetup) const {
0448 std::vector<TrajectorySeed> result;
0449
0450 const std::string category = "Muon|RecoMuon|CosmicMuonSeedGenerator";
0451
0452 MuonPatternRecoDumper dumper;
0453
0454 float dphi = deltaPhi((hitpair.first)->globalDirection().barePhi(), (hitpair.second)->globalDirection().barePhi());
0455
0456 LogTrace(category) << "hitpair.type " << hitpair.type;
0457
0458 map<string, float>::const_iterator iterPar = theParameters.find(hitpair.type);
0459 if (iterPar == theParameters.end()) {
0460 return result;
0461 }
0462
0463
0464 int charge = (dphi > 0) ? -1 : 1;
0465
0466 double pt = 999.0;
0467 float paraC = (iterPar->second);
0468
0469 if (fabs(dphi) > 1e-5) {
0470 pt = paraC / fabs(dphi);
0471 }
0472
0473 if (pt < 10.0) {
0474 return result;
0475 }
0476
0477 AlgebraicVector t(4);
0478 AlgebraicSymMatrix mat(5, 0);
0479
0480 MuonTransientTrackingRecHit::MuonRecHitPointer hit = hitpair.first;
0481 if (hit->dimension() < (hitpair.second)->dimension())
0482 hit = hitpair.second;
0483
0484
0485 LocalPoint segPos = hit->localPosition();
0486
0487 GlobalVector polar(GlobalVector::Spherical(hit->globalDirection().theta(), hit->globalDirection().phi(), 1.));
0488
0489 if (theForcePointDownFlag) {
0490 if (hit->geographicalId().subdetId() == MuonSubdetId::DT && fabs(hit->globalDirection().eta()) < 4.0 &&
0491 hit->globalDirection().phi() > 0)
0492 polar = -polar;
0493
0494 if (hit->geographicalId().subdetId() == MuonSubdetId::CSC && fabs(hit->globalDirection().eta()) > 2.3)
0495 polar = -polar;
0496 }
0497
0498 polar *= fabs(pt) / polar.perp();
0499
0500 LocalVector segDir = hit->det()->toLocal(polar);
0501
0502 LocalTrajectoryParameters param(segPos, segDir, charge);
0503
0504 mat = hit->parametersError().similarityT(hit->projectionMatrix());
0505
0506 float p_err = 0.004 / paraC;
0507 if (pt < 10.01)
0508 p_err = 0.1;
0509 mat[0][0] = p_err;
0510
0511 LocalTrajectoryError error(asSMatrix<5>(mat));
0512
0513
0514 TrajectoryStateOnSurface tsos(param, error, hit->det()->surface(), &*theField);
0515
0516 LogTrace(category) << "Trajectory State on Surface of Seed";
0517 LogTrace(category) << "mom: " << tsos.globalMomentum() << " phi: " << tsos.globalMomentum().phi();
0518 LogTrace(category) << "pos: " << tsos.globalPosition();
0519 LogTrace(category) << "The RecSegment relies on: ";
0520 LogTrace(category) << dumper.dumpMuonId(hit->geographicalId());
0521
0522 edm::OwnVector<TrackingRecHit> container;
0523 container.push_back(hitpair.first->hit()->clone());
0524 container.push_back(hitpair.second->hit()->clone());
0525
0526 result.push_back(tsosToSeed(tsos, hit->geographicalId().rawId(), container));
0527
0528 return result;
0529 }
0530
0531 TrajectorySeed CosmicMuonSeedGenerator::tsosToSeed(const TrajectoryStateOnSurface& tsos, uint32_t id) const {
0532 edm::OwnVector<TrackingRecHit> container;
0533 return tsosToSeed(tsos, id, container);
0534 }
0535
0536 TrajectorySeed CosmicMuonSeedGenerator::tsosToSeed(const TrajectoryStateOnSurface& tsos,
0537 uint32_t id,
0538 edm::OwnVector<TrackingRecHit>& container) const {
0539 PTrajectoryStateOnDet const& seedTSOS = trajectoryStateTransform::persistentState(tsos, id);
0540 TrajectorySeed seed(seedTSOS, container, alongMomentum);
0541 return seed;
0542 }