File indexing completed on 2025-01-09 23:33:59
0001
0002
0003
0004
0005 #include "RecoMuon/MuonSeedGenerator/interface/SETPatternRecognition.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0011 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0012 #include "TMath.h"
0013
0014 using namespace edm;
0015 using namespace std;
0016
0017 SETPatternRecognition::SETPatternRecognition(const ParameterSet& parameterSet, edm::ConsumesCollector& iC)
0018 : MuonSeedVPatternRecognition(parameterSet.getParameter<ParameterSet>("SETTrajBuilderParameters")
0019 .getParameter<ParameterSet>("FilterParameters")) {
0020 const string metname = "Muon|RecoMuon|SETPatternRecognition";
0021
0022 ParameterSet trajectoryBuilderParameters = parameterSet.getParameter<ParameterSet>("SETTrajBuilderParameters");
0023
0024 ParameterSet filterPSet = trajectoryBuilderParameters.getParameter<ParameterSet>("FilterParameters");
0025 maxActiveChambers = filterPSet.getParameter<int>("maxActiveChambers");
0026 useRPCs = filterPSet.getParameter<bool>("EnableRPCMeasurement");
0027 DTRecSegmentLabel = filterPSet.getParameter<edm::InputTag>("DTRecSegmentLabel");
0028 CSCRecSegmentLabel = filterPSet.getParameter<edm::InputTag>("CSCRecSegmentLabel");
0029 RPCRecSegmentLabel = filterPSet.getParameter<edm::InputTag>("RPCRecSegmentLabel");
0030
0031 outsideChamberErrorScale = filterPSet.getParameter<double>("OutsideChamberErrorScale");
0032 minLocalSegmentAngle = filterPSet.getParameter<double>("MinLocalSegmentAngle");
0033
0034 dtToken = iC.consumes<DTRecSegment4DCollection>(DTRecSegmentLabel);
0035 cscToken = iC.consumes<CSCSegmentCollection>(CSCRecSegmentLabel);
0036 rpcToken = iC.consumes<RPCRecHitCollection>(RPCRecSegmentLabel);
0037 }
0038
0039
0040 void SETPatternRecognition::produce(const edm::Event& event,
0041 const edm::EventSetup& eventSetup,
0042 std::vector<MuonRecHitContainer>& segments_clusters) {
0043 const string metname = "Muon|RecoMuon|SETMuonSeedSeed";
0044
0045
0046 MuonRecHitContainer muonRecHits;
0047 MuonRecHitContainer muonRecHits_DT2D_hasPhi;
0048 MuonRecHitContainer muonRecHits_DT2D_hasZed;
0049 MuonRecHitContainer muonRecHits_RPC;
0050
0051
0052
0053
0054
0055 edm::Handle<DTRecSegment4DCollection> dtRecHits;
0056 event.getByToken(dtToken, dtRecHits);
0057 std::vector<DTChamberId> chambers_DT;
0058 std::vector<DTChamberId>::const_iterator chIt_DT;
0059 for (DTRecSegment4DCollection::const_iterator rechit = dtRecHits->begin(); rechit != dtRecHits->end(); ++rechit) {
0060 bool insert = true;
0061 for (chIt_DT = chambers_DT.begin(); chIt_DT != chambers_DT.end(); ++chIt_DT) {
0062 if (((*rechit).chamberId().wheel()) == ((*chIt_DT).wheel()) &&
0063 ((*rechit).chamberId().station() == (*chIt_DT).station()) &&
0064 ((*rechit).chamberId().sector() == (*chIt_DT).sector())) {
0065 insert = false;
0066 }
0067 }
0068 if (insert) {
0069 chambers_DT.push_back((*rechit).chamberId());
0070 }
0071 if (segmentCleaning((*rechit).geographicalId(),
0072 rechit->localPosition(),
0073 rechit->localPositionError(),
0074 rechit->localDirection(),
0075 rechit->localDirectionError(),
0076 rechit->chi2(),
0077 rechit->degreesOfFreedom())) {
0078 continue;
0079 }
0080 if ((rechit->hasZed() && rechit->hasPhi())) {
0081 muonRecHits.push_back(MuonTransientTrackingRecHit::specificBuild(
0082 theService->trackingGeometry()->idToDet((*rechit).geographicalId()), &*rechit));
0083 } else if (rechit->hasZed()) {
0084 muonRecHits_DT2D_hasZed.push_back(MuonTransientTrackingRecHit::specificBuild(
0085 theService->trackingGeometry()->idToDet((*rechit).geographicalId()), &*rechit));
0086 } else if (rechit->hasPhi()) {
0087 muonRecHits_DT2D_hasPhi.push_back(MuonTransientTrackingRecHit::specificBuild(
0088 theService->trackingGeometry()->idToDet((*rechit).geographicalId()), &*rechit));
0089 } else {
0090
0091 }
0092 }
0093
0094
0095
0096
0097
0098
0099 edm::Handle<CSCSegmentCollection> cscSegments;
0100 event.getByToken(cscToken, cscSegments);
0101 std::vector<CSCDetId> chambers_CSC;
0102 std::vector<CSCDetId>::const_iterator chIt_CSC;
0103 for (CSCSegmentCollection::const_iterator rechit = cscSegments->begin(); rechit != cscSegments->end(); ++rechit) {
0104 bool insert = true;
0105 for (chIt_CSC = chambers_CSC.begin(); chIt_CSC != chambers_CSC.end(); ++chIt_CSC) {
0106 if (((*rechit).cscDetId().chamber() == (*chIt_CSC).chamber()) &&
0107 ((*rechit).cscDetId().station() == (*chIt_CSC).station()) &&
0108 ((*rechit).cscDetId().ring() == (*chIt_CSC).ring()) &&
0109 ((*rechit).cscDetId().endcap() == (*chIt_CSC).endcap())) {
0110 insert = false;
0111 }
0112 }
0113 if (insert) {
0114 chambers_CSC.push_back((*rechit).cscDetId().chamberId());
0115 }
0116 if (segmentCleaning((*rechit).geographicalId(),
0117 rechit->localPosition(),
0118 rechit->localPositionError(),
0119 rechit->localDirection(),
0120 rechit->localDirectionError(),
0121 rechit->chi2(),
0122 rechit->degreesOfFreedom())) {
0123 continue;
0124 }
0125 muonRecHits.push_back(MuonTransientTrackingRecHit::specificBuild(
0126 theService->trackingGeometry()->idToDet((*rechit).geographicalId()), &*rechit));
0127 }
0128
0129
0130
0131
0132
0133
0134 edm::Handle<RPCRecHitCollection> rpcRecHits;
0135 event.getByToken(rpcToken, rpcRecHits);
0136 if (useRPCs) {
0137 for (RPCRecHitCollection::const_iterator rechit = rpcRecHits->begin(); rechit != rpcRecHits->end(); ++rechit) {
0138
0139 const LocalVector localDirection(0., 0., 1.);
0140 const LocalError localDirectionError(0., 0., 0.);
0141 const double chi2 = 1.;
0142 const int ndf = 1;
0143 if (segmentCleaning((*rechit).geographicalId(),
0144 rechit->localPosition(),
0145 rechit->localPositionError(),
0146 localDirection,
0147 localDirectionError,
0148 chi2,
0149 ndf)) {
0150 continue;
0151 }
0152 muonRecHits_RPC.push_back(MuonTransientTrackingRecHit::specificBuild(
0153 theService->trackingGeometry()->idToDet((*rechit).geographicalId()), &*rechit));
0154 }
0155 }
0156
0157
0158 if (int(chambers_DT.size() + chambers_CSC.size()) > maxActiveChambers) {
0159
0160
0161 edm::LogWarning("tooManyActiveChambers") << " Too many active chambers : nDT = " << chambers_DT.size()
0162 << " nCSC = " << chambers_CSC.size() << " Skip them all.";
0163 muonRecHits.clear();
0164 muonRecHits_DT2D_hasPhi.clear();
0165 muonRecHits_DT2D_hasZed.clear();
0166 muonRecHits_RPC.clear();
0167 }
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177 bool useDT2D_hasPhi = true;
0178 bool useDT2D_hasZed = true;
0179 double dXclusBoxMax = 0.60;
0180 double dYclusBoxMax = 0.;
0181
0182
0183
0184
0185
0186
0187
0188 dYclusBoxMax = 0.02;
0189
0190
0191
0192 float dXclus = 0.0;
0193 float dXclus_box = 0.0;
0194 float dYclus_box = 0.0;
0195
0196 MuonRecHitContainer temp;
0197
0198 std::vector<MuonRecHitContainer> seeds;
0199
0200 std::vector<float> running_meanX;
0201 std::vector<float> running_meanY;
0202
0203 std::vector<float> seed_minX;
0204 std::vector<float> seed_maxX;
0205 std::vector<float> seed_minY;
0206 std::vector<float> seed_maxY;
0207
0208
0209
0210
0211 for (MuonRecHitContainer::const_iterator it = muonRecHits.begin(); it != muonRecHits.end(); ++it) {
0212
0213
0214
0215
0216
0217
0218 temp.clear();
0219
0220 temp.push_back((*it));
0221
0222 seeds.push_back(temp);
0223
0224
0225
0226
0227 running_meanX.push_back((*it)->globalPosition().phi());
0228 running_meanY.push_back((*it)->globalPosition().theta());
0229
0230
0231 seed_minX.push_back((*it)->globalPosition().phi());
0232 seed_maxX.push_back((*it)->globalPosition().phi());
0233 seed_minY.push_back((*it)->globalPosition().theta());
0234 seed_maxY.push_back((*it)->globalPosition().theta());
0235 }
0236
0237
0238
0239 for (unsigned int NNN = 0; NNN < seeds.size(); ++NNN) {
0240 for (unsigned int MMM = NNN + 1; MMM < seeds.size(); ++MMM) {
0241 if (running_meanX[MMM] == 999999. || running_meanX[NNN] == 999999.) {
0242
0243
0244 continue;
0245 }
0246
0247
0248
0249 double temp_meanX = running_meanX[NNN];
0250 double temp_minX = seed_minX[NNN];
0251 double temp_maxX = seed_maxX[NNN];
0252
0253
0254
0255 dXclus = running_meanX[NNN] - running_meanX[MMM];
0256 if (dXclus > TMath::Pi()) {
0257 temp_meanX = temp_meanX - 2. * TMath::Pi();
0258 temp_minX = temp_minX - 2. * TMath::Pi();
0259 temp_maxX = temp_maxX - 2. * TMath::Pi();
0260 }
0261 if (dXclus < -TMath::Pi()) {
0262 temp_meanX = temp_meanX + 2. * TMath::Pi();
0263 temp_minX = temp_minX + 2. * TMath::Pi();
0264 temp_maxX = temp_maxX + 2. * TMath::Pi();
0265 }
0266
0267
0268
0269
0270
0271
0272 if (temp_meanX > running_meanX[MMM])
0273 dXclus_box = temp_minX - seed_maxX[MMM];
0274 else
0275 dXclus_box = seed_minX[MMM] - temp_maxX;
0276 if (running_meanY[NNN] > running_meanY[MMM])
0277 dYclus_box = seed_minY[NNN] - seed_maxY[MMM];
0278 else
0279 dYclus_box = seed_minY[MMM] - seed_maxY[NNN];
0280
0281 if (dXclus_box < dXclusBoxMax && dYclus_box < dYclusBoxMax) {
0282
0283
0284
0285
0286
0287 running_meanX[MMM] = (temp_meanX * seeds[NNN].size() + running_meanX[MMM] * seeds[MMM].size()) /
0288 (seeds[NNN].size() + seeds[MMM].size());
0289 running_meanY[MMM] = (running_meanY[NNN] * seeds[NNN].size() + running_meanY[MMM] * seeds[MMM].size()) /
0290 (seeds[NNN].size() + seeds[MMM].size());
0291
0292
0293
0294 if (temp_minX <= seed_minX[MMM])
0295 seed_minX[MMM] = temp_minX;
0296 if (temp_maxX > seed_maxX[MMM])
0297 seed_maxX[MMM] = temp_maxX;
0298 if (seed_minY[NNN] <= seed_minY[MMM])
0299 seed_minY[MMM] = seed_minY[NNN];
0300 if (seed_maxY[NNN] > seed_maxY[MMM])
0301 seed_maxY[MMM] = seed_maxY[NNN];
0302
0303
0304
0305 if (running_meanX[MMM] > TMath::Pi()) {
0306 running_meanX[MMM] = running_meanX[MMM] - 2. * TMath::Pi();
0307 seed_minX[MMM] = seed_minX[MMM] - 2. * TMath::Pi();
0308 seed_maxX[MMM] = seed_maxX[MMM] - 2. * TMath::Pi();
0309 }
0310 if (running_meanX[MMM] < -TMath::Pi()) {
0311 running_meanX[MMM] = running_meanX[MMM] + 2. * TMath::Pi();
0312 seed_minX[MMM] = seed_minX[MMM] + 2. * TMath::Pi();
0313 seed_maxX[MMM] = seed_maxX[MMM] + 2. * TMath::Pi();
0314 }
0315
0316
0317 seeds[MMM].insert(seeds[MMM].end(), seeds[NNN].begin(), seeds[NNN].end());
0318
0319
0320 running_meanX[NNN] = 999999.;
0321 running_meanY[NNN] = 999999.;
0322
0323
0324 break;
0325 }
0326 }
0327 }
0328 bool tooCloseClusters = false;
0329 if (seeds.size() > 1) {
0330 std::vector<double> seedTheta(seeds.size());
0331 for (unsigned int iSeed = 0; iSeed < seeds.size(); ++iSeed) {
0332 seedTheta[iSeed] = seeds[iSeed][0]->globalPosition().theta();
0333 if (iSeed) {
0334 double dTheta = fabs(seedTheta[iSeed] - seedTheta[iSeed - 1]);
0335 if (dTheta < 0.5) {
0336 tooCloseClusters = true;
0337 break;
0338 }
0339 }
0340 }
0341 }
0342
0343
0344
0345
0346
0347 for (unsigned int NNN = 0; NNN < seeds.size(); ++NNN) {
0348 if (running_meanX[NNN] == 999999.)
0349 continue;
0350
0351
0352 if (useDT2D_hasZed) {
0353 for (MuonRecHitContainer::const_iterator it2 = muonRecHits_DT2D_hasZed.begin();
0354 it2 != muonRecHits_DT2D_hasZed.end();
0355 ++it2) {
0356
0357 if (((*it2)->globalPosition().theta() < seed_maxY[NNN] + dYclusBoxMax) &&
0358 ((*it2)->globalPosition().theta() > seed_minY[NNN] - dYclusBoxMax)) {
0359
0360
0361
0362 if (!((((*it2)->globalPosition().phi() + 0.09) < (seed_minX[NNN] - dXclusBoxMax) &&
0363 ((*it2)->globalPosition().phi() - 0.09) < (seed_minX[NNN] - dXclusBoxMax)) ||
0364 (((*it2)->globalPosition().phi() + 0.09) > (seed_maxX[NNN] + dXclusBoxMax) &&
0365
0366 ((*it2)->globalPosition().phi() - 0.09) > (seed_maxX[NNN] + dXclusBoxMax)))) {
0367 seeds[NNN].push_back((*it2));
0368 }
0369 }
0370 }
0371 }
0372
0373
0374 if (useDT2D_hasPhi) {
0375 for (MuonRecHitContainer::const_iterator it2 = muonRecHits_DT2D_hasPhi.begin();
0376 it2 != muonRecHits_DT2D_hasPhi.end();
0377 ++it2) {
0378 if (((*it2)->globalPosition().phi() < seed_maxX[NNN] + dXclusBoxMax) &&
0379 ((*it2)->globalPosition().phi() > seed_minX[NNN] - dXclusBoxMax)) {
0380 if (!((((*it2)->globalPosition().theta() + 0.3) < (seed_minY[NNN] - dYclusBoxMax) &&
0381 ((*it2)->globalPosition().theta() - 0.3) < (seed_minY[NNN] - dYclusBoxMax)) ||
0382 (((*it2)->globalPosition().theta() + 0.3) > (seed_maxY[NNN] + dYclusBoxMax) &&
0383
0384 ((*it2)->globalPosition().theta() - 0.3) > (seed_maxY[NNN] + dYclusBoxMax)))) {
0385 seeds[NNN].push_back((*it2));
0386 }
0387 }
0388 }
0389 }
0390
0391
0392 int secondCh = 0;
0393 DetId detId_prev;
0394 if (seeds[NNN].size() > 1) {
0395 for (unsigned int iRH = 0; iRH < seeds[NNN].size(); ++iRH) {
0396 if (iRH && detId_prev != seeds[NNN][iRH]->hit()->geographicalId()) {
0397 ++secondCh;
0398 break;
0399 }
0400 detId_prev = seeds[NNN][iRH]->hit()->geographicalId();
0401 }
0402 }
0403
0404 if (useRPCs && !secondCh && !tooCloseClusters) {
0405 for (MuonRecHitContainer::const_iterator it2 = muonRecHits_RPC.begin(); it2 != muonRecHits_RPC.end(); ++it2) {
0406 if (((*it2)->globalPosition().phi() < seed_maxX[NNN] + dXclusBoxMax) &&
0407 ((*it2)->globalPosition().phi() > seed_minX[NNN] - dXclusBoxMax)) {
0408 if (!((((*it2)->globalPosition().theta() + 0.3) < (seed_minY[NNN] - dYclusBoxMax) &&
0409 ((*it2)->globalPosition().theta() - 0.3) < (seed_minY[NNN] - dYclusBoxMax)) ||
0410 (((*it2)->globalPosition().theta() + 0.3) > (seed_maxY[NNN] + dYclusBoxMax) &&
0411
0412 ((*it2)->globalPosition().theta() - 0.3) > (seed_maxY[NNN] + dYclusBoxMax)))) {
0413 seeds[NNN].push_back((*it2));
0414 }
0415 }
0416 }
0417 }
0418 }
0419
0420
0421
0422
0423 for (unsigned int NNN = 0; NNN < seeds.size(); ++NNN) {
0424 if (running_meanX[NNN] == 999999.)
0425 continue;
0426
0427 segments_clusters.push_back(seeds[NNN]);
0428 }
0429 }
0430
0431 bool SETPatternRecognition::segmentCleaning(const DetId& detId,
0432 const LocalPoint& localPosition,
0433 const LocalError& localError,
0434 const LocalVector& localDirection,
0435 const LocalError& localDirectionError,
0436 const double& chi2,
0437 const int& ndf) {
0438
0439 bool dropTheSegment = true;
0440 const GeomDet* geomDet = theService->trackingGeometry()->idToDet(detId);
0441
0442 bool insideCh = geomDet->surface().bounds().inside(localPosition, localError, outsideChamberErrorScale);
0443
0444
0445
0446 bool parallelSegment = fabs(localDirection.z()) > minLocalSegmentAngle ? false : true;
0447
0448 if (insideCh && !parallelSegment) {
0449 dropTheSegment = false;
0450 }
0451
0452
0453 return dropTheSegment;
0454 }