File indexing completed on 2025-01-09 23:33:56
0001
0002
0003
0004
0005
0006
0007 #include <RecoMuon/MuonSeedGenerator/interface/MuonSeedCleaner.h>
0008
0009
0010 #include <DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h>
0011 #include <DataFormats/TrajectorySeed/interface/TrajectorySeed.h>
0012 #include <DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h>
0013 #include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
0014 #include <DataFormats/CSCRecHit/interface/CSCSegmentCollection.h>
0015 #include <DataFormats/CSCRecHit/interface/CSCSegment.h>
0016 #include <DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h>
0017 #include <DataFormats/DTRecHit/interface/DTRecSegment4D.h>
0018
0019
0020 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0021 #include <TrackingTools/DetLayers/interface/DetLayer.h>
0022 #include <RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h>
0023 #include <RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h>
0024 #include <RecoMuon/Records/interface/MuonRecoGeometryRecord.h>
0025
0026
0027 #include <TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h>
0028 #include <DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h>
0029
0030
0031 #include <FWCore/ParameterSet/interface/ParameterSet.h>
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include <FWCore/Framework/interface/ESHandle.h>
0034 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0035 #include <DataFormats/Common/interface/Handle.h>
0036
0037
0038 #include <vector>
0039 #include <deque>
0040 #include <utility>
0041
0042
0043
0044 static bool lengthSorting(const TrajectorySeed& s1, const TrajectorySeed& s2) { return (s1.nHits() > s2.nHits()); }
0045
0046
0047
0048
0049 MuonSeedCleaner::MuonSeedCleaner(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC) {
0050
0051 debug = pset.getParameter<bool>("DebugMuonSeed");
0052
0053
0054 edm::ParameterSet serviceParameters = pset.getParameter<edm::ParameterSet>("ServiceParameters");
0055 theService = new MuonServiceProxy(serviceParameters, std::move(iC));
0056 }
0057
0058
0059
0060
0061 MuonSeedCleaner::~MuonSeedCleaner() {
0062 if (theService)
0063 delete theService;
0064 }
0065
0066
0067
0068
0069
0070
0071
0072 std::vector<TrajectorySeed> MuonSeedCleaner::seedCleaner(const edm::EventSetup& eventSetup,
0073 std::vector<TrajectorySeed>& seeds) {
0074 theService->update(eventSetup);
0075
0076 std::vector<TrajectorySeed> FinalSeeds;
0077
0078
0079 std::vector<SeedContainer> theCollection = GroupSeeds(seeds);
0080
0081
0082 for (size_t i = 0; i < theCollection.size(); i++) {
0083
0084 SeedContainer goodSeeds = SeedCandidates(theCollection[i], true);
0085 SeedContainer otherSeeds = SeedCandidates(theCollection[i], false);
0086 if (MomentumFilter(goodSeeds)) {
0087
0088 TrajectorySeed bestSeed = Chi2LengthSelection(goodSeeds);
0089 FinalSeeds.push_back(bestSeed);
0090
0091 GlobalPoint seedgp = SeedPosition(bestSeed);
0092 double eta = fabs(seedgp.eta());
0093 if (goodSeeds.size() > 2 && eta > 1.5) {
0094 TrajectorySeed anotherSeed = MoreRecHits(goodSeeds);
0095 FinalSeeds.push_back(anotherSeed);
0096 }
0097 } else if (MomentumFilter(otherSeeds)) {
0098
0099 TrajectorySeed bestSeed = MoreRecHits(otherSeeds);
0100 FinalSeeds.push_back(bestSeed);
0101
0102 GlobalPoint seedgp = SeedPosition(bestSeed);
0103 double eta = fabs(seedgp.eta());
0104 if (otherSeeds.size() > 2 && eta > 1.5) {
0105 TrajectorySeed anotherSeed = LeanHighMomentum(otherSeeds);
0106 FinalSeeds.push_back(anotherSeed);
0107 }
0108 } else {
0109
0110 TrajectorySeed bestSeed = LeanHighMomentum(theCollection[i]);
0111 FinalSeeds.push_back(bestSeed);
0112
0113 GlobalPoint seedgp = SeedPosition(bestSeed);
0114 double eta = fabs(seedgp.eta());
0115 if (theCollection.size() > 2 && eta > 1.5) {
0116 TrajectorySeed anotherSeed = BiggerCone(theCollection[i]);
0117 FinalSeeds.push_back(anotherSeed);
0118 }
0119 }
0120 }
0121 return FinalSeeds;
0122 }
0123
0124 TrajectorySeed MuonSeedCleaner::Chi2LengthSelection(std::vector<TrajectorySeed>& seeds) {
0125 if (seeds.size() == 1)
0126 return seeds[0];
0127
0128 int winner = 0;
0129 int moreHits = 0;
0130 double bestChi2 = 99999.;
0131 for (size_t i = 0; i < seeds.size(); i++) {
0132
0133
0134
0135
0136
0137 double theChi2 = SeedChi2(seeds[i]);
0138 double dChi2 = fabs(1. - (theChi2 / bestChi2));
0139 int theHits = seeds[i].nHits();
0140 int dHits = theHits - moreHits;
0141
0142
0143
0144 if (theChi2 < bestChi2 && dChi2 > 0.05) {
0145 winner = static_cast<int>(i);
0146 bestChi2 = theChi2;
0147 moreHits = theHits;
0148 }
0149
0150 if (theChi2 >= bestChi2 && dChi2 < 0.05 && dHits > 0) {
0151 winner = static_cast<int>(i);
0152 bestChi2 = theChi2;
0153 moreHits = theHits;
0154 }
0155 }
0156
0157 TrajectorySeed theSeed = seeds[winner];
0158 seeds.erase(seeds.begin() + winner);
0159 return theSeed;
0160 }
0161
0162 TrajectorySeed MuonSeedCleaner::BiggerCone(std::vector<TrajectorySeed>& seeds) {
0163 if (seeds.size() == 1)
0164 return seeds[0];
0165
0166 float biggerProjErr = 9999.;
0167 int winner = 0;
0168 AlgebraicSymMatrix mat(5, 0);
0169 for (size_t i = 0; i < seeds.size(); i++) {
0170 auto r1 = seeds[i].recHits().begin();
0171 mat = r1->parametersError().similarityT(r1->projectionMatrix());
0172
0173 int NRecHits = NRecHitsFromSegment(*r1);
0174
0175 float ddx = mat[1][1];
0176 float ddy = mat[2][2];
0177 float dxx = mat[3][3];
0178 float dyy = mat[4][4];
0179 float projectErr = sqrt((ddx * 10000.) + (ddy * 10000.) + dxx + dyy);
0180
0181 if (NRecHits < 5)
0182 continue;
0183 if (projectErr < biggerProjErr)
0184 continue;
0185
0186 winner = static_cast<int>(i);
0187 biggerProjErr = projectErr;
0188 }
0189 TrajectorySeed theSeed = seeds[winner];
0190 seeds.erase(seeds.begin() + winner);
0191 return theSeed;
0192 }
0193
0194 TrajectorySeed MuonSeedCleaner::LeanHighMomentum(std::vector<TrajectorySeed>& seeds) {
0195 if (seeds.size() == 1)
0196 return seeds[0];
0197
0198 double highestPt = 0.;
0199 int winner = 0;
0200 for (size_t i = 0; i < seeds.size(); i++) {
0201 GlobalVector mom = SeedMomentum(seeds[i]);
0202 double pt = sqrt((mom.x() * mom.x()) + (mom.y() * mom.y()));
0203 if (pt > highestPt) {
0204 winner = static_cast<int>(i);
0205 highestPt = pt;
0206 }
0207 }
0208 TrajectorySeed theSeed = seeds[winner];
0209 seeds.erase(seeds.begin() + winner);
0210 return theSeed;
0211 }
0212
0213 TrajectorySeed MuonSeedCleaner::MoreRecHits(std::vector<TrajectorySeed>& seeds) {
0214 if (seeds.size() == 1)
0215 return seeds[0];
0216
0217 int winner = 0;
0218 int moreHits = 0;
0219 double betterChi2 = 99999.;
0220 for (size_t i = 0; i < seeds.size(); i++) {
0221 int theHits = 0;
0222 for (auto const& r1 : seeds[i].recHits()) {
0223 theHits += NRecHitsFromSegment(r1);
0224 }
0225
0226 double theChi2 = SeedChi2(seeds[i]);
0227
0228 if (theHits == moreHits && theChi2 < betterChi2) {
0229 betterChi2 = theChi2;
0230 winner = static_cast<int>(i);
0231 }
0232 if (theHits > moreHits) {
0233 moreHits = theHits;
0234 betterChi2 = theChi2;
0235 winner = static_cast<int>(i);
0236 }
0237 }
0238 TrajectorySeed theSeed = seeds[winner];
0239 seeds.erase(seeds.begin() + winner);
0240 return theSeed;
0241 }
0242
0243 SeedContainer MuonSeedCleaner::LengthFilter(std::vector<TrajectorySeed>& seeds) {
0244 SeedContainer longSeeds;
0245 int NSegs = 0;
0246 for (size_t i = 0; i < seeds.size(); i++) {
0247 int theLength = static_cast<int>(seeds[i].nHits());
0248 if (theLength > NSegs) {
0249 NSegs = theLength;
0250 longSeeds.clear();
0251 longSeeds.push_back(seeds[i]);
0252 } else if (theLength == NSegs) {
0253 longSeeds.push_back(seeds[i]);
0254 } else {
0255 continue;
0256 }
0257 }
0258
0259
0260 return longSeeds;
0261 }
0262
0263 bool MuonSeedCleaner::MomentumFilter(std::vector<TrajectorySeed>& seeds) {
0264 bool findgoodMomentum = false;
0265 SeedContainer goodMomentumSeeds = seeds;
0266 seeds.clear();
0267 for (size_t i = 0; i < goodMomentumSeeds.size(); i++) {
0268 GlobalVector mom = SeedMomentum(goodMomentumSeeds[i]);
0269 double pt = sqrt((mom.x() * mom.x()) + (mom.y() * mom.y()));
0270 if (pt < 6. || pt > 2000.)
0271 continue;
0272
0273
0274 seeds.push_back(goodMomentumSeeds[i]);
0275 findgoodMomentum = true;
0276 }
0277 if (seeds.empty())
0278 seeds = goodMomentumSeeds;
0279
0280 return findgoodMomentum;
0281 }
0282
0283 SeedContainer MuonSeedCleaner::SeedCandidates(std::vector<TrajectorySeed>& seeds, bool good) {
0284 SeedContainer theCandidate;
0285 theCandidate.clear();
0286
0287 bool longSeed = false;
0288 bool withFirstLayer = false;
0289
0290
0291 for (size_t i = 0; i < seeds.size(); i++) {
0292 if (seeds[i].nHits() > 1)
0293 longSeed = true;
0294
0295
0296 for (auto const& r1 : seeds[i].recHits()) {
0297 const GeomDet* gdet = theService->trackingGeometry()->idToDet(r1.geographicalId());
0298 DetId geoId = gdet->geographicalId();
0299
0300 if (geoId.subdetId() == MuonSubdetId::DT) {
0301 DTChamberId DT_Id(r1.geographicalId());
0302
0303 if (DT_Id.station() != 1)
0304 continue;
0305 withFirstLayer = true;
0306 }
0307 if (geoId.subdetId() == MuonSubdetId::CSC) {
0308 CSCDetId CSC_Id = CSCDetId(r1.geographicalId());
0309
0310 if (CSC_Id.station() != 1)
0311 continue;
0312 withFirstLayer = true;
0313 }
0314 }
0315 bool goodseed = (longSeed && withFirstLayer) ? true : false;
0316
0317 if (goodseed == good)
0318 theCandidate.push_back(seeds[i]);
0319 }
0320 return theCandidate;
0321 }
0322
0323 std::vector<SeedContainer> MuonSeedCleaner::GroupSeeds(std::vector<TrajectorySeed>& seeds) {
0324 std::vector<SeedContainer> seedCollection;
0325 seedCollection.clear();
0326 std::vector<TrajectorySeed> theGroup;
0327 std::vector<bool> usedSeed(seeds.size(), false);
0328
0329
0330 for (size_t i = 0; i < seeds.size(); i++) {
0331 if (usedSeed[i])
0332 continue;
0333 theGroup.push_back(seeds[i]);
0334 usedSeed[i] = true;
0335
0336 GlobalPoint pos1 = SeedPosition(seeds[i]);
0337
0338 for (size_t j = i + 1; j < seeds.size(); j++) {
0339
0340 unsigned int overlapping = OverlapSegments(seeds[i], seeds[j]);
0341 if (!usedSeed[j] && overlapping > 0) {
0342
0343 if (seeds[i].nHits() == overlapping && seeds[j].nHits() == overlapping) {
0344 usedSeed[j] = true;
0345 continue;
0346 }
0347 theGroup.push_back(seeds[j]);
0348 usedSeed[j] = true;
0349 }
0350 if (usedSeed[j])
0351 continue;
0352
0353
0354 GlobalPoint pos2 = SeedPosition(seeds[j]);
0355 double dh = pos1.eta() - pos2.eta();
0356 double df = pos1.phi() - pos2.phi();
0357 double dR = sqrt((dh * dh) + (df * df));
0358
0359 if (dR > 0.3 && seeds[j].nHits() == 1)
0360 continue;
0361 if (dR > 0.2 && seeds[j].nHits() > 1)
0362 continue;
0363 theGroup.push_back(seeds[j]);
0364 usedSeed[j] = true;
0365 }
0366 sort(theGroup.begin(), theGroup.end(), lengthSorting);
0367 seedCollection.push_back(theGroup);
0368
0369 theGroup.clear();
0370 }
0371 return seedCollection;
0372 }
0373
0374 unsigned int MuonSeedCleaner::OverlapSegments(const TrajectorySeed& seed1, const TrajectorySeed& seed2) {
0375 unsigned int overlapping = 0;
0376 for (auto const& r1 : seed1.recHits()) {
0377 DetId id1 = r1.geographicalId();
0378 const GeomDet* gdet1 = theService->trackingGeometry()->idToDet(id1);
0379 GlobalPoint gp1 = gdet1->toGlobal(r1.localPosition());
0380
0381 for (auto const& r2 : seed2.recHits()) {
0382 DetId id2 = r2.geographicalId();
0383 if (id1 != id2)
0384 continue;
0385
0386 const GeomDet* gdet2 = theService->trackingGeometry()->idToDet(id2);
0387 GlobalPoint gp2 = gdet2->toGlobal(r2.localPosition());
0388
0389 double dx = gp1.x() - gp2.x();
0390 double dy = gp1.y() - gp2.y();
0391 double dz = gp1.z() - gp2.z();
0392 double dL = sqrt(dx * dx + dy * dy + dz * dz);
0393
0394 if (dL < 1.)
0395 overlapping++;
0396 }
0397 }
0398 return overlapping;
0399 }
0400
0401 double MuonSeedCleaner::SeedChi2(const TrajectorySeed& seed) {
0402 double theChi2 = 0.;
0403 for (auto const& r1 : seed.recHits()) {
0404
0405 theChi2 += NChi2OfSegment(r1);
0406 }
0407 theChi2 = theChi2 / seed.nHits();
0408
0409
0410 return theChi2;
0411 }
0412
0413 int MuonSeedCleaner::SeedLength(const TrajectorySeed& seed) {
0414 int theHits = 0;
0415 for (auto const& recHit : seed.recHits()) {
0416
0417 theHits += NRecHitsFromSegment(recHit);
0418 }
0419
0420
0421 return theHits;
0422 }
0423
0424 GlobalPoint MuonSeedCleaner::SeedPosition(const TrajectorySeed& seed) {
0425 PTrajectoryStateOnDet pTSOD = seed.startingState();
0426 DetId SeedDetId(pTSOD.detId());
0427 const GeomDet* geoDet = theService->trackingGeometry()->idToDet(SeedDetId);
0428 TrajectoryStateOnSurface SeedTSOS =
0429 trajectoryStateTransform::transientState(pTSOD, &(geoDet->surface()), &*theService->magneticField());
0430 GlobalPoint pos = SeedTSOS.globalPosition();
0431
0432 return pos;
0433 }
0434
0435 GlobalVector MuonSeedCleaner::SeedMomentum(const TrajectorySeed& seed) {
0436 PTrajectoryStateOnDet pTSOD = seed.startingState();
0437 DetId SeedDetId(pTSOD.detId());
0438 const GeomDet* geoDet = theService->trackingGeometry()->idToDet(SeedDetId);
0439 TrajectoryStateOnSurface SeedTSOS =
0440 trajectoryStateTransform::transientState(pTSOD, &(geoDet->surface()), &*theService->magneticField());
0441 GlobalVector mom = SeedTSOS.globalMomentum();
0442
0443 return mom;
0444 }
0445
0446 int MuonSeedCleaner::NRecHitsFromSegment(const TrackingRecHit& rhit) {
0447 int NRechits = 0;
0448 const GeomDet* gdet = theService->trackingGeometry()->idToDet(rhit.geographicalId());
0449 MuonTransientTrackingRecHit::MuonRecHitPointer theSeg =
0450 MuonTransientTrackingRecHit::specificBuild(gdet, rhit.clone());
0451
0452 DetId geoId = gdet->geographicalId();
0453 if (geoId.subdetId() == MuonSubdetId::DT) {
0454 DTChamberId DT_Id(rhit.geographicalId());
0455 std::vector<TrackingRecHit*> DThits = theSeg->recHits();
0456 int dt1DHits = 0;
0457 for (size_t j = 0; j < DThits.size(); j++) {
0458 dt1DHits += (DThits[j]->recHits()).size();
0459 }
0460 NRechits = dt1DHits;
0461 }
0462
0463 if (geoId.subdetId() == MuonSubdetId::CSC) {
0464 NRechits = (theSeg->recHits()).size();
0465 }
0466 return NRechits;
0467 }
0468
0469 int MuonSeedCleaner::NRecHitsFromSegment(MuonTransientTrackingRecHit* rhit) {
0470 int NRechits = 0;
0471 DetId geoId = rhit->geographicalId();
0472 if (geoId.subdetId() == MuonSubdetId::DT) {
0473 DTChamberId DT_Id(geoId);
0474 std::vector<TrackingRecHit*> DThits = rhit->recHits();
0475 int dt1DHits = 0;
0476 for (size_t j = 0; j < DThits.size(); j++) {
0477 dt1DHits += (DThits[j]->recHits()).size();
0478 }
0479 NRechits = dt1DHits;
0480
0481 }
0482 if (geoId.subdetId() == MuonSubdetId::CSC) {
0483 NRechits = (rhit->recHits()).size();
0484
0485 }
0486 return NRechits;
0487 }
0488
0489 double MuonSeedCleaner::NChi2OfSegment(const TrackingRecHit& rhit) {
0490 double NChi2 = 999999.;
0491 const GeomDet* gdet = theService->trackingGeometry()->idToDet(rhit.geographicalId());
0492 MuonTransientTrackingRecHit::MuonRecHitPointer theSeg =
0493 MuonTransientTrackingRecHit::specificBuild(gdet, rhit.clone());
0494
0495 double dof = static_cast<double>(theSeg->degreesOfFreedom());
0496 NChi2 = theSeg->chi2() / dof;
0497
0498
0499 return NChi2;
0500 }