TkTrackingRegions/plugins/L1MuonSeededTrackingRegionsProducer.h

0001 #ifndef RecoTracker_TkTrackingRegions_L1MuonSeededTrackingRegionsProducer_h
0002 #define RecoTracker_TkTrackingRegions_L1MuonSeededTrackingRegionsProducer_h
0003
0004 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegionProducer.h"
0005 #include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
0006 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Framework/interface/ConsumesCollector.h"
0013 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0014 #include "DataFormats/Common/interface/Handle.h"
0015 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0016 #include "DataFormats/VertexReco/interface/Vertex.h"
0017 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0018 #include "DataFormats/L1Trigger/interface/Muon.h"
0019 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0020 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0021 #include "MagneticField/Engine/interface/MagneticField.h"
0022 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0023 #include "RecoTracker/Record/interface/TrackerMultipleScatteringRecord.h"
0024 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisationMaker.h"
0025 #include "CLHEP/Vector/ThreeVector.h"
0026
0027 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0028
0029 /** class L1MuonSeededTrackingRegionsProducer
0030  *
0031  * eta-phi TrackingRegions producer in directions defined by L1 muon objects of interest
0032  * from a collection defined by the "input" parameter.
0033  *
0034  * Four operational modes are supported ("mode" parameter):
0035  *
0036  *   BeamSpotFixed:
0037  *     origin is defined by the beam spot
0038  *     z-half-length is defined by a fixed zErrorBeamSpot parameter
0039  *   BeamSpotSigma:
0040  *     origin is defined by the beam spot
0041  *     z-half-length is defined by nSigmaZBeamSpot * beamSpot.sigmaZ
0042  *   VerticesFixed:
0043  *     origins are defined by vertices from VertexCollection (use maximum MaxNVertices of them)
0044  *     z-half-length is defined by a fixed zErrorVetex parameter
0045  *   VerticesSigma:
0046  *     origins are defined by vertices from VertexCollection (use maximum MaxNVertices of them)
0047  *     z-half-length is defined by nSigmaZVertex * vetex.zError
0048  *
0049  *   If, while using one of the "Vertices" modes, there's no vertices in an event, we fall back into
0050  *   either BeamSpotSigma or BeamSpotFixed mode, depending on the positiveness of nSigmaZBeamSpot.
0051  *
0052  *   \author M. Oh.
0053  *       based on RecoTracker/TkTrackingRegions/plugins/CandidateSeededTrackingRegionsProducer.h
0054  */
0055 class L1MuonSeededTrackingRegionsProducer : public TrackingRegionProducer {
0056 public:
0057   typedef enum { BEAM_SPOT_FIXED, BEAM_SPOT_SIGMA, VERTICES_FIXED, VERTICES_SIGMA } Mode;
0058
0059   explicit L1MuonSeededTrackingRegionsProducer(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC)
0060       : token_field(iC.esConsumes()),
0061         l1MinPt_(iConfig.getParameter<double>("L1MinPt")),
0062         l1MaxEta_(iConfig.getParameter<double>("L1MaxEta")),
0063         l1MinQuality_(iConfig.getParameter<unsigned int>("L1MinQuality")),
0064         minPtBarrel_(iConfig.getParameter<double>("SetMinPtBarrelTo")),
0065         minPtEndcap_(iConfig.getParameter<double>("SetMinPtEndcapTo")),
0066         centralBxOnly_(iConfig.getParameter<bool>("CentralBxOnly")),
0067         propagatorName_(iConfig.getParameter<std::string>("Propagator")) {
0068     edm::ParameterSet regPSet = iConfig.getParameter<edm::ParameterSet>("RegionPSet");
0069
0070     // operation mode
0071     std::string modeString = regPSet.getParameter<std::string>("mode");
0072     if (modeString == "BeamSpotFixed")
0073       m_mode = BEAM_SPOT_FIXED;
0074     else if (modeString == "BeamSpotSigma")
0075       m_mode = BEAM_SPOT_SIGMA;
0076     else if (modeString == "VerticesFixed")
0077       m_mode = VERTICES_FIXED;
0078     else if (modeString == "VerticesSigma")
0079       m_mode = VERTICES_SIGMA;
0080     else
0081       edm::LogError("L1MuonSeededTrackingRegionsProducer") << "Unknown mode string: " << modeString;
0082
0083     // basic inputs
0084     token_input = iC.consumes<l1t::MuonBxCollection>(regPSet.getParameter<edm::InputTag>("input"));
0085     m_maxNRegions = regPSet.getParameter<int>("maxNRegions");
0086     token_beamSpot = iC.consumes<reco::BeamSpot>(regPSet.getParameter<edm::InputTag>("beamSpot"));
0087     m_maxNVertices = 1;
0088     if (m_mode == VERTICES_FIXED || m_mode == VERTICES_SIGMA) {
0089       token_vertex = iC.consumes<reco::VertexCollection>(regPSet.getParameter<edm::InputTag>("vertexCollection"));
0090       m_maxNVertices = regPSet.getParameter<int>("maxNVertices");
0091     }
0092
0093     // RectangularEtaPhiTrackingRegion parameters:
0094     m_ptMin = regPSet.getParameter<double>("ptMin");
0095     m_originRadius = regPSet.getParameter<double>("originRadius");
0096     m_zErrorBeamSpot = regPSet.getParameter<double>("zErrorBeamSpot");
0097     m_ptRanges = regPSet.getParameter<std::vector<double>>("ptRanges");
0098     if (m_ptRanges.size() < 2) {
0099       edm::LogError("L1MuonSeededTrackingRegionsProducer") << "Size of ptRanges does not be less than 2" << std::endl;
0100     }
0101     m_deltaEtas = regPSet.getParameter<std::vector<double>>("deltaEtas");
0102     if (m_deltaEtas.size() != m_ptRanges.size() - 1) {
0103       edm::LogError("L1MuonSeededTrackingRegionsProducer")
0104           << "Size of deltaEtas does not match number of pt bins." << std::endl;
0105     }
0106     m_deltaPhis = regPSet.getParameter<std::vector<double>>("deltaPhis");
0107     if (m_deltaPhis.size() != m_ptRanges.size() - 1) {
0108       edm::LogError("L1MuonSeededTrackingRegionsProducer")
0109           << "Size of deltaPhis does not match number of pt bins." << std::endl;
0110     }
0111
0112     m_precise = regPSet.getParameter<bool>("precise");
0113     m_whereToUseMeasurementTracker = RectangularEtaPhiTrackingRegion::stringToUseMeasurementTracker(
0114         regPSet.getParameter<std::string>("whereToUseMeasurementTracker"));
0115     if (m_whereToUseMeasurementTracker != RectangularEtaPhiTrackingRegion::UseMeasurementTracker::kNever) {
0116       token_measurementTracker =
0117           iC.consumes<MeasurementTrackerEvent>(regPSet.getParameter<edm::InputTag>("measurementTrackerName"));
0118     }
0119
0120     m_searchOpt = regPSet.getParameter<bool>("searchOpt");
0121
0122     // mode-dependent z-halflength of tracking regions
0123     if (m_mode == VERTICES_SIGMA)
0124       m_nSigmaZVertex = regPSet.getParameter<double>("nSigmaZVertex");
0125     if (m_mode == VERTICES_FIXED)
0126       m_zErrorVetex = regPSet.getParameter<double>("zErrorVetex");
0127     m_nSigmaZBeamSpot = -1.;
0128     if (m_mode == BEAM_SPOT_SIGMA) {
0129       m_nSigmaZBeamSpot = regPSet.getParameter<double>("nSigmaZBeamSpot");
0130       if (m_nSigmaZBeamSpot < 0.)
0131         edm::LogError("L1MuonSeededTrackingRegionsProducer")
0132             << "nSigmaZBeamSpot must be positive for BeamSpotSigma mode!";
0133     }
0134     if (m_precise) {
0135       token_msmaker = iC.esConsumes();
0136     }
0137
0138     // MuonServiceProxy
0139     edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
0140     service_ = std::make_unique<MuonServiceProxy>(serviceParameters, std::move(iC));
0141   }
0142
0143   ~L1MuonSeededTrackingRegionsProducer() override = default;
0144
0145   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0146     edm::ParameterSetDescription desc;
0147
0148     // L1 muon selection parameters
0149     desc.add<std::string>("Propagator", "");
0150     desc.add<double>("L1MinPt", -1.);
0151     desc.add<double>("L1MaxEta", 5.0);
0152     desc.add<unsigned int>("L1MinQuality", 0);
0153     desc.add<double>("SetMinPtBarrelTo", 3.5);
0154     desc.add<double>("SetMinPtEndcapTo", 1.0);
0155     desc.add<bool>("CentralBxOnly", true);
0156
0157     // Tracking region parameters
0158     edm::ParameterSetDescription descRegion;
0159     descRegion.add<std::string>("mode", "BeamSpotSigma");
0160     descRegion.add<edm::InputTag>("input", edm::InputTag(""));
0161     descRegion.add<int>("maxNRegions", 10);
0162     descRegion.add<edm::InputTag>("beamSpot", edm::InputTag("hltOnlineBeamSpot"));
0163     descRegion.add<edm::InputTag>("vertexCollection", edm::InputTag("notUsed"));
0164     descRegion.add<int>("maxNVertices", 1);
0165     descRegion.add<double>("ptMin", 0.0);
0166     descRegion.add<double>("originRadius", 0.2);
0167     descRegion.add<double>("zErrorBeamSpot", 24.2);
0168     descRegion.add<std::vector<double>>("ptRanges", {0., 1.e9});
0169     descRegion.add<std::vector<double>>("deltaEtas", {0.35});
0170     descRegion.add<std::vector<double>>("deltaPhis", {0.2});
0171     descRegion.add<bool>("precise", true);
0172     descRegion.add<double>("nSigmaZVertex", 3.);
0173     descRegion.add<double>("zErrorVetex", 0.2);
0174     descRegion.add<double>("nSigmaZBeamSpot", 4.);
0175     descRegion.add<std::string>("whereToUseMeasurementTracker", "Never");
0176     descRegion.add<edm::InputTag>("measurementTrackerName", edm::InputTag(""));
0177     descRegion.add<bool>("searchOpt", false);
0178     desc.add<edm::ParameterSetDescription>("RegionPSet", descRegion);
0179
0180     // MuonServiceProxy for the propagation
0181     edm::ParameterSetDescription psd0;
0182     psd0.addUntracked<std::vector<std::string>>("Propagators", {"SteppingHelixPropagatorAny"});
0183     psd0.add<bool>("RPCLayers", false);
0184     psd0.addUntracked<bool>("UseMuonNavigation", false);
0185     desc.add<edm::ParameterSetDescription>("ServiceParameters", psd0);
0186
0187     descriptions.add("hltIterL3MuonPixelTracksTrackingRegions", desc);
0188   }
0189
0190   std::vector<std::unique_ptr<TrackingRegion>> regions(const edm::Event& iEvent,
0191                                                        const edm::EventSetup& iSetup) const override {
0192     std::vector<std::unique_ptr<TrackingRegion>> result;
0193
0194     // pick up the candidate objects of interest
0195     edm::Handle<l1t::MuonBxCollection> muColl;
0196     iEvent.getByToken(token_input, muColl);
0197     if (muColl->size() == 0)
0198       return result;
0199
0200     // always need the beam spot (as a fall back strategy for vertex modes)
0201     edm::Handle<reco::BeamSpot> bs;
0202     iEvent.getByToken(token_beamSpot, bs);
0203     if (!bs.isValid())
0204       return result;
0205
0206     // this is a default origin for all modes
0207     GlobalPoint default_origin(bs->x0(), bs->y0(), bs->z0());
0208
0209     // vector of origin & halfLength pairs:
0210     std::vector<std::pair<GlobalPoint, float>> origins;
0211
0212     // fill the origins and halfLengths depending on the mode
0213     if (m_mode == BEAM_SPOT_FIXED || m_mode == BEAM_SPOT_SIGMA) {
0214       origins.push_back(std::make_pair(
0215           default_origin, (m_mode == BEAM_SPOT_FIXED) ? m_zErrorBeamSpot : m_nSigmaZBeamSpot * bs->sigmaZ()));
0216     } else if (m_mode == VERTICES_FIXED || m_mode == VERTICES_SIGMA) {
0217       edm::Handle<reco::VertexCollection> vertices;
0218       iEvent.getByToken(token_vertex, vertices);
0219       int n_vert = 0;
0220       for (reco::VertexCollection::const_iterator v = vertices->begin();
0221            v != vertices->end() && n_vert < m_maxNVertices;
0222            ++v) {
0223         if (v->isFake() || !v->isValid())
0224           continue;
0225
0226         origins.push_back(std::make_pair(GlobalPoint(v->x(), v->y(), v->z()),
0227                                          (m_mode == VERTICES_FIXED) ? m_zErrorVetex : m_nSigmaZVertex * v->zError()));
0228         ++n_vert;
0229       }
0230       // no-vertex fall-back case:
0231       if (origins.empty()) {
0232         origins.push_back(std::make_pair(
0233             default_origin, (m_nSigmaZBeamSpot > 0.) ? m_nSigmaZBeamSpot * bs->z0Error() : m_zErrorBeamSpot));
0234       }
0235     }
0236
0237     const MeasurementTrackerEvent* measurementTracker = nullptr;
0238     if (!token_measurementTracker.isUninitialized()) {
0239       measurementTracker = &iEvent.get(token_measurementTracker);
0240     }
0241
0242     const auto& field = iSetup.getData(token_field);
0243     const MultipleScatteringParametrisationMaker* msmaker = nullptr;
0244     if (m_precise) {
0245       msmaker = &iSetup.getData(token_msmaker);
0246     }
0247
0248     // create tracking regions (maximum MaxNRegions of them) in directions of the
0249     // objects of interest (we expect that the collection was sorted in decreasing pt order)
0250     int n_regions = 0;
0251     for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX() && n_regions < m_maxNRegions; ++ibx) {
0252       if (centralBxOnly_ && (ibx != 0))
0253         continue;
0254
0255       for (auto it = muColl->begin(ibx); it != muColl->end(ibx) && n_regions < m_maxNRegions; it++) {
0256         unsigned int quality = it->hwQual();
0257         if (quality <= l1MinQuality_)
0258           continue;
0259
0260         float pt = it->pt();
0261         float eta = it->eta();
0262         if (pt < l1MinPt_ || std::abs(eta) > l1MaxEta_)
0263           continue;
0264
0265         float theta = 2 * atan(exp(-eta));
0266         float phi = it->phi();
0267
0268         int valid_charge = it->hwChargeValid();
0269         int charge = it->charge();
0270         if (!valid_charge)
0271           charge = 0;
0272
0273         int link = l1MuonTF_link_EMTFP_i_ + (int)(it->tfMuonIndex() / 3.);
0274         bool barrel = true;
0275         if ((link >= l1MuonTF_link_EMTFP_i_ && link <= l1MuonTF_link_EMTFP_f_) ||
0276             (link >= l1MuonTF_link_EMTFN_i_ && link <= l1MuonTF_link_EMTFN_f_))
0277           barrel = false;
0278
0279         // propagate L1 FTS to BS
0280         service_->update(iSetup);
0281         const DetLayer* detLayer = nullptr;
0282         float radius = 0.;
0283
0284         CLHEP::Hep3Vector vec(0., 1., 0.);
0285         vec.setTheta(theta);
0286         vec.setPhi(phi);
0287
0288         DetId theid;
0289         // Get the det layer on which the state should be put
0290         if (barrel) {
0291           // MB2
0292           theid = DTChamberId(0, 2, 0);
0293           detLayer = service_->detLayerGeometry()->idToLayer(theid);
0294
0295           const BoundSurface* sur = &(detLayer->surface());
0296           const BoundCylinder* bc = dynamic_cast<const BoundCylinder*>(sur);
0297
0298           radius = std::abs(bc->radius() / sin(theta));
0299
0300           if (pt < minPtBarrel_)
0301             pt = minPtBarrel_;
0302         } else {
0303           // ME2
0304           theid = theta < M_PI / 2. ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
0305
0306           detLayer = service_->detLayerGeometry()->idToLayer(theid);
0307
0308           radius = std::abs(detLayer->position().z() / cos(theta));
0309
0310           if (pt < minPtEndcap_)
0311             pt = minPtEndcap_;
0312         }
0313         vec.setMag(radius);
0314         GlobalPoint pos(vec.x(), vec.y(), vec.z());
0315         GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));
0316         GlobalTrajectoryParameters param(pos, mom, charge, &*service_->magneticField());
0317
0318         AlgebraicSymMatrix55 mat;
0319         mat[0][0] = (sigma_qbpt_barrel_ / pt) * (sigma_qbpt_barrel_ / pt);  // sigma^2(charge/abs_momentum)
0320         if (!barrel)
0321           mat[0][0] = (sigma_qbpt_endcap_ / pt) * (sigma_qbpt_endcap_ / pt);
0322
0323         //Assign q/pt = 0 +- 1/pt if charge has been declared invalid
0324         if (!valid_charge)
0325           mat[0][0] = (sigma_qbpt_invalid_charge_ / pt) * (sigma_qbpt_invalid_charge_ / pt);
0326
0327         mat[1][1] = sigma_lambda_ * sigma_lambda_;  // sigma^2(lambda)
0328         mat[2][2] = sigma_phi_ * sigma_phi_;        // sigma^2(phi)
0329         mat[3][3] = sigma_x_ * sigma_x_;            // sigma^2(x_transverse))
0330         mat[4][4] = sigma_y_ * sigma_y_;            // sigma^2(y_transverse))
0331
0332         CurvilinearTrajectoryError error(mat);
0333
0334         const FreeTrajectoryState state(param, error);
0335
0336         FreeTrajectoryState state_bs = service_->propagator(propagatorName_)->propagate(state, *bs.product());
0337
0338         GlobalVector direction(state_bs.momentum().x(), state_bs.momentum().y(), state_bs.momentum().z());
0339
0340         // set deltaEta and deltaPhi from L1 muon pt
0341         auto deltaEta = m_deltaEtas.at(0);
0342         auto deltaPhi = m_deltaPhis.at(0);
0343         if (it->pt() < m_ptRanges.back()) {
0344           auto lowEdge = std::upper_bound(m_ptRanges.begin(), m_ptRanges.end(), it->pt());
0345           deltaEta = m_deltaEtas.at(lowEdge - m_ptRanges.begin() - 1);
0346           deltaPhi = m_deltaPhis.at(lowEdge - m_ptRanges.begin() - 1);
0347         }
0348
0349         for (size_t j = 0; j < origins.size() && n_regions < m_maxNRegions; ++j) {
0350           result.push_back(std::make_unique<RectangularEtaPhiTrackingRegion>(direction,
0351                                                                              origins[j].first,
0352                                                                              m_ptMin,
0353                                                                              m_originRadius,
0354                                                                              origins[j].second,
0355                                                                              deltaEta,
0356                                                                              deltaPhi,
0357                                                                              field,
0358                                                                              msmaker,
0359                                                                              m_precise,
0360                                                                              m_whereToUseMeasurementTracker,
0361                                                                              measurementTracker,
0362                                                                              m_searchOpt));
0363           ++n_regions;
0364         }
0365       }
0366     }
0367     edm::LogInfo("L1MuonSeededTrackingRegionsProducer") << "produced " << n_regions << " regions";
0368
0369     return result;
0370   }
0371
0372 private:
0373   Mode m_mode;
0374
0375   int m_maxNRegions;
0376   edm::EDGetTokenT<reco::VertexCollection> token_vertex;
0377   edm::EDGetTokenT<reco::BeamSpot> token_beamSpot;
0378   edm::EDGetTokenT<l1t::MuonBxCollection> token_input;
0379   int m_maxNVertices;
0380
0381   float m_ptMin;
0382   float m_originRadius;
0383   float m_zErrorBeamSpot;
0384   std::vector<double> m_ptRanges;
0385   std::vector<double> m_deltaEtas;
0386   std::vector<double> m_deltaPhis;
0387   bool m_precise;
0388   edm::EDGetTokenT<MeasurementTrackerEvent> token_measurementTracker;
0389   RectangularEtaPhiTrackingRegion::UseMeasurementTracker m_whereToUseMeasurementTracker;
0390   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> token_field;
0391   edm::ESGetToken<MultipleScatteringParametrisationMaker, TrackerMultipleScatteringRecord> token_msmaker;
0392   bool m_searchOpt;
0393
0394   float m_nSigmaZVertex;
0395   float m_zErrorVetex;
0396   float m_nSigmaZBeamSpot;
0397
0398   const double l1MinPt_;
0399   const double l1MaxEta_;
0400   const unsigned l1MinQuality_;
0401   const double minPtBarrel_;
0402   const double minPtEndcap_;
0403   const bool centralBxOnly_;
0404
0405   const std::string propagatorName_;
0406   std::unique_ptr<MuonServiceProxy> service_;
0407
0408   // link number indices of the optical fibres that connect the uGMT with the track finders
0409   // EMTF+ : 36-41, OMTF+ : 42-47, BMTF : 48-59, OMTF- : 60-65, EMTF- : 66-71
0410   static constexpr int l1MuonTF_link_EMTFP_i_{36};
0411   static constexpr int l1MuonTF_link_EMTFP_f_{41};
0412   static constexpr int l1MuonTF_link_EMTFN_i_{66};
0413   static constexpr int l1MuonTF_link_EMTFN_f_{71};
0414
0415   // fixed error matrix parameters for L1 muon FTS
0416   static constexpr double sigma_qbpt_barrel_{0.25};
0417   static constexpr double sigma_qbpt_endcap_{0.4};
0418   static constexpr double sigma_qbpt_invalid_charge_{1.0};
0419   static constexpr double sigma_lambda_{0.05};
0420   static constexpr double sigma_phi_{0.2};
0421   static constexpr double sigma_x_{20.0};
0422   static constexpr double sigma_y_{20.0};
0423 };
0424
0425 #endif