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