File indexing completed on 2024-04-06 12:26:57
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 #include "RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGeneratorFromL1T.h"
0022
0023
0024 #include "FWCore/Framework/interface/ConsumesCollector.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0030
0031 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0032 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
0033 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0034 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0035 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0036 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0037 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0038
0039 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0040 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0041
0042 using namespace std;
0043 using namespace edm;
0044 using namespace l1t;
0045
0046
0047 bool sortByL1Pt(L2MuonTrajectorySeed &A, L2MuonTrajectorySeed &B) {
0048 l1t::MuonRef Ref_L1A = A.l1tParticle();
0049 l1t::MuonRef Ref_L1B = B.l1tParticle();
0050 return (Ref_L1A->pt() > Ref_L1B->pt());
0051 };
0052
0053 bool sortByL1QandPt(L2MuonTrajectorySeed &A, L2MuonTrajectorySeed &B) {
0054 l1t::MuonRef Ref_L1A = A.l1tParticle();
0055 l1t::MuonRef Ref_L1B = B.l1tParticle();
0056
0057
0058 if (Ref_L1A->hwQual() > Ref_L1B->hwQual())
0059 return true;
0060 if (Ref_L1A->hwQual() < Ref_L1B->hwQual())
0061 return false;
0062
0063
0064 return (Ref_L1A->pt() > Ref_L1B->pt());
0065 };
0066
0067
0068 L2MuonSeedGeneratorFromL1T::L2MuonSeedGeneratorFromL1T(const edm::ParameterSet &iConfig)
0069 : theSource(iConfig.getParameter<InputTag>("InputObjects")),
0070 theL1GMTReadoutCollection(iConfig.getParameter<InputTag>("GMTReadoutCollection")),
0071 thePropagatorName(iConfig.getParameter<string>("Propagator")),
0072 theL1MinPt(iConfig.getParameter<double>("L1MinPt")),
0073 theL1MaxEta(iConfig.getParameter<double>("L1MaxEta")),
0074 theL1MinQuality(iConfig.getParameter<unsigned int>("L1MinQuality")),
0075 theMinPtBarrel(iConfig.getParameter<double>("SetMinPtBarrelTo")),
0076 theMinPtEndcap(iConfig.getParameter<double>("SetMinPtEndcapTo")),
0077 useOfflineSeed(iConfig.getUntrackedParameter<bool>("UseOfflineSeed", false)),
0078 useUnassociatedL1(iConfig.getParameter<bool>("UseUnassociatedL1")),
0079 matchingDR(iConfig.getParameter<std::vector<double>>("MatchDR")),
0080 etaBins(iConfig.getParameter<std::vector<double>>("EtaMatchingBins")),
0081 centralBxOnly_(iConfig.getParameter<bool>("CentralBxOnly")),
0082 matchType(iConfig.getParameter<unsigned int>("MatchType")),
0083 sortType(iConfig.getParameter<unsigned int>("SortType")) {
0084 muCollToken_ = consumes<MuonBxCollection>(theSource);
0085
0086 if (useOfflineSeed) {
0087 theOfflineSeedLabel = iConfig.getUntrackedParameter<InputTag>("OfflineSeedLabel");
0088 offlineSeedToken_ = consumes<edm::View<TrajectorySeed>>(theOfflineSeedLabel);
0089
0090
0091 if (matchingDR.size() != etaBins.size() - 1) {
0092 throw cms::Exception("Configuration") << "Size of MatchDR "
0093 << "does not match number of eta bins." << endl;
0094 }
0095 }
0096
0097 if (matchType > 4 || sortType > 3) {
0098 throw cms::Exception("Configuration") << "Wrong MatchType or SortType" << endl;
0099 }
0100
0101
0102 ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
0103
0104
0105 theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0106
0107
0108 theEstimator = new Chi2MeasurementEstimator(10000.);
0109
0110 produces<L2MuonTrajectorySeedCollection>();
0111 }
0112
0113
0114 L2MuonSeedGeneratorFromL1T::~L2MuonSeedGeneratorFromL1T() {
0115 if (theService)
0116 delete theService;
0117 if (theEstimator)
0118 delete theEstimator;
0119 }
0120
0121 void L2MuonSeedGeneratorFromL1T::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0122 edm::ParameterSetDescription desc;
0123 desc.add<edm::InputTag>("GMTReadoutCollection", edm::InputTag(""));
0124 desc.add<edm::InputTag>("InputObjects", edm::InputTag("hltGmtStage2Digis"));
0125 desc.add<string>("Propagator", "");
0126 desc.add<double>("L1MinPt", -1.);
0127 desc.add<double>("L1MaxEta", 5.0);
0128 desc.add<unsigned int>("L1MinQuality", 0);
0129 desc.add<double>("SetMinPtBarrelTo", 3.5);
0130 desc.add<double>("SetMinPtEndcapTo", 1.0);
0131 desc.addUntracked<bool>("UseOfflineSeed", false);
0132 desc.add<bool>("UseUnassociatedL1", true);
0133 desc.add<std::vector<double>>("MatchDR", {0.3});
0134 desc.add<std::vector<double>>("EtaMatchingBins", {0., 2.5});
0135 desc.add<bool>("CentralBxOnly", true);
0136 desc.add<unsigned int>("MatchType", 0)
0137 ->setComment(
0138 "MatchType : 0 Old matching, 1 L1 Order(1to1), 2 Min dR(1to1), 3 Higher Q(1to1), 4 All matched L1");
0139 desc.add<unsigned int>("SortType", 0)->setComment("SortType : 0 not sort, 1 Pt, 2 Q and Pt");
0140 desc.addUntracked<edm::InputTag>("OfflineSeedLabel", edm::InputTag(""));
0141
0142 edm::ParameterSetDescription psd0;
0143 psd0.addUntracked<std::vector<std::string>>("Propagators", {"SteppingHelixPropagatorAny"});
0144 psd0.add<bool>("RPCLayers", true);
0145 psd0.addUntracked<bool>("UseMuonNavigation", true);
0146 desc.add<edm::ParameterSetDescription>("ServiceParameters", psd0);
0147 descriptions.add("L2MuonSeedGeneratorFromL1T", desc);
0148 }
0149
0150 void L2MuonSeedGeneratorFromL1T::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0151 const std::string metname = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
0152 MuonPatternRecoDumper debug;
0153
0154 auto output = std::make_unique<L2MuonTrajectorySeedCollection>();
0155
0156 edm::Handle<MuonBxCollection> muColl;
0157 iEvent.getByToken(muCollToken_, muColl);
0158 LogDebug(metname) << "Number of muons " << muColl->size() << endl;
0159
0160
0161 if (matchType == 0) {
0162 edm::Handle<edm::View<TrajectorySeed>> offlineSeedHandle;
0163 vector<int> offlineSeedMap;
0164 if (useOfflineSeed) {
0165 iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
0166 LogDebug(metname) << "Number of offline seeds " << offlineSeedHandle->size() << endl;
0167 offlineSeedMap = vector<int>(offlineSeedHandle->size(), 0);
0168 }
0169
0170 for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) {
0171 if (centralBxOnly_ && (ibx != 0))
0172 continue;
0173 for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++) {
0174 unsigned int quality = it->hwQual();
0175 int valid_charge = it->hwChargeValid();
0176
0177 float pt = it->pt();
0178 float eta = it->eta();
0179 float theta = 2 * atan(exp(-eta));
0180 float phi = it->phi();
0181 int charge = it->charge();
0182
0183 if (!valid_charge)
0184 charge = 0;
0185
0186
0187
0188 int link = 36 + (int)(it->tfMuonIndex() / 3.);
0189 bool barrel = true;
0190 if ((link >= 36 && link <= 41) || (link >= 66 && link <= 71))
0191 barrel = false;
0192
0193 if (pt < theL1MinPt || fabs(eta) > theL1MaxEta)
0194 continue;
0195
0196 LogDebug(metname) << "New L2 Muon Seed";
0197 LogDebug(metname) << "Pt = " << pt << " GeV/c";
0198 LogDebug(metname) << "eta = " << eta;
0199 LogDebug(metname) << "theta = " << theta << " rad";
0200 LogDebug(metname) << "phi = " << phi << " rad";
0201 LogDebug(metname) << "charge = " << charge;
0202 LogDebug(metname) << "In Barrel? = " << barrel;
0203
0204 if (quality <= theL1MinQuality)
0205 continue;
0206 LogDebug(metname) << "quality = " << quality;
0207
0208
0209 theService->update(iSetup);
0210
0211 const DetLayer *detLayer = nullptr;
0212 float radius = 0.;
0213
0214 CLHEP::Hep3Vector vec(0., 1., 0.);
0215 vec.setTheta(theta);
0216 vec.setPhi(phi);
0217
0218 DetId theid;
0219
0220 if (barrel) {
0221 LogDebug(metname) << "The seed is in the barrel";
0222
0223
0224 theid = DTChamberId(0, 2, 0);
0225 detLayer = theService->detLayerGeometry()->idToLayer(theid);
0226 LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0227
0228 const BoundSurface *sur = &(detLayer->surface());
0229 const BoundCylinder *bc = dynamic_cast<const BoundCylinder *>(sur);
0230
0231 radius = fabs(bc->radius() / sin(theta));
0232
0233 LogDebug(metname) << "radius " << radius;
0234
0235 if (pt < theMinPtBarrel)
0236 pt = theMinPtBarrel;
0237 } else {
0238 LogDebug(metname) << "The seed is in the endcap";
0239
0240
0241 theid = theta < Geom::pi() / 2. ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
0242
0243 detLayer = theService->detLayerGeometry()->idToLayer(theid);
0244 LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0245
0246 radius = fabs(detLayer->position().z() / cos(theta));
0247
0248 if (pt < theMinPtEndcap)
0249 pt = theMinPtEndcap;
0250 }
0251
0252 vec.setMag(radius);
0253
0254 GlobalPoint pos(vec.x(), vec.y(), vec.z());
0255
0256 GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));
0257
0258 GlobalTrajectoryParameters param(pos, mom, charge, &*theService->magneticField());
0259 AlgebraicSymMatrix55 mat;
0260
0261 mat[0][0] = (0.25 / pt) * (0.25 / pt);
0262 if (!barrel)
0263 mat[0][0] = (0.4 / pt) * (0.4 / pt);
0264
0265
0266 if (!valid_charge)
0267 mat[0][0] = (1. / pt) * (1. / pt);
0268
0269 mat[1][1] = 0.05 * 0.05;
0270 mat[2][2] = 0.2 * 0.2;
0271 mat[3][3] = 20. * 20.;
0272 mat[4][4] = 20. * 20.;
0273
0274 CurvilinearTrajectoryError error(mat);
0275
0276 const FreeTrajectoryState state(param, error);
0277
0278 LogDebug(metname) << "Free trajectory State from the parameters";
0279 LogDebug(metname) << debug.dumpFTS(state);
0280
0281
0282 TrajectoryStateOnSurface tsos =
0283 theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());
0284
0285 LogDebug(metname) << "State after the propagation on the layer";
0286 LogDebug(metname) << debug.dumpLayer(detLayer);
0287 LogDebug(metname) << debug.dumpTSOS(tsos);
0288
0289 double dRcone = matchingDR[0];
0290 if (fabs(eta) < etaBins.back()) {
0291 std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(eta));
0292 dRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0293 }
0294
0295 if (tsos.isValid()) {
0296 edm::OwnVector<TrackingRecHit> container;
0297
0298 if (useOfflineSeed && (!valid_charge || charge == 0)) {
0299 const TrajectorySeed *assoOffseed =
0300 associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, tsos, dRcone);
0301
0302 if (assoOffseed != nullptr) {
0303 PTrajectoryStateOnDet const &seedTSOS = assoOffseed->startingState();
0304 for (auto const &recHit : assoOffseed->recHits()) {
0305 container.push_back(recHit);
0306 }
0307 output->push_back(
0308 L2MuonTrajectorySeed(seedTSOS,
0309 container,
0310 alongMomentum,
0311 MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0312 } else {
0313 if (useUnassociatedL1) {
0314
0315 PTrajectoryStateOnDet const &seedTSOS = trajectoryStateTransform::persistentState(tsos, theid.rawId());
0316 output->push_back(
0317 L2MuonTrajectorySeed(seedTSOS,
0318 container,
0319 alongMomentum,
0320 MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0321 }
0322 }
0323 } else if (useOfflineSeed && valid_charge) {
0324
0325 std::vector<pair<const GeomDet *, TrajectoryStateOnSurface>> detsWithStates =
0326 detLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
0327
0328 if (detsWithStates.empty() && barrel) {
0329
0330 DetId fallback_id;
0331 theta < Geom::pi() / 2. ? fallback_id = CSCDetId(1, 2, 0, 0, 0) : fallback_id = CSCDetId(2, 2, 0, 0, 0);
0332 const DetLayer *ME2DetLayer = theService->detLayerGeometry()->idToLayer(fallback_id);
0333
0334 tsos = theService->propagator(thePropagatorName)->propagate(state, ME2DetLayer->surface());
0335 detsWithStates =
0336 ME2DetLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
0337 }
0338
0339 if (!detsWithStates.empty()) {
0340 TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
0341 const GeomDet *newTSOSDet = detsWithStates.front().first;
0342
0343 LogDebug(metname) << "Most compatible det";
0344 LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
0345
0346 if (newTSOS.isValid()) {
0347 LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag()
0348 << ", phi=" << newTSOS.globalPosition().phi()
0349 << ", eta=" << newTSOS.globalPosition().eta() << ")";
0350 LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge() * newTSOS.globalMomentum().perp()
0351 << ", phi=" << newTSOS.globalMomentum().phi()
0352 << ", eta=" << newTSOS.globalMomentum().eta() << ")";
0353
0354 const TrajectorySeed *assoOffseed =
0355 associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, newTSOS, dRcone);
0356
0357 if (assoOffseed != nullptr) {
0358 PTrajectoryStateOnDet const &seedTSOS = assoOffseed->startingState();
0359 for (auto const &recHit : assoOffseed->recHits()) {
0360 container.push_back(recHit);
0361 }
0362 output->push_back(
0363 L2MuonTrajectorySeed(seedTSOS,
0364 container,
0365 alongMomentum,
0366 MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0367 } else {
0368 if (useUnassociatedL1) {
0369
0370 PTrajectoryStateOnDet const &seedTSOS =
0371 trajectoryStateTransform::persistentState(newTSOS, newTSOSDet->geographicalId().rawId());
0372 output->push_back(
0373 L2MuonTrajectorySeed(seedTSOS,
0374 container,
0375 alongMomentum,
0376 MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0377 }
0378 }
0379 }
0380 }
0381 } else {
0382
0383 PTrajectoryStateOnDet const &seedTSOS = trajectoryStateTransform::persistentState(tsos, theid.rawId());
0384 output->push_back(L2MuonTrajectorySeed(
0385 seedTSOS, container, alongMomentum, MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0386 }
0387 }
0388 }
0389 }
0390
0391 }
0392
0393
0394 else if (matchType > 0) {
0395 unsigned int nMuColl = muColl->size();
0396
0397 vector<vector<double>> dRmtx;
0398 vector<vector<const TrajectorySeed *>> selOffseeds;
0399
0400 edm::Handle<edm::View<TrajectorySeed>> offlineSeedHandle;
0401 if (useOfflineSeed) {
0402 iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
0403 unsigned int nOfflineSeed = offlineSeedHandle->size();
0404 LogDebug(metname) << "Number of offline seeds " << nOfflineSeed << endl;
0405
0406
0407 dRmtx = vector<vector<double>>(nMuColl, vector<double>(nOfflineSeed, 999.0));
0408 selOffseeds =
0409 vector<vector<const TrajectorySeed *>>(nMuColl, vector<const TrajectorySeed *>(nOfflineSeed, nullptr));
0410 }
0411
0412
0413 bool isAsso = false;
0414
0415
0416 for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) {
0417 if (centralBxOnly_ && (ibx != 0))
0418 continue;
0419 for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++) {
0420
0421 unsigned int imu = distance(muColl->begin(muColl->getFirstBX()), it);
0422
0423 unsigned int quality = it->hwQual();
0424 int valid_charge = it->hwChargeValid();
0425
0426 float pt = it->pt();
0427 float eta = it->eta();
0428 float theta = 2 * atan(exp(-eta));
0429 float phi = it->phi();
0430 int charge = it->charge();
0431
0432 if (!valid_charge)
0433 charge = 0;
0434
0435
0436
0437 int link = 36 + (int)(it->tfMuonIndex() / 3.);
0438 bool barrel = true;
0439 if ((link >= 36 && link <= 41) || (link >= 66 && link <= 71))
0440 barrel = false;
0441
0442 if (pt < theL1MinPt || fabs(eta) > theL1MaxEta)
0443 continue;
0444
0445 LogDebug(metname) << "New L2 Muon Seed";
0446 LogDebug(metname) << "Pt = " << pt << " GeV/c";
0447 LogDebug(metname) << "eta = " << eta;
0448 LogDebug(metname) << "theta = " << theta << " rad";
0449 LogDebug(metname) << "phi = " << phi << " rad";
0450 LogDebug(metname) << "charge = " << charge;
0451 LogDebug(metname) << "In Barrel? = " << barrel;
0452
0453 if (quality <= theL1MinQuality)
0454 continue;
0455 LogDebug(metname) << "quality = " << quality;
0456
0457
0458 theService->update(iSetup);
0459
0460 const DetLayer *detLayer = nullptr;
0461 float radius = 0.;
0462
0463 CLHEP::Hep3Vector vec(0., 1., 0.);
0464 vec.setTheta(theta);
0465 vec.setPhi(phi);
0466
0467 DetId theid;
0468
0469 if (barrel) {
0470 LogDebug(metname) << "The seed is in the barrel";
0471
0472
0473 theid = DTChamberId(0, 2, 0);
0474 detLayer = theService->detLayerGeometry()->idToLayer(theid);
0475 LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0476
0477 const BoundSurface *sur = &(detLayer->surface());
0478 const BoundCylinder *bc = dynamic_cast<const BoundCylinder *>(sur);
0479
0480 radius = fabs(bc->radius() / sin(theta));
0481
0482 LogDebug(metname) << "radius " << radius;
0483
0484 if (pt < theMinPtBarrel)
0485 pt = theMinPtBarrel;
0486 } else {
0487 LogDebug(metname) << "The seed is in the endcap";
0488
0489
0490 theid = theta < Geom::pi() / 2. ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
0491
0492 detLayer = theService->detLayerGeometry()->idToLayer(theid);
0493 LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0494
0495 radius = fabs(detLayer->position().z() / cos(theta));
0496
0497 if (pt < theMinPtEndcap)
0498 pt = theMinPtEndcap;
0499 }
0500
0501 vec.setMag(radius);
0502
0503 GlobalPoint pos(vec.x(), vec.y(), vec.z());
0504
0505 GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));
0506
0507 GlobalTrajectoryParameters param(pos, mom, charge, &*theService->magneticField());
0508 AlgebraicSymMatrix55 mat;
0509
0510 mat[0][0] = (0.25 / pt) * (0.25 / pt);
0511 if (!barrel)
0512 mat[0][0] = (0.4 / pt) * (0.4 / pt);
0513
0514
0515 if (!valid_charge)
0516 mat[0][0] = (1. / pt) * (1. / pt);
0517
0518 mat[1][1] = 0.05 * 0.05;
0519 mat[2][2] = 0.2 * 0.2;
0520 mat[3][3] = 20. * 20.;
0521 mat[4][4] = 20. * 20.;
0522
0523 CurvilinearTrajectoryError error(mat);
0524
0525 const FreeTrajectoryState state(param, error);
0526
0527 LogDebug(metname) << "Free trajectory State from the parameters";
0528 LogDebug(metname) << debug.dumpFTS(state);
0529
0530
0531 TrajectoryStateOnSurface tsos =
0532 theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());
0533
0534 LogDebug(metname) << "State after the propagation on the layer";
0535 LogDebug(metname) << debug.dumpLayer(detLayer);
0536 LogDebug(metname) << debug.dumpTSOS(tsos);
0537
0538 double dRcone = matchingDR[0];
0539 if (fabs(eta) < etaBins.back()) {
0540 std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(eta));
0541 dRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0542 }
0543
0544 if (tsos.isValid()) {
0545 edm::OwnVector<TrackingRecHit> container;
0546
0547 if (useOfflineSeed) {
0548 if ((!valid_charge || charge == 0)) {
0549 bool isAssoOffseed = isAssociateOfflineSeedToL1(offlineSeedHandle, dRmtx, tsos, imu, selOffseeds, dRcone);
0550
0551 if (isAssoOffseed) {
0552 isAsso = true;
0553 }
0554
0555
0556 else {
0557 if (useUnassociatedL1) {
0558
0559 PTrajectoryStateOnDet const &seedTSOS =
0560 trajectoryStateTransform::persistentState(tsos, theid.rawId());
0561 output->push_back(
0562 L2MuonTrajectorySeed(seedTSOS,
0563 container,
0564 alongMomentum,
0565 MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0566 }
0567 }
0568
0569 }
0570
0571 else if (valid_charge) {
0572
0573 std::vector<pair<const GeomDet *, TrajectoryStateOnSurface>> detsWithStates =
0574 detLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
0575
0576 if (detsWithStates.empty() && barrel) {
0577
0578 DetId fallback_id;
0579 theta < Geom::pi() / 2. ? fallback_id = CSCDetId(1, 2, 0, 0, 0) : fallback_id = CSCDetId(2, 2, 0, 0, 0);
0580 const DetLayer *ME2DetLayer = theService->detLayerGeometry()->idToLayer(fallback_id);
0581
0582 tsos = theService->propagator(thePropagatorName)->propagate(state, ME2DetLayer->surface());
0583 detsWithStates =
0584 ME2DetLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
0585 }
0586
0587 if (!detsWithStates.empty()) {
0588 TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
0589 const GeomDet *newTSOSDet = detsWithStates.front().first;
0590
0591 LogDebug(metname) << "Most compatible det";
0592 LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
0593
0594 if (newTSOS.isValid()) {
0595 LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag()
0596 << ", phi=" << newTSOS.globalPosition().phi()
0597 << ", eta=" << newTSOS.globalPosition().eta() << ")";
0598 LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge() * newTSOS.globalMomentum().perp()
0599 << ", phi=" << newTSOS.globalMomentum().phi()
0600 << ", eta=" << newTSOS.globalMomentum().eta() << ")";
0601
0602 bool isAssoOffseed =
0603 isAssociateOfflineSeedToL1(offlineSeedHandle, dRmtx, newTSOS, imu, selOffseeds, dRcone);
0604
0605 if (isAssoOffseed) {
0606 isAsso = true;
0607 }
0608
0609
0610 else {
0611 if (useUnassociatedL1) {
0612
0613 PTrajectoryStateOnDet const &seedTSOS =
0614 trajectoryStateTransform::persistentState(newTSOS, newTSOSDet->geographicalId().rawId());
0615 output->push_back(
0616 L2MuonTrajectorySeed(seedTSOS,
0617 container,
0618 alongMomentum,
0619 MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0620 }
0621 }
0622 }
0623 }
0624 }
0625 }
0626
0627 else {
0628
0629 PTrajectoryStateOnDet const &seedTSOS = trajectoryStateTransform::persistentState(tsos, theid.rawId());
0630 output->push_back(L2MuonTrajectorySeed(
0631 seedTSOS, container, alongMomentum, MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0632 }
0633 }
0634 }
0635 }
0636
0637
0638 if (useOfflineSeed && isAsso) {
0639 unsigned int nOfflineSeed1 = offlineSeedHandle->size();
0640 unsigned int nL1;
0641 unsigned int i, j;
0642
0643 if (matchType == 1) {
0644 vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);
0645
0646 for (nL1 = 0; nL1 < nMuColl; ++nL1) {
0647 double tempDR = 99.;
0648 unsigned int theOffs = 0;
0649
0650 for (j = 0; j < nOfflineSeed1; ++j) {
0651 if (removed_col[j])
0652 continue;
0653 if (tempDR > dRmtx[nL1][j]) {
0654 tempDR = dRmtx[nL1][j];
0655 theOffs = j;
0656 }
0657 }
0658
0659 auto it = muColl->begin(muColl->getFirstBX()) + nL1;
0660
0661 double newDRcone = matchingDR[0];
0662 if (fabs(it->eta()) < etaBins.back()) {
0663 std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
0664 newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0665 }
0666
0667 if (!(tempDR < newDRcone))
0668 continue;
0669
0670
0671 removed_col[theOffs] = true;
0672
0673 if (selOffseeds[nL1][theOffs] != nullptr) {
0674
0675 edm::OwnVector<TrackingRecHit> newContainer;
0676
0677 PTrajectoryStateOnDet const &seedTSOS = selOffseeds[nL1][theOffs]->startingState();
0678 for (auto const &recHit : selOffseeds[nL1][theOffs]->recHits()) {
0679 newContainer.push_back(recHit);
0680 }
0681 output->push_back(L2MuonTrajectorySeed(seedTSOS,
0682 newContainer,
0683 alongMomentum,
0684 MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0685 }
0686 }
0687 }
0688
0689 else if (matchType == 2) {
0690 vector<bool> removed_row = vector<bool>(nMuColl, false);
0691 vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);
0692
0693 for (nL1 = 0; nL1 < nMuColl; ++nL1) {
0694 double tempDR = 99.;
0695 unsigned int theL1 = 0;
0696 unsigned int theOffs = 0;
0697
0698 for (i = 0; i < nMuColl; ++i) {
0699 if (removed_row[i])
0700 continue;
0701 for (j = 0; j < nOfflineSeed1; ++j) {
0702 if (removed_col[j])
0703 continue;
0704 if (tempDR > dRmtx[i][j]) {
0705 tempDR = dRmtx[i][j];
0706 theL1 = i;
0707 theOffs = j;
0708 }
0709 }
0710 }
0711
0712 auto it = muColl->begin(muColl->getFirstBX()) + theL1;
0713
0714 double newDRcone = matchingDR[0];
0715 if (fabs(it->eta()) < etaBins.back()) {
0716 std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
0717 newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0718 }
0719
0720 if (!(tempDR < newDRcone))
0721 continue;
0722
0723
0724 removed_col[theOffs] = true;
0725 removed_row[theL1] = true;
0726
0727 if (selOffseeds[theL1][theOffs] != nullptr) {
0728
0729 edm::OwnVector<TrackingRecHit> newContainer;
0730
0731 PTrajectoryStateOnDet const &seedTSOS = selOffseeds[theL1][theOffs]->startingState();
0732 for (auto const &recHit : selOffseeds[theL1][theOffs]->recHits()) {
0733 newContainer.push_back(recHit);
0734 }
0735 output->push_back(L2MuonTrajectorySeed(seedTSOS,
0736 newContainer,
0737 alongMomentum,
0738 MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0739 }
0740 }
0741 }
0742
0743 else if (matchType == 3) {
0744 vector<bool> removed_row = vector<bool>(nMuColl, false);
0745 vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);
0746
0747 for (nL1 = 0; nL1 < nMuColl; ++nL1) {
0748 double tempDR = 99.;
0749 unsigned int theL1 = 0;
0750 unsigned int theOffs = 0;
0751 auto theit = muColl->begin(muColl->getFirstBX());
0752
0753
0754 for (i = 0; i < nMuColl; ++i) {
0755 if (removed_row[i])
0756 continue;
0757 theit = muColl->begin(muColl->getFirstBX()) + i;
0758 if (theit->hwQual() > 10) {
0759 for (j = 0; j < nOfflineSeed1; ++j) {
0760 if (removed_col[j])
0761 continue;
0762 if (tempDR > dRmtx[i][j]) {
0763 tempDR = dRmtx[i][j];
0764 theL1 = i;
0765 theOffs = j;
0766 }
0767 }
0768 }
0769 }
0770
0771
0772 if (tempDR == 99.) {
0773 for (i = 0; i < nMuColl; ++i) {
0774 if (removed_row[i])
0775 continue;
0776 theit = muColl->begin(muColl->getFirstBX()) + i;
0777 if ((theit->hwQual() <= 10) && (theit->hwQual() > 6)) {
0778 for (j = 0; j < nOfflineSeed1; ++j) {
0779 if (removed_col[j])
0780 continue;
0781 if (tempDR > dRmtx[i][j]) {
0782 tempDR = dRmtx[i][j];
0783 theL1 = i;
0784 theOffs = j;
0785 }
0786 }
0787 }
0788 }
0789 }
0790
0791
0792 if (tempDR == 99.) {
0793 for (i = 0; i < nMuColl; ++i) {
0794 if (removed_row[i])
0795 continue;
0796 theit = muColl->begin(muColl->getFirstBX()) + i;
0797 if (theit->hwQual() <= 6) {
0798 for (j = 0; j < nOfflineSeed1; ++j) {
0799 if (removed_col[j])
0800 continue;
0801 if (tempDR > dRmtx[i][j]) {
0802 tempDR = dRmtx[i][j];
0803 theL1 = i;
0804 theOffs = j;
0805 }
0806 }
0807 }
0808 }
0809 }
0810
0811 auto it = muColl->begin(muColl->getFirstBX()) + theL1;
0812
0813 double newDRcone = matchingDR[0];
0814 if (fabs(it->eta()) < etaBins.back()) {
0815 std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
0816 newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0817 }
0818
0819 if (!(tempDR < newDRcone))
0820 continue;
0821
0822
0823 removed_col[theOffs] = true;
0824 removed_row[theL1] = true;
0825
0826 if (selOffseeds[theL1][theOffs] != nullptr) {
0827
0828 edm::OwnVector<TrackingRecHit> newContainer;
0829
0830 PTrajectoryStateOnDet const &seedTSOS = selOffseeds[theL1][theOffs]->startingState();
0831 for (auto const &recHit : selOffseeds[theL1][theOffs]->recHits()) {
0832 newContainer.push_back(recHit);
0833 }
0834 output->push_back(L2MuonTrajectorySeed(seedTSOS,
0835 newContainer,
0836 alongMomentum,
0837 MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0838 }
0839 }
0840 }
0841
0842 else if (matchType == 4) {
0843 for (i = 0; i < nMuColl; ++i) {
0844 auto it = muColl->begin(muColl->getFirstBX()) + i;
0845
0846 double tempDR = 99.;
0847 unsigned int theOffs = 0;
0848
0849 double newDRcone = matchingDR[0];
0850 if (fabs(it->eta()) < etaBins.back()) {
0851 std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
0852 newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
0853 }
0854
0855 for (j = 0; j < nOfflineSeed1; ++j) {
0856 if (tempDR > dRmtx[i][j]) {
0857 tempDR = dRmtx[i][j];
0858 theOffs = j;
0859 }
0860 }
0861
0862 if (!(tempDR < newDRcone))
0863 continue;
0864
0865 if (selOffseeds[i][theOffs] != nullptr) {
0866
0867 edm::OwnVector<TrackingRecHit> newContainer;
0868
0869 PTrajectoryStateOnDet const &seedTSOS = selOffseeds[i][theOffs]->startingState();
0870 for (auto const &recHit : selOffseeds[i][theOffs]->recHits()) {
0871 newContainer.push_back(recHit);
0872 }
0873 output->push_back(L2MuonTrajectorySeed(seedTSOS,
0874 newContainer,
0875 alongMomentum,
0876 MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
0877 }
0878 }
0879 }
0880
0881 }
0882
0883 }
0884
0885
0886 if (sortType == 1) {
0887 sort(output->begin(), output->end(), sortByL1Pt);
0888 }
0889
0890
0891 else if (sortType == 2) {
0892 sort(output->begin(), output->end(), sortByL1QandPt);
0893 }
0894
0895 iEvent.put(std::move(output));
0896 }
0897
0898
0899 const TrajectorySeed *L2MuonSeedGeneratorFromL1T::associateOfflineSeedToL1(
0900 edm::Handle<edm::View<TrajectorySeed>> &offseeds,
0901 std::vector<int> &offseedMap,
0902 TrajectoryStateOnSurface &newTsos,
0903 double dRcone) {
0904 const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
0905 MuonPatternRecoDumper debugtmp;
0906
0907 edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
0908 const TrajectorySeed *selOffseed = nullptr;
0909 double bestDr = 99999.;
0910 unsigned int nOffseed(0);
0911 int lastOffseed(-1);
0912
0913 for (offseed = offseeds->begin(); offseed != endOffseed; ++offseed, ++nOffseed) {
0914 if (offseedMap[nOffseed] != 0)
0915 continue;
0916 GlobalPoint glbPos = theService->trackingGeometry()
0917 ->idToDet(offseed->startingState().detId())
0918 ->surface()
0919 .toGlobal(offseed->startingState().parameters().position());
0920 GlobalVector glbMom = theService->trackingGeometry()
0921 ->idToDet(offseed->startingState().detId())
0922 ->surface()
0923 .toGlobal(offseed->startingState().parameters().momentum());
0924
0925
0926 double preDr = deltaR(newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi());
0927 if (preDr > 1.0)
0928 continue;
0929
0930 const FreeTrajectoryState offseedFTS(
0931 glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField());
0932 TrajectoryStateOnSurface offseedTsos =
0933 theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
0934 LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
0935 LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
0936 LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" << offseedFTS.position().phi()
0937 << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
0938 LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge() * offseedFTS.momentum().perp()
0939 << ", phi=" << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")"
0940 << std::endl
0941 << std::endl;
0942
0943 if (offseedTsos.isValid()) {
0944 LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
0945 LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag()
0946 << ", phi=" << offseedTsos.globalPosition().phi()
0947 << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
0948 LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge() * offseedTsos.globalMomentum().perp()
0949 << ", phi=" << offseedTsos.globalMomentum().phi()
0950 << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl
0951 << std::endl;
0952 double newDr = deltaR(newTsos.globalPosition().eta(),
0953 newTsos.globalPosition().phi(),
0954 offseedTsos.globalPosition().eta(),
0955 offseedTsos.globalPosition().phi());
0956 LogDebug(metlabel) << " -- DR = " << newDr << std::endl;
0957 if (newDr < dRcone && newDr < bestDr) {
0958 LogDebug(metlabel) << " --> OK! " << newDr << std::endl << std::endl;
0959 selOffseed = &*offseed;
0960 bestDr = newDr;
0961 offseedMap[nOffseed] = 1;
0962 if (lastOffseed > -1)
0963 offseedMap[lastOffseed] = 0;
0964 lastOffseed = nOffseed;
0965 } else {
0966 LogDebug(metlabel) << " --> Rejected. " << newDr << std::endl << std::endl;
0967 }
0968 } else {
0969 LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
0970 }
0971 }
0972
0973 return selOffseed;
0974 }
0975
0976 bool L2MuonSeedGeneratorFromL1T::isAssociateOfflineSeedToL1(
0977 edm::Handle<edm::View<TrajectorySeed>> &offseeds,
0978 std::vector<std::vector<double>> &dRmtx,
0979 TrajectoryStateOnSurface &newTsos,
0980 unsigned int imu,
0981 std::vector<std::vector<const TrajectorySeed *>> &selOffseeds,
0982 double dRcone) {
0983 const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
0984 bool isAssociated = false;
0985 MuonPatternRecoDumper debugtmp;
0986
0987 edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
0988 unsigned int nOffseed(0);
0989
0990 for (offseed = offseeds->begin(); offseed != endOffseed; ++offseed, ++nOffseed) {
0991 GlobalPoint glbPos = theService->trackingGeometry()
0992 ->idToDet(offseed->startingState().detId())
0993 ->surface()
0994 .toGlobal(offseed->startingState().parameters().position());
0995 GlobalVector glbMom = theService->trackingGeometry()
0996 ->idToDet(offseed->startingState().detId())
0997 ->surface()
0998 .toGlobal(offseed->startingState().parameters().momentum());
0999
1000
1001 double preDr = deltaR(newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi());
1002 if (preDr > 1.0)
1003 continue;
1004
1005 const FreeTrajectoryState offseedFTS(
1006 glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField());
1007 TrajectoryStateOnSurface offseedTsos =
1008 theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
1009 LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
1010 LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
1011 LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" << offseedFTS.position().phi()
1012 << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
1013 LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge() * offseedFTS.momentum().perp()
1014 << ", phi=" << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")"
1015 << std::endl
1016 << std::endl;
1017
1018 if (offseedTsos.isValid()) {
1019 LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
1020 LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag()
1021 << ", phi=" << offseedTsos.globalPosition().phi()
1022 << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
1023 LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge() * offseedTsos.globalMomentum().perp()
1024 << ", phi=" << offseedTsos.globalMomentum().phi()
1025 << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl
1026 << std::endl;
1027 double newDr = deltaR(newTsos.globalPosition().eta(),
1028 newTsos.globalPosition().phi(),
1029 offseedTsos.globalPosition().eta(),
1030 offseedTsos.globalPosition().phi());
1031
1032 LogDebug(metlabel) << " -- DR = " << newDr << std::endl;
1033 if (newDr < dRcone) {
1034 LogDebug(metlabel) << " --> OK! " << newDr << std::endl << std::endl;
1035
1036 dRmtx[imu][nOffseed] = newDr;
1037 selOffseeds[imu][nOffseed] = &*offseed;
1038
1039 isAssociated = true;
1040 } else {
1041 LogDebug(metlabel) << " --> Rejected. " << newDr << std::endl << std::endl;
1042 }
1043 } else {
1044 LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
1045 }
1046 }
1047
1048 return isAssociated;
1049 }