Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-04-08 03:34:54 

0001 /**
0002   \class    pat::OutsideInMuonSeeder MuonReSeeder.h "MuonAnalysis/MuonAssociators/interface/MuonReSeeder.h"
0003   \brief    Matcher of reconstructed objects to other reconstructed objects using the tracks inside them
0004 
0005   \author   Giovanni Petrucciani
0006 */
0007 
0008 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0009 #include "DataFormats/Math/interface/deltaR.h"
0010 #include "DataFormats/MuonReco/interface/Muon.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0018 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0019 #include "MagneticField/Engine/interface/MagneticField.h"
0020 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0021 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0022 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0023 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0024 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0025 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
0026 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0027 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
0028 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0029 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0030 #include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
0031 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0032 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
0033 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0034 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0035 
0036 class OutsideInMuonSeeder final : public edm::stream::EDProducer<> {
0037 public:
0038   explicit OutsideInMuonSeeder(const edm::ParameterSet &iConfig);
0039   ~OutsideInMuonSeeder() override = default;
0040 
0041   void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override;
0042 
0043   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0044 
0045 private:
0046   /// Labels for input collections
0047   edm::EDGetTokenT<edm::View<reco::Muon>> src_;
0048 
0049   /// Muon selection
0050   StringCutObjectSelector<reco::Muon> selector_;
0051 
0052   /// How many layers to try
0053   const int layersToTry_;
0054 
0055   /// How many hits to try on same layer
0056   const int hitsToTry_;
0057 
0058   /// Do inside-out
0059   const bool fromVertex_;
0060 
0061   /// How much to rescale errors from STA
0062   const double errorRescaling_;
0063 
0064   const edm::ESGetToken<Propagator, TrackingComponentsRecord> trackerPropagatorToken_;
0065   const edm::ESGetToken<Propagator, TrackingComponentsRecord> muonPropagatorToken_;
0066   const edm::EDGetTokenT<MeasurementTrackerEvent> measurementTrackerTag_;
0067   const edm::ESGetToken<Chi2MeasurementEstimatorBase, TrackingComponentsRecord> estimatorToken_;
0068   const edm::ESGetToken<TrajectoryStateUpdator, TrackingComponentsRecord> updatorToken_;
0069   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magfieldToken_;
0070   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> geometryToken_;
0071   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeometryToken_;
0072 
0073   float const minEtaForTEC_;
0074   float const maxEtaForTOB_;
0075 
0076   const MagneticField *magfield_;
0077   const Propagator *muonPropagator_;
0078   const GlobalTrackingGeometry *geometry_;
0079   const Chi2MeasurementEstimatorBase *estimator_;
0080   const TrajectoryStateUpdator *updator_;
0081 
0082   /// Dump deug information
0083   const bool debug_;
0084 
0085   int doLayer(const GeometricSearchDet &layer,
0086               const TrajectoryStateOnSurface &state,
0087               std::vector<TrajectorySeed> &out,
0088               const Propagator &muon_propagator,
0089               const Propagator &tracker_propagator,
0090               const MeasurementTrackerEvent &mte) const;
0091   void doDebug(const reco::Track &tk) const;
0092 };
0093 
0094 OutsideInMuonSeeder::OutsideInMuonSeeder(const edm::ParameterSet &iConfig)
0095     : src_(consumes<edm::View<reco::Muon>>(iConfig.getParameter<edm::InputTag>("src"))),
0096       selector_(iConfig.existsAs<std::string>("cut") ? iConfig.getParameter<std::string>("cut") : "", true),
0097       layersToTry_(iConfig.getParameter<int32_t>("layersToTry")),
0098       hitsToTry_(iConfig.getParameter<int32_t>("hitsToTry")),
0099       fromVertex_(iConfig.getParameter<bool>("fromVertex")),
0100       errorRescaling_(iConfig.getParameter<double>("errorRescaleFactor")),
0101       trackerPropagatorToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("trackerPropagator")))),
0102       muonPropagatorToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("muonPropagator")))),
0103       measurementTrackerTag_(
0104           consumes<MeasurementTrackerEvent>(iConfig.getParameter<edm::InputTag>("measurementTkEvent"))),
0105       estimatorToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("hitCollector")))),
0106       updatorToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("updatorLabel")))),
0107       magfieldToken_(esConsumes()),
0108       geometryToken_(esConsumes()),
0109       tkGeometryToken_(esConsumes()),
0110       minEtaForTEC_(iConfig.getParameter<double>("minEtaForTEC")),
0111       maxEtaForTOB_(iConfig.getParameter<double>("maxEtaForTOB")),
0112       debug_(iConfig.getUntrackedParameter<bool>("debug", false)) {
0113   produces<std::vector<TrajectorySeed>>();
0114 }
0115 
0116 void OutsideInMuonSeeder::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0117   using namespace edm;
0118   using namespace std;
0119 
0120   magfield_ = &iSetup.getData(magfieldToken_);
0121   auto const &trackerPropagator = iSetup.getData(trackerPropagatorToken_);
0122   muonPropagator_ = &iSetup.getData(muonPropagatorToken_);
0123   geometry_ = &iSetup.getData(geometryToken_);
0124   estimator_ = &iSetup.getData(estimatorToken_);
0125   updator_ = &iSetup.getData(updatorToken_);
0126 
0127   Handle<MeasurementTrackerEvent> measurementTracker;
0128   iEvent.getByToken(measurementTrackerTag_, measurementTracker);
0129 
0130   const auto &tmpTkGeometry = iSetup.getData(tkGeometryToken_);
0131 
0132   Handle<View<reco::Muon>> src;
0133   iEvent.getByToken(src_, src);
0134 
0135   auto out = std::make_unique<std::vector<TrajectorySeed>>();
0136 
0137   for (auto const &mu : *src) {
0138     if (mu.outerTrack().isNull() || !selector_(mu))
0139       continue;
0140     if (debug_ && mu.innerTrack().isNonnull())
0141       doDebug(*mu.innerTrack());
0142 
0143     // Better clone here and not directly into doLayer to avoid
0144     // useless clone/destroy operations to set, in the end, the
0145     // very same direction every single time.
0146     std::unique_ptr<Propagator> pmuon_cloned =
0147         SetPropagationDirection(*muonPropagator_, fromVertex_ ? alongMomentum : oppositeToMomentum);
0148     std::unique_ptr<Propagator> ptracker_cloned = SetPropagationDirection(trackerPropagator, alongMomentum);
0149 
0150     int sizeBefore = out->size();
0151     if (debug_)
0152       LogDebug("OutsideInMuonSeeder") << "\n\n\nSeeding for muon of pt " << mu.pt() << ", eta " << mu.eta() << ", phi "
0153                                       << mu.phi() << std::endl;
0154     const reco::Track &tk = *mu.outerTrack();
0155 
0156     TrajectoryStateOnSurface state =
0157         fromVertex_ ? TrajectoryStateOnSurface(trajectoryStateTransform::initialFreeState(tk, magfield_))
0158                     : trajectoryStateTransform::innerStateOnSurface(tk, *geometry_, magfield_);
0159 
0160     if (std::abs(tk.eta()) < maxEtaForTOB_) {
0161       std::vector<BarrelDetLayer const *> const &tob = measurementTracker->geometricSearchTracker()->tobLayers();
0162       int found = 0;
0163       int iLayer = tob.size();
0164       if (iLayer == 0)
0165         LogError("OutsideInMuonSeeder") << "TOB has no layers.";
0166 
0167       for (auto it = tob.rbegin(), ed = tob.rend(); it != ed; ++it, --iLayer) {
0168         if (debug_)
0169           LogDebug("OutsideInMuonSeeder") << "\n ==== Trying TOB " << iLayer << " ====" << std::endl;
0170         if (doLayer(**it, state, *out, *(pmuon_cloned.get()), *(ptracker_cloned.get()), *measurementTracker)) {
0171           if (++found == layersToTry_)
0172             break;
0173         }
0174       }
0175     }
0176     if (tk.eta() > minEtaForTEC_) {
0177       const auto &forwLayers = tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC)
0178                                    ? measurementTracker->geometricSearchTracker()->posTidLayers()
0179                                    : measurementTracker->geometricSearchTracker()->posTecLayers();
0180       if (tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC)) {
0181         LogDebug("OutsideInMuonSeeder") << "\n We are using the Phase2 Outer Tracker (defined as a TID+). ";
0182       }
0183       LogTrace("OutsideInMuonSeeder") << "\n ==== TEC+ tot layers " << forwLayers.size() << " ====" << std::endl;
0184       int found = 0;
0185       int iLayer = forwLayers.size();
0186       if (iLayer == 0)
0187         LogError("OutsideInMuonSeeder") << "TEC+ has no layers.";
0188 
0189       if (debug_)
0190         LogDebug("OutsideInMuonSeeder") << "\n ==== Tot layers " << forwLayers.size() << " ====" << std::endl;
0191       for (auto it = forwLayers.rbegin(), ed = forwLayers.rend(); it != ed; ++it, --iLayer) {
0192         if (debug_)
0193           LogDebug("OutsideInMuonSeeder") << "\n ==== Trying Forward Layer +" << +iLayer << " ====" << std::endl;
0194         if (doLayer(**it, state, *out, *(pmuon_cloned.get()), *(ptracker_cloned.get()), *measurementTracker)) {
0195           if (++found == layersToTry_)
0196             break;
0197         }
0198       }
0199     }
0200     if (tk.eta() < -minEtaForTEC_) {
0201       const auto &forwLayers = tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC)
0202                                    ? measurementTracker->geometricSearchTracker()->negTidLayers()
0203                                    : measurementTracker->geometricSearchTracker()->negTecLayers();
0204       if (tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC)) {
0205         LogDebug("OutsideInMuonSeeder") << "\n We are using the Phase2 Outer Tracker (defined as a TID-). ";
0206       }
0207       LogTrace("OutsideInMuonSeeder") << "\n ==== TEC- tot layers " << forwLayers.size() << " ====" << std::endl;
0208       int found = 0;
0209       int iLayer = forwLayers.size();
0210       if (iLayer == 0)
0211         LogError("OutsideInMuonSeeder") << "TEC- has no layers.";
0212 
0213       if (debug_)
0214         LogDebug("OutsideInMuonSeeder") << "\n ==== Tot layers " << forwLayers.size() << " ====" << std::endl;
0215       for (auto it = forwLayers.rbegin(), ed = forwLayers.rend(); it != ed; ++it, --iLayer) {
0216         if (debug_)
0217           LogDebug("OutsideInMuonSeeder") << "\n ==== Trying Forward Layer -" << -iLayer << " ====" << std::endl;
0218         if (doLayer(**it, state, *out, *(pmuon_cloned.get()), *(ptracker_cloned.get()), *measurementTracker)) {
0219           if (++found == layersToTry_)
0220             break;
0221         }
0222       }
0223     }
0224     if (debug_)
0225       LogDebug("OutsideInMuonSeeder") << "Outcome of seeding for muon of pt " << mu.pt() << ", eta " << mu.eta()
0226                                       << ", phi " << mu.phi() << ": found " << (out->size() - sizeBefore) << " seeds."
0227                                       << std::endl;
0228   }
0229 
0230   iEvent.put(std::move(out));
0231 }
0232 
0233 int OutsideInMuonSeeder::doLayer(const GeometricSearchDet &layer,
0234                                  const TrajectoryStateOnSurface &state,
0235                                  std::vector<TrajectorySeed> &out,
0236                                  const Propagator &muon_propagator,
0237                                  const Propagator &tracker_propagator,
0238                                  const MeasurementTrackerEvent &measurementTracker) const {
0239   TrajectoryStateOnSurface onLayer(state);
0240   onLayer.rescaleError(errorRescaling_);
0241   std::vector<GeometricSearchDet::DetWithState> dets;
0242   layer.compatibleDetsV(onLayer, muon_propagator, *estimator_, dets);
0243 
0244   if (debug_) {
0245     LogDebug("OutsideInMuonSeeder") << "Query on layer around x = " << onLayer.globalPosition()
0246                                     << " with local pos error " << sqrt(onLayer.localError().positionError().xx())
0247                                     << " ,  " << sqrt(onLayer.localError().positionError().yy()) << " ,  "
0248                                     << " returned " << dets.size() << " compatible detectors" << std::endl;
0249   }
0250 
0251   std::vector<TrajectoryMeasurement> meas;
0252   for (std::vector<GeometricSearchDet::DetWithState>::const_iterator it = dets.begin(), ed = dets.end(); it != ed;
0253        ++it) {
0254     MeasurementDetWithData det = measurementTracker.idToDet(it->first->geographicalId());
0255     if (det.isNull()) {
0256       std::cerr << "BOGUS detid " << it->first->geographicalId().rawId() << std::endl;
0257       continue;
0258     }
0259     if (!it->second.isValid())
0260       continue;
0261     std::vector<TrajectoryMeasurement> mymeas =
0262         det.fastMeasurements(it->second, state, tracker_propagator, *estimator_);
0263     if (debug_)
0264       LogDebug("OutsideInMuonSeeder") << "Query on detector " << it->first->geographicalId().rawId() << " returned "
0265                                       << mymeas.size() << " measurements." << std::endl;
0266     for (std::vector<TrajectoryMeasurement>::const_iterator it2 = mymeas.begin(), ed2 = mymeas.end(); it2 != ed2;
0267          ++it2) {
0268       if (it2->recHit()->isValid())
0269         meas.push_back(*it2);
0270     }
0271   }
0272   int found = 0;
0273   std::sort(meas.begin(), meas.end(), TrajMeasLessEstim());
0274   for (std::vector<TrajectoryMeasurement>::const_iterator it2 = meas.begin(), ed2 = meas.end(); it2 != ed2; ++it2) {
0275     if (debug_) {
0276       LogDebug("OutsideInMuonSeeder") << "  inspecting Hit with chi2 = " << it2->estimate() << std::endl;
0277       LogDebug("OutsideInMuonSeeder") << "        track state     " << it2->forwardPredictedState().globalPosition()
0278                                       << std::endl;
0279       LogDebug("OutsideInMuonSeeder") << "        rechit position " << it2->recHit()->globalPosition() << std::endl;
0280     }
0281     TrajectoryStateOnSurface updated = updator_->update(it2->forwardPredictedState(), *it2->recHit());
0282     if (updated.isValid()) {
0283       if (debug_)
0284         LogDebug("OutsideInMuonSeeder") << "          --> updated state: x = " << updated.globalPosition()
0285                                         << ", p = " << updated.globalMomentum() << std::endl;
0286       edm::OwnVector<TrackingRecHit> seedHits;
0287       seedHits.push_back(*it2->recHit()->hit());
0288       PTrajectoryStateOnDet const &pstate =
0289           trajectoryStateTransform::persistentState(updated, it2->recHit()->geographicalId().rawId());
0290       TrajectorySeed seed(pstate, std::move(seedHits), oppositeToMomentum);
0291       out.push_back(seed);
0292       found++;
0293       if (found == hitsToTry_)
0294         break;
0295     }
0296   }
0297   return found;
0298 }
0299 
0300 void OutsideInMuonSeeder::doDebug(const reco::Track &tk) const {
0301   TrajectoryStateOnSurface tsos = trajectoryStateTransform::innerStateOnSurface(tk, *geometry_, magfield_);
0302   std::unique_ptr<Propagator> pmuon_cloned = SetPropagationDirection(*muonPropagator_, alongMomentum);
0303   for (unsigned int i = 0; i < tk.recHitsSize(); ++i) {
0304     const TrackingRecHit *hit = &*tk.recHit(i);
0305     const GeomDet *det = geometry_->idToDet(hit->geographicalId());
0306     if (det == nullptr)
0307       continue;
0308     if (i != 0)
0309       tsos = pmuon_cloned->propagate(tsos, det->surface());
0310     if (!tsos.isValid())
0311       continue;
0312     LogDebug("OutsideInMuonSeeder") << "  state " << i << " at x = " << tsos.globalPosition()
0313                                     << ", p = " << tsos.globalMomentum() << std::endl;
0314     if (hit->isValid()) {
0315       LogDebug("OutsideInMuonSeeder") << "         valid   rechit on detid " << hit->geographicalId().rawId()
0316                                       << std::endl;
0317     } else {
0318       LogDebug("OutsideInMuonSeeder") << "         invalid rechit on detid " << hit->geographicalId().rawId()
0319                                       << std::endl;
0320     }
0321   }
0322 }
0323 
0324 void OutsideInMuonSeeder::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0325   edm::ParameterSetDescription desc;
0326   desc.add<edm::InputTag>("src", edm::InputTag("muons"));
0327   desc.add<std::string>("cut", {});
0328   desc.add<int32_t>("layersToTry", 3);
0329   desc.add<int32_t>("hitsToTry", 3);
0330   desc.add<bool>("fromVertex", true);
0331   desc.add<double>("errorRescaleFactor", 2.0);
0332   desc.add<std::string>("trackerPropagator", {});
0333   desc.add<std::string>("muonPropagator", {});
0334   desc.add<edm::InputTag>("measurementTkEvent", edm::InputTag("MeasurementTrackerEvent"));
0335   desc.add<std::string>("hitCollector", {});
0336   desc.add<std::string>("updatorLabel", "KFUpdator");
0337   desc.add<double>("minEtaForTEC", 0.7);
0338   desc.add<double>("maxEtaForTOB", 1.8);
0339   desc.addUntracked<bool>("debug", false);
0340   descriptions.addWithDefaultLabel(desc);
0341 }
0342 
0343 #include "FWCore/Framework/interface/MakerMacros.h"
0344 DEFINE_FWK_MODULE(OutsideInMuonSeeder);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!