SpecialSeedGenerators/src/MuonReSeeder.cc

0001
0002 //
0003 //
0004
0005 /**
0006   \class    pat::MuonReSeeder MuonReSeeder.h "MuonAnalysis/MuonAssociators/interface/MuonReSeeder.h"
0007   \brief    Matcher of reconstructed objects to other reconstructed objects using the tracks inside them
0008
0009   \author   Giovanni Petrucciani
0010 */
0011
0012 #include "FWCore/Framework/interface/ConsumesCollector.h"
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017
0018 #include "DataFormats/Math/interface/deltaR.h"
0019
0020 #include "DataFormats/MuonReco/interface/Muon.h"
0021 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
0022 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0023 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0024 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0025 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0026 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0027 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0028 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0029
0030 class MuonReSeeder : public edm::stream::EDProducer<> {
0031 public:
0032   explicit MuonReSeeder(const edm::ParameterSet &iConfig);
0033   ~MuonReSeeder() override {}
0034
0035   void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override;
0036
0037 private:
0038   /// Labels for input collections
0039   edm::EDGetTokenT<edm::View<reco::Muon>> src_;
0040   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0041
0042   /// Muon selection
0043   StringCutObjectSelector<reco::Muon> selector_;
0044
0045   /// How many hits to keep from the muon trajectory
0046   int layersToKeep_;
0047
0048   /// Do inside-out
0049   bool insideOut_;
0050
0051   /// Dump deug information
0052   bool debug_;
0053
0054   /// Track Transformer
0055   TrackTransformer refitter_;
0056 };
0057
0058 MuonReSeeder::MuonReSeeder(const edm::ParameterSet &iConfig)
0059     : src_(consumes<edm::View<reco::Muon>>(iConfig.getParameter<edm::InputTag>("src"))),
0060       tTopoToken_(esConsumes()),
0061       selector_(iConfig.existsAs<std::string>("cut") ? iConfig.getParameter<std::string>("cut") : "", true),
0062       layersToKeep_(iConfig.getParameter<int32_t>("layersToKeep")),
0063       insideOut_(iConfig.getParameter<bool>("insideOut")),
0064       debug_(iConfig.getUntrackedParameter<bool>("debug", false)),
0065       refitter_(iConfig, consumesCollector()) {
0066   produces<std::vector<TrajectorySeed>>();
0067 }
0068
0069 void MuonReSeeder::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0070   using namespace edm;
0071   using namespace std;
0072
0073   refitter_.setServices(iSetup);
0074
0075   Handle<View<reco::Muon>> src;
0076   iEvent.getByToken(src_, src);
0077
0078   //Retrieve tracker topology from geometry
0079   const TrackerTopology &tTopo = iSetup.getData(tTopoToken_);
0080
0081   auto out = std::make_unique<std::vector<TrajectorySeed>>();
0082   unsigned int nsrc = src->size();
0083   out->reserve(nsrc);
0084
0085   for (View<reco::Muon>::const_iterator it = src->begin(), ed = src->end(); it != ed; ++it) {
0086     const reco::Muon &mu = *it;
0087     if (mu.track().isNull() || !selector_(mu))
0088       continue;
0089     std::vector<Trajectory> traj = refitter_.transform(*mu.track());
0090     if (traj.size() != 1)
0091       continue;
0092     edm::OwnVector<TrackingRecHit> seedHits;
0093     const std::vector<TrajectoryMeasurement> &tms = traj.front().measurements();
0094     TrajectoryStateOnSurface tsos;
0095     const TrackingRecHit *hit = nullptr;
0096     bool fromInside = (insideOut_ == (traj.front().direction() == alongMomentum));
0097     if (debug_) {
0098       std::cout << "Considering muon of pt " << mu.pt() << ", eta = " << mu.eta() << ", phi = " << mu.phi()
0099                 << std::endl;
0100       std::cout << "Trajectory is " << (traj.front().direction() == alongMomentum ? "along" : "opposite")
0101                 << " to momentum, so will start from " << (fromInside ? "inside" : "outside") << std::endl;
0102     }
0103     const TrajectoryMeasurement &tin = (fromInside ? tms.front() : tms.back());
0104     const TrajectoryMeasurement &tou = (fromInside ? tms.front() : tms.back());
0105     if (debug_) {
0106       std::cout << "IN state: subdetector   = " << tin.recHit()->geographicalId().subdetId() << std::endl;
0107       std::cout << "          global pos Rho  " << tin.updatedState().globalPosition().perp() << ", Z   "
0108                 << tin.updatedState().globalPosition().z() << std::endl;
0109       std::cout << "OU state: subdetector   = " << tou.recHit()->geographicalId().subdetId() << std::endl;
0110       std::cout << "          global pos Rho  " << tou.updatedState().globalPosition().perp() << ", Z   "
0111                 << tou.updatedState().globalPosition().z() << std::endl;
0112     }
0113     int lastSubdet = 0, lastLayer = -1;
0114     for (int i = (fromInside ? 0 : tms.size() - 1),
0115              end = (fromInside ? tms.size() : -1),
0116              step = (fromInside ? +1 : -1),
0117              taken = 0;
0118          (end - i) * step > 0;
0119          i += step) {
0120       const TrackingRecHit *lastHit = hit;
0121       hit = tms[i].recHit()->hit();
0122       if (debug_)
0123         std::cout << "  considering hit " << i << ": rechit on " << (hit ? hit->geographicalId().rawId() : -1)
0124                   << std::endl;
0125       if (!hit)
0126         continue;
0127       int subdet = hit->geographicalId().subdetId();
0128       int lay = tTopo.layer(hit->geographicalId());
0129       if (subdet != lastSubdet || lay != lastLayer) {
0130         // I'm on a new layer
0131         if (lastHit != nullptr && taken == layersToKeep_) {
0132           // I've had enough layers, I can stop here
0133           hit = lastHit;
0134           break;
0135         }
0136         lastSubdet = subdet;
0137         lastLayer = lay;
0138         taken++;
0139       }
0140       seedHits.push_back(*hit);
0141       tsos = tms[i].updatedState().isValid()
0142                  ? tms[i].updatedState()
0143                  : (abs(i - end) < abs(i) ? tms[i].forwardPredictedState() : tms[i].backwardPredictedState());
0144       if (debug_) {
0145         std::cout << "     hit  : subdetector   = " << tms[i].recHit()->geographicalId().subdetId() << std::endl;
0146         if (hit->isValid()) {
0147           std::cout << "            global pos Rho  " << tms[i].recHit()->globalPosition().perp() << ", Z   "
0148                     << tms[i].recHit()->globalPosition().z() << std::endl;
0149         } else {
0150           std::cout << "            invalid tracking rec hit, so no global position" << std::endl;
0151         }
0152         std::cout << "     state: global pos Rho  " << tsos.globalPosition().perp() << ", Z   "
0153                   << tsos.globalPosition().z() << std::endl;
0154       }
0155     }
0156     if (!tsos.isValid())
0157       continue;
0158     PTrajectoryStateOnDet const &PTraj = trajectoryStateTransform::persistentState(tsos, hit->geographicalId().rawId());
0159     TrajectorySeed seed(PTraj, std::move(seedHits), insideOut_ ? alongMomentum : oppositeToMomentum);
0160     out->push_back(seed);
0161   }
0162
0163   iEvent.put(std::move(out));
0164 }
0165
0166 #include "FWCore/Framework/interface/MakerMacros.h"
0167 DEFINE_FWK_MODULE(MuonReSeeder);