Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:22:47 

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 "FWCore/Framework/interface/stream/EDProducer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013 
0014 #include "DataFormats/Math/interface/deltaR.h"
0015 
0016 #include "DataFormats/MuonReco/interface/Muon.h"
0017 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0018 #include "MagneticField/Engine/interface/MagneticField.h"
0019 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0020 
0021 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
0022 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0023 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0024 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0025 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0026 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0027 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0028 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0029 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
0030 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0031 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0032 #include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
0033 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
0034 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0035 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0036 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0037 
0038 class OutsideInMuonSeeder final : public edm::stream::EDProducer<> {
0039 public:
0040   explicit OutsideInMuonSeeder(const edm::ParameterSet &iConfig);
0041   ~OutsideInMuonSeeder() override {}
0042 
0043   void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override;
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_(consumes<MeasurementTrackerEvent>(edm::InputTag("MeasurementTrackerEvent"))),
0104       estimatorToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("hitCollector")))),
0105       updatorToken_(esConsumes(edm::ESInputTag("", "KFUpdator"))),
0106       magfieldToken_(esConsumes()),
0107       geometryToken_(esConsumes()),
0108       tkGeometryToken_(esConsumes()),
0109       minEtaForTEC_(iConfig.getParameter<double>("minEtaForTEC")),
0110       maxEtaForTOB_(iConfig.getParameter<double>("maxEtaForTOB")),
0111       debug_(iConfig.getUntrackedParameter<bool>("debug", false)) {
0112   produces<std::vector<TrajectorySeed>>();
0113 }
0114 
0115 void OutsideInMuonSeeder::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0116   using namespace edm;
0117   using namespace std;
0118 
0119   magfield_ = &iSetup.getData(magfieldToken_);
0120   auto const &trackerPropagator = iSetup.getData(trackerPropagatorToken_);
0121   muonPropagator_ = &iSetup.getData(muonPropagatorToken_);
0122   geometry_ = &iSetup.getData(geometryToken_);
0123   estimator_ = &iSetup.getData(estimatorToken_);
0124   updator_ = &iSetup.getData(updatorToken_);
0125 
0126   Handle<MeasurementTrackerEvent> measurementTracker;
0127   iEvent.getByToken(measurementTrackerTag_, measurementTracker);
0128 
0129   const auto &tmpTkGeometry = iSetup.getData(tkGeometryToken_);
0130 
0131   Handle<View<reco::Muon>> src;
0132   iEvent.getByToken(src_, src);
0133 
0134   auto out = std::make_unique<std::vector<TrajectorySeed>>();
0135 
0136   for (auto const &mu : *src) {
0137     if (mu.outerTrack().isNull() || !selector_(mu))
0138       continue;
0139     if (debug_ && mu.innerTrack().isNonnull())
0140       doDebug(*mu.innerTrack());
0141 
0142     // Better clone here and not directly into doLayer to avoid
0143     // useless clone/destroy operations to set, in the end, the
0144     // very same direction every single time.
0145     std::unique_ptr<Propagator> pmuon_cloned =
0146         SetPropagationDirection(*muonPropagator_, fromVertex_ ? alongMomentum : oppositeToMomentum);
0147     std::unique_ptr<Propagator> ptracker_cloned = SetPropagationDirection(trackerPropagator, alongMomentum);
0148 
0149     int sizeBefore = out->size();
0150     if (debug_)
0151       LogDebug("OutsideInMuonSeeder") << "\n\n\nSeeding for muon of pt " << mu.pt() << ", eta " << mu.eta() << ", phi "
0152                                       << mu.phi() << std::endl;
0153     const reco::Track &tk = *mu.outerTrack();
0154 
0155     TrajectoryStateOnSurface state =
0156         fromVertex_ ? TrajectoryStateOnSurface(trajectoryStateTransform::initialFreeState(tk, magfield_))
0157                     : trajectoryStateTransform::innerStateOnSurface(tk, *geometry_, magfield_);
0158 
0159     if (std::abs(tk.eta()) < maxEtaForTOB_) {
0160       std::vector<BarrelDetLayer const *> const &tob = measurementTracker->geometricSearchTracker()->tobLayers();
0161       int found = 0;
0162       int iLayer = tob.size();
0163       if (iLayer == 0)
0164         LogError("OutsideInMuonSeeder") << "TOB has no layers.";
0165 
0166       for (auto it = tob.rbegin(), ed = tob.rend(); it != ed; ++it, --iLayer) {
0167         if (debug_)
0168           LogDebug("OutsideInMuonSeeder") << "\n ==== Trying TOB " << iLayer << " ====" << std::endl;
0169         if (doLayer(**it, state, *out, *(pmuon_cloned.get()), *(ptracker_cloned.get()), *measurementTracker)) {
0170           if (++found == layersToTry_)
0171             break;
0172         }
0173       }
0174     }
0175     if (tk.eta() > minEtaForTEC_) {
0176       const auto &forwLayers = tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC)
0177                                    ? measurementTracker->geometricSearchTracker()->posTidLayers()
0178                                    : measurementTracker->geometricSearchTracker()->posTecLayers();
0179       if (tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC)) {
0180         LogDebug("OutsideInMuonSeeder") << "\n We are using the Phase2 Outer Tracker (defined as a TID+). ";
0181       }
0182       LogTrace("OutsideInMuonSeeder") << "\n ==== TEC+ tot layers " << forwLayers.size() << " ====" << std::endl;
0183       int found = 0;
0184       int iLayer = forwLayers.size();
0185       if (iLayer == 0)
0186         LogError("OutsideInMuonSeeder") << "TEC+ has no layers.";
0187 
0188       if (debug_)
0189         LogDebug("OutsideInMuonSeeder") << "\n ==== Tot layers " << forwLayers.size() << " ====" << std::endl;
0190       for (auto it = forwLayers.rbegin(), ed = forwLayers.rend(); it != ed; ++it, --iLayer) {
0191         if (debug_)
0192           LogDebug("OutsideInMuonSeeder") << "\n ==== Trying Forward Layer +" << +iLayer << " ====" << std::endl;
0193         if (doLayer(**it, state, *out, *(pmuon_cloned.get()), *(ptracker_cloned.get()), *measurementTracker)) {
0194           if (++found == layersToTry_)
0195             break;
0196         }
0197       }
0198     }
0199     if (tk.eta() < -minEtaForTEC_) {
0200       const auto &forwLayers = tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC)
0201                                    ? measurementTracker->geometricSearchTracker()->negTidLayers()
0202                                    : measurementTracker->geometricSearchTracker()->negTecLayers();
0203       if (tmpTkGeometry.isThere(GeomDetEnumerators::P2OTEC)) {
0204         LogDebug("OutsideInMuonSeeder") << "\n We are using the Phase2 Outer Tracker (defined as a TID-). ";
0205       }
0206       LogTrace("OutsideInMuonSeeder") << "\n ==== TEC- tot layers " << forwLayers.size() << " ====" << std::endl;
0207       int found = 0;
0208       int iLayer = forwLayers.size();
0209       if (iLayer == 0)
0210         LogError("OutsideInMuonSeeder") << "TEC- has no layers.";
0211 
0212       if (debug_)
0213         LogDebug("OutsideInMuonSeeder") << "\n ==== Tot layers " << forwLayers.size() << " ====" << std::endl;
0214       for (auto it = forwLayers.rbegin(), ed = forwLayers.rend(); it != ed; ++it, --iLayer) {
0215         if (debug_)
0216           LogDebug("OutsideInMuonSeeder") << "\n ==== Trying Forward Layer -" << -iLayer << " ====" << std::endl;
0217         if (doLayer(**it, state, *out, *(pmuon_cloned.get()), *(ptracker_cloned.get()), *measurementTracker)) {
0218           if (++found == layersToTry_)
0219             break;
0220         }
0221       }
0222     }
0223     if (debug_)
0224       LogDebug("OutsideInMuonSeeder") << "Outcome of seeding for muon of pt " << mu.pt() << ", eta " << mu.eta()
0225                                       << ", phi " << mu.phi() << ": found " << (out->size() - sizeBefore) << " seeds."
0226                                       << std::endl;
0227   }
0228 
0229   iEvent.put(std::move(out));
0230 }
0231 
0232 int OutsideInMuonSeeder::doLayer(const GeometricSearchDet &layer,
0233                                  const TrajectoryStateOnSurface &state,
0234                                  std::vector<TrajectorySeed> &out,
0235                                  const Propagator &muon_propagator,
0236                                  const Propagator &tracker_propagator,
0237                                  const MeasurementTrackerEvent &measurementTracker) const {
0238   TrajectoryStateOnSurface onLayer(state);
0239   onLayer.rescaleError(errorRescaling_);
0240   std::vector<GeometricSearchDet::DetWithState> dets;
0241   layer.compatibleDetsV(onLayer, muon_propagator, *estimator_, dets);
0242 
0243   if (debug_) {
0244     LogDebug("OutsideInMuonSeeder") << "Query on layer around x = " << onLayer.globalPosition()
0245                                     << " with local pos error " << sqrt(onLayer.localError().positionError().xx())
0246                                     << " ,  " << sqrt(onLayer.localError().positionError().yy()) << " ,  "
0247                                     << " returned " << dets.size() << " compatible detectors" << std::endl;
0248   }
0249 
0250   std::vector<TrajectoryMeasurement> meas;
0251   for (std::vector<GeometricSearchDet::DetWithState>::const_iterator it = dets.begin(), ed = dets.end(); it != ed;
0252        ++it) {
0253     MeasurementDetWithData det = measurementTracker.idToDet(it->first->geographicalId());
0254     if (det.isNull()) {
0255       std::cerr << "BOGUS detid " << it->first->geographicalId().rawId() << std::endl;
0256       continue;
0257     }
0258     if (!it->second.isValid())
0259       continue;
0260     std::vector<TrajectoryMeasurement> mymeas =
0261         det.fastMeasurements(it->second, state, tracker_propagator, *estimator_);
0262     if (debug_)
0263       LogDebug("OutsideInMuonSeeder") << "Query on detector " << it->first->geographicalId().rawId() << " returned "
0264                                       << mymeas.size() << " measurements." << std::endl;
0265     for (std::vector<TrajectoryMeasurement>::const_iterator it2 = mymeas.begin(), ed2 = mymeas.end(); it2 != ed2;
0266          ++it2) {
0267       if (it2->recHit()->isValid())
0268         meas.push_back(*it2);
0269     }
0270   }
0271   int found = 0;
0272   std::sort(meas.begin(), meas.end(), TrajMeasLessEstim());
0273   for (std::vector<TrajectoryMeasurement>::const_iterator it2 = meas.begin(), ed2 = meas.end(); it2 != ed2; ++it2) {
0274     if (debug_) {
0275       LogDebug("OutsideInMuonSeeder") << "  inspecting Hit with chi2 = " << it2->estimate() << std::endl;
0276       LogDebug("OutsideInMuonSeeder") << "        track state     " << it2->forwardPredictedState().globalPosition()
0277                                       << std::endl;
0278       LogDebug("OutsideInMuonSeeder") << "        rechit position " << it2->recHit()->globalPosition() << std::endl;
0279     }
0280     TrajectoryStateOnSurface updated = updator_->update(it2->forwardPredictedState(), *it2->recHit());
0281     if (updated.isValid()) {
0282       if (debug_)
0283         LogDebug("OutsideInMuonSeeder") << "          --> updated state: x = " << updated.globalPosition()
0284                                         << ", p = " << updated.globalMomentum() << std::endl;
0285       edm::OwnVector<TrackingRecHit> seedHits;
0286       seedHits.push_back(*it2->recHit()->hit());
0287       PTrajectoryStateOnDet const &pstate =
0288           trajectoryStateTransform::persistentState(updated, it2->recHit()->geographicalId().rawId());
0289       TrajectorySeed seed(pstate, std::move(seedHits), oppositeToMomentum);
0290       out.push_back(seed);
0291       found++;
0292       if (found == hitsToTry_)
0293         break;
0294     }
0295   }
0296   return found;
0297 }
0298 
0299 void OutsideInMuonSeeder::doDebug(const reco::Track &tk) const {
0300   TrajectoryStateOnSurface tsos = trajectoryStateTransform::innerStateOnSurface(tk, *geometry_, magfield_);
0301   std::unique_ptr<Propagator> pmuon_cloned = SetPropagationDirection(*muonPropagator_, alongMomentum);
0302   for (unsigned int i = 0; i < tk.recHitsSize(); ++i) {
0303     const TrackingRecHit *hit = &*tk.recHit(i);
0304     const GeomDet *det = geometry_->idToDet(hit->geographicalId());
0305     if (det == nullptr)
0306       continue;
0307     if (i != 0)
0308       tsos = pmuon_cloned->propagate(tsos, det->surface());
0309     if (!tsos.isValid())
0310       continue;
0311     LogDebug("OutsideInMuonSeeder") << "  state " << i << " at x = " << tsos.globalPosition()
0312                                     << ", p = " << tsos.globalMomentum() << std::endl;
0313     if (hit->isValid()) {
0314       LogDebug("OutsideInMuonSeeder") << "         valid   rechit on detid " << hit->geographicalId().rawId()
0315                                       << std::endl;
0316     } else {
0317       LogDebug("OutsideInMuonSeeder") << "         invalid rechit on detid " << hit->geographicalId().rawId()
0318                                       << std::endl;
0319     }
0320   }
0321 }
0322 
0323 #include "FWCore/Framework/interface/MakerMacros.h"
0324 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!