File indexing completed on 2025-01-27 02:50:31
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/Utilities/interface/InputTag.h"
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "FWCore/Framework/interface/ConsumesCollector.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0023
0024
0025 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeed.h"
0026 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeedCollection.h"
0027 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0028 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0029 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0030 #include "DataFormats/Common/interface/Handle.h"
0031 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
0032 #include "DataFormats/Math/interface/deltaR.h"
0033 #include "DataFormats/L1Trigger/interface/Muon.h"
0034 #include "DataFormats/L1TMuonPhase2/interface/TrackerMuon.h"
0035
0036 #include "CLHEP/Vector/ThreeVector.h"
0037 #include "Geometry/CommonDetUnit/interface/GeomDetEnumerators.h"
0038 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0039 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
0040 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0041 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0042 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0043 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0044 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0045
0046 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0047 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0048
0049 using namespace std;
0050 using namespace edm;
0051 using namespace l1t;
0052
0053 class L2MuonSeedGeneratorFromL1TkMu : public edm::stream::EDProducer<> {
0054 public:
0055
0056 explicit L2MuonSeedGeneratorFromL1TkMu(const edm::ParameterSet &);
0057
0058
0059 ~L2MuonSeedGeneratorFromL1TkMu() override = default;
0060
0061 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0062 void produce(edm::Event &, const edm::EventSetup &) override;
0063
0064 private:
0065 edm::InputTag source_;
0066 edm::EDGetTokenT<l1t::TrackerMuonCollection> muCollToken_;
0067
0068 edm::InputTag offlineSeedLabel_;
0069 edm::EDGetTokenT<edm::View<TrajectorySeed>> offlineSeedToken_;
0070
0071 std::string propagatorName_;
0072
0073 const double l1MinPt_;
0074 const double l1MaxEta_;
0075 const double minPtBarrel_;
0076 const double minPtEndcap_;
0077 const double minPL1Tk_;
0078 const double minPtL1TkBarrel_;
0079 const bool useOfflineSeed_;
0080 const bool useUnassociatedL1_;
0081 std::vector<double> matchingDR_;
0082 std::vector<double> etaBins_;
0083
0084
0085
0086
0087 static constexpr float etaBoundary_{1.1};
0088 static constexpr float distMB2_{550.};
0089 static constexpr float distME2_{850.};
0090 static constexpr float phiCorr0_{1.464};
0091 static constexpr float phiCorr1_{1.7};
0092 static constexpr float phiCorr2_{144.};
0093
0094
0095 std::unique_ptr<MuonServiceProxy> service_;
0096 std::unique_ptr<MeasurementEstimator> estimator_;
0097
0098 const TrajectorySeed *associateOfflineSeedToL1(edm::Handle<edm::View<TrajectorySeed>> &,
0099 std::vector<int> &,
0100 TrajectoryStateOnSurface &,
0101 double);
0102 };
0103
0104
0105 L2MuonSeedGeneratorFromL1TkMu::L2MuonSeedGeneratorFromL1TkMu(const edm::ParameterSet &iConfig)
0106 : source_(iConfig.getParameter<InputTag>("InputObjects")),
0107 muCollToken_(consumes(source_)),
0108 propagatorName_(iConfig.getParameter<string>("Propagator")),
0109 l1MinPt_(iConfig.getParameter<double>("L1MinPt")),
0110 l1MaxEta_(iConfig.getParameter<double>("L1MaxEta")),
0111 minPtBarrel_(iConfig.getParameter<double>("SetMinPtBarrelTo")),
0112 minPtEndcap_(iConfig.getParameter<double>("SetMinPtEndcapTo")),
0113 minPL1Tk_(iConfig.getParameter<double>("MinPL1Tk")),
0114 minPtL1TkBarrel_(iConfig.getParameter<double>("MinPtL1TkBarrel")),
0115 useOfflineSeed_(iConfig.getUntrackedParameter<bool>("UseOfflineSeed", false)),
0116 useUnassociatedL1_(iConfig.getParameter<bool>("UseUnassociatedL1")),
0117 matchingDR_(iConfig.getParameter<std::vector<double>>("MatchDR")),
0118 etaBins_(iConfig.getParameter<std::vector<double>>("EtaMatchingBins")) {
0119 if (useOfflineSeed_) {
0120 offlineSeedLabel_ = iConfig.getUntrackedParameter<InputTag>("OfflineSeedLabel");
0121 offlineSeedToken_ = consumes<edm::View<TrajectorySeed>>(offlineSeedLabel_);
0122
0123
0124 if (matchingDR_.size() != etaBins_.size() - 1) {
0125 throw cms::Exception("Configuration") << "Size of MatchDR "
0126 << "does not match number of eta bins." << endl;
0127 }
0128 }
0129
0130
0131 ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
0132
0133
0134 service_ = std::make_unique<MuonServiceProxy>(serviceParameters, consumesCollector());
0135
0136
0137 estimator_ = std::make_unique<Chi2MeasurementEstimator>(10000.);
0138
0139 produces<L2MuonTrajectorySeedCollection>();
0140 }
0141
0142 void L2MuonSeedGeneratorFromL1TkMu::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0143 edm::ParameterSetDescription desc;
0144 desc.add<edm::InputTag>("InputObjects", edm::InputTag("hltGmtStage2Digis"));
0145 desc.add<string>("Propagator", "");
0146 desc.add<double>("L1MinPt", -1.);
0147 desc.add<double>("L1MaxEta", 5.0);
0148 desc.add<double>("SetMinPtBarrelTo", 3.5);
0149 desc.add<double>("SetMinPtEndcapTo", 1.0);
0150 desc.add<double>("MinPL1Tk", 3.5);
0151 desc.add<double>("MinPtL1TkBarrel", 3.5);
0152 desc.addUntracked<bool>("UseOfflineSeed", false);
0153 desc.add<bool>("UseUnassociatedL1", true);
0154 desc.add<std::vector<double>>("MatchDR", {0.3});
0155 desc.add<std::vector<double>>("EtaMatchingBins", {0., 2.5});
0156 desc.addUntracked<edm::InputTag>("OfflineSeedLabel", edm::InputTag(""));
0157
0158 edm::ParameterSetDescription psd0;
0159 psd0.addUntracked<std::vector<std::string>>("Propagators", {"SteppingHelixPropagatorAny"});
0160 psd0.add<bool>("RPCLayers", true);
0161 psd0.addUntracked<bool>("UseMuonNavigation", true);
0162 desc.add<edm::ParameterSetDescription>("ServiceParameters", psd0);
0163 descriptions.add("L2MuonSeedGeneratorFromL1TkMu", desc);
0164 }
0165
0166 void L2MuonSeedGeneratorFromL1TkMu::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0167 const std::string metname = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1TkMu";
0168 MuonPatternRecoDumper debug;
0169
0170 auto output = std::make_unique<L2MuonTrajectorySeedCollection>();
0171
0172 auto const muColl = iEvent.getHandle(muCollToken_);
0173 LogDebug(metname) << "Number of muons " << muColl->size() << endl;
0174
0175 edm::Handle<edm::View<TrajectorySeed>> offlineSeedHandle;
0176 vector<int> offlineSeedMap;
0177 if (useOfflineSeed_) {
0178 iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
0179 offlineSeedMap = vector<int>(offlineSeedHandle->size(), 0);
0180 }
0181
0182 for (auto const &tkmu : *muColl) {
0183
0184 auto const &it = tkmu.trkPtr();
0185
0186
0187 auto p3 = it->momentum();
0188
0189 float tk_p = p3.mag();
0190 if (tk_p < minPL1Tk_)
0191 continue;
0192
0193 float tk_pt = p3.perp();
0194 float tk_eta = p3.eta();
0195 float tk_aeta = std::abs(tk_eta);
0196
0197 bool barrel = tk_aeta < etaBoundary_;
0198 if (barrel && tk_pt < minPtL1TkBarrel_)
0199 continue;
0200
0201 float tk_phi = p3.phi();
0202 float tk_q = it->rInv() > 0 ? 1. : -1.;
0203 float tk_z = it->POCA().z();
0204
0205 float dzCorrPhi = 1.;
0206 float deta = 0;
0207 float etaProp = tk_aeta;
0208
0209 if (barrel) {
0210 etaProp = etaBoundary_;
0211 deta = tk_z / distMB2_ / cosh(tk_aeta);
0212 } else {
0213 float delta = tk_z / distME2_;
0214 if (tk_eta > 0)
0215 delta *= -1;
0216 dzCorrPhi = 1. + delta;
0217
0218 float zOzs = tk_z / distME2_;
0219 if (tk_eta > 0)
0220 deta = zOzs / (1. - zOzs);
0221 else
0222 deta = zOzs / (1. + zOzs);
0223 deta = deta * tanh(tk_eta);
0224 }
0225 float resPhi = tk_phi - phiCorr0_ * tk_q * cosh(phiCorr1_) / cosh(etaProp) / tk_pt * dzCorrPhi - M_PI / phiCorr2_;
0226 resPhi = reco::reducePhiRange(resPhi);
0227
0228 float pt = tk_pt;
0229 float eta = tk_eta + deta;
0230 float theta = 2 * atan(exp(-eta));
0231 float phi = resPhi;
0232 int charge = it->rInv() > 0 ? 1 : -1;
0233
0234 if (pt < l1MinPt_ || std::abs(eta) > l1MaxEta_)
0235 continue;
0236
0237 LogDebug(metname) << "New L2 Muon Seed";
0238 LogDebug(metname) << "Pt = " << pt << " GeV/c";
0239 LogDebug(metname) << "eta = " << eta;
0240 LogDebug(metname) << "theta = " << theta << " rad";
0241 LogDebug(metname) << "phi = " << phi << " rad";
0242 LogDebug(metname) << "charge = " << charge;
0243 LogDebug(metname) << "In Barrel? = " << barrel;
0244
0245
0246 service_->update(iSetup);
0247
0248 const DetLayer *detLayer = nullptr;
0249 float radius = 0.;
0250
0251 CLHEP::Hep3Vector vec(0., 1., 0.);
0252 vec.setTheta(theta);
0253 vec.setPhi(phi);
0254
0255 DetId theid;
0256
0257 if (barrel) {
0258 LogDebug(metname) << "The seed is in the barrel";
0259
0260
0261 theid = DTChamberId(0, 2, 0);
0262 detLayer = service_->detLayerGeometry()->idToLayer(theid);
0263 LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0264
0265 const BoundSurface *sur = &(detLayer->surface());
0266 const BoundCylinder *bc = dynamic_cast<const BoundCylinder *>(sur);
0267
0268 radius = std::abs(bc->radius() / sin(theta));
0269
0270 LogDebug(metname) << "radius " << radius;
0271
0272 if (pt < minPtBarrel_)
0273 pt = minPtBarrel_;
0274 } else {
0275 LogDebug(metname) << "The seed is in the endcap";
0276
0277
0278 theid = theta < Geom::pi() / 2. ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
0279
0280 detLayer = service_->detLayerGeometry()->idToLayer(theid);
0281 LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0282
0283 radius = std::abs(detLayer->position().z() / cos(theta));
0284
0285 if (pt < minPtEndcap_)
0286 pt = minPtEndcap_;
0287 }
0288
0289 vec.setMag(radius);
0290
0291 GlobalPoint pos(vec.x(), vec.y(), vec.z());
0292
0293 GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));
0294
0295 GlobalTrajectoryParameters param(pos, mom, charge, &*service_->magneticField());
0296 AlgebraicSymMatrix55 mat;
0297
0298 mat[0][0] = (0.25 / pt) * (0.25 / pt);
0299 if (!barrel)
0300 mat[0][0] = (0.4 / pt) * (0.4 / pt);
0301
0302 mat[1][1] = 0.05 * 0.05;
0303 mat[2][2] = 0.2 * 0.2;
0304 mat[3][3] = 20. * 20.;
0305 mat[4][4] = 20. * 20.;
0306
0307 CurvilinearTrajectoryError error(mat);
0308
0309 const FreeTrajectoryState state(param, error);
0310
0311 LogDebug(metname) << "Free trajectory State from the parameters";
0312 LogDebug(metname) << debug.dumpFTS(state);
0313
0314
0315 TrajectoryStateOnSurface tsos = service_->propagator(propagatorName_)->propagate(state, detLayer->surface());
0316
0317 LogDebug(metname) << "State after the propagation on the layer";
0318 LogDebug(metname) << debug.dumpLayer(detLayer);
0319 LogDebug(metname) << debug.dumpTSOS(tsos);
0320
0321 double dRcone = matchingDR_[0];
0322 if (std::abs(eta) < etaBins_.back()) {
0323 std::vector<double>::iterator lowEdge = std::upper_bound(etaBins_.begin(), etaBins_.end(), std::abs(eta));
0324 dRcone = matchingDR_.at(lowEdge - etaBins_.begin() - 1);
0325 }
0326
0327 if (tsos.isValid()) {
0328 edm::OwnVector<TrackingRecHit> container;
0329
0330 if (useOfflineSeed_) {
0331
0332 std::vector<pair<const GeomDet *, TrajectoryStateOnSurface>> detsWithStates =
0333 detLayer->compatibleDets(tsos, *service_->propagator(propagatorName_), *estimator_);
0334
0335 if (detsWithStates.empty() && barrel) {
0336
0337 DetId fallback_id;
0338 theta < Geom::pi() / 2. ? fallback_id = CSCDetId(1, 2, 0, 0, 0) : fallback_id = CSCDetId(2, 2, 0, 0, 0);
0339 const DetLayer *ME2DetLayer = service_->detLayerGeometry()->idToLayer(fallback_id);
0340
0341 tsos = service_->propagator(propagatorName_)->propagate(state, ME2DetLayer->surface());
0342 detsWithStates = ME2DetLayer->compatibleDets(tsos, *service_->propagator(propagatorName_), *estimator_);
0343 }
0344
0345 if (!detsWithStates.empty()) {
0346 TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
0347 const GeomDet *newTSOSDet = detsWithStates.front().first;
0348
0349 LogDebug(metname) << "Most compatible det";
0350 LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
0351
0352 if (newTSOS.isValid()) {
0353 LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag()
0354 << ", phi=" << newTSOS.globalPosition().phi()
0355 << ", eta=" << newTSOS.globalPosition().eta() << ")";
0356 LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge() * newTSOS.globalMomentum().perp()
0357 << ", phi=" << newTSOS.globalMomentum().phi()
0358 << ", eta=" << newTSOS.globalMomentum().eta() << ")";
0359
0360 const TrajectorySeed *assoOffseed =
0361 associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, newTSOS, dRcone);
0362
0363 if (assoOffseed != nullptr) {
0364 PTrajectoryStateOnDet const &seedTSOS = assoOffseed->startingState();
0365 for (auto const &recHit : assoOffseed->recHits()) {
0366 container.push_back(recHit);
0367 }
0368 auto dummyRef = edm::Ref<MuonBxCollection>();
0369 output->emplace_back(L2MuonTrajectorySeed(seedTSOS, container, alongMomentum, dummyRef));
0370 } else {
0371 if (useUnassociatedL1_) {
0372
0373 PTrajectoryStateOnDet const &seedTSOS =
0374 trajectoryStateTransform::persistentState(newTSOS, newTSOSDet->geographicalId().rawId());
0375 auto dummyRef = edm::Ref<MuonBxCollection>();
0376 output->emplace_back(L2MuonTrajectorySeed(seedTSOS, container, alongMomentum, dummyRef));
0377 }
0378 }
0379 }
0380 }
0381 } else {
0382
0383 PTrajectoryStateOnDet const &seedTSOS = trajectoryStateTransform::persistentState(tsos, theid.rawId());
0384 auto dummyRef = edm::Ref<MuonBxCollection>();
0385 output->emplace_back(L2MuonTrajectorySeed(seedTSOS, container, alongMomentum, dummyRef));
0386 }
0387 }
0388 }
0389
0390 iEvent.put(std::move(output));
0391 }
0392
0393
0394 const TrajectorySeed *L2MuonSeedGeneratorFromL1TkMu::associateOfflineSeedToL1(
0395 edm::Handle<edm::View<TrajectorySeed>> &offseeds,
0396 std::vector<int> &offseedMap,
0397 TrajectoryStateOnSurface &newTsos,
0398 double dRcone) {
0399 if (dRcone < 0.)
0400 return nullptr;
0401
0402 const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1TkMu";
0403 MuonPatternRecoDumper debugtmp;
0404
0405 edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
0406 const TrajectorySeed *selOffseed = nullptr;
0407 double bestDr2 = 99999.;
0408 unsigned int nOffseed(0);
0409 int lastOffseed(-1);
0410
0411 for (offseed = offseeds->begin(); offseed != endOffseed; ++offseed, ++nOffseed) {
0412 if (offseedMap[nOffseed] != 0)
0413 continue;
0414 GlobalPoint glbPos = service_->trackingGeometry()
0415 ->idToDet(offseed->startingState().detId())
0416 ->surface()
0417 .toGlobal(offseed->startingState().parameters().position());
0418 GlobalVector glbMom = service_->trackingGeometry()
0419 ->idToDet(offseed->startingState().detId())
0420 ->surface()
0421 .toGlobal(offseed->startingState().parameters().momentum());
0422
0423
0424 double preDr2 = deltaR2(newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi());
0425 if (preDr2 > 1.0)
0426 continue;
0427
0428 const FreeTrajectoryState offseedFTS(
0429 glbPos, glbMom, offseed->startingState().parameters().charge(), &*service_->magneticField());
0430 TrajectoryStateOnSurface offseedTsos =
0431 service_->propagator(propagatorName_)->propagate(offseedFTS, newTsos.surface());
0432 LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
0433 LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
0434 LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" << offseedFTS.position().phi()
0435 << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
0436 LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge() * offseedFTS.momentum().perp()
0437 << ", phi=" << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")"
0438 << std::endl
0439 << std::endl;
0440
0441 if (offseedTsos.isValid()) {
0442 LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
0443 LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag()
0444 << ", phi=" << offseedTsos.globalPosition().phi()
0445 << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
0446 LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge() * offseedTsos.globalMomentum().perp()
0447 << ", phi=" << offseedTsos.globalMomentum().phi()
0448 << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl
0449 << std::endl;
0450 double newDr2 = deltaR2(newTsos.globalPosition().eta(),
0451 newTsos.globalPosition().phi(),
0452 offseedTsos.globalPosition().eta(),
0453 offseedTsos.globalPosition().phi());
0454 LogDebug(metlabel) << " -- DR = " << newDr2 << std::endl;
0455 if (newDr2 < bestDr2 && newDr2 < dRcone * dRcone) {
0456 LogDebug(metlabel) << " --> OK! " << newDr2 << std::endl << std::endl;
0457 selOffseed = &*offseed;
0458 bestDr2 = newDr2;
0459 offseedMap[nOffseed] = 1;
0460 if (lastOffseed > -1)
0461 offseedMap[lastOffseed] = 0;
0462 lastOffseed = nOffseed;
0463 } else {
0464 LogDebug(metlabel) << " --> Rejected. " << newDr2 << std::endl << std::endl;
0465 }
0466 } else {
0467 LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
0468 }
0469 }
0470
0471 return selOffseed;
0472 }
0473
0474
0475 DEFINE_FWK_MODULE(L2MuonSeedGeneratorFromL1TkMu);