Muons/plugins/FastTSGFromL2Muon.cc

0001 #include "FastSimulation/Muons/plugins/FastTSGFromL2Muon.h"
0002
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeedCollection.h"
0007 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
0008 #include "DataFormats/TrackReco/interface/Track.h"
0009 #include "DataFormats/Math/interface/deltaPhi.h"
0010
0011 #include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
0012 #include "RecoMuon/GlobalTrackingTools/interface/MuonTrackingRegionBuilder.h"
0013
0014 #include <set>
0015
0016 FastTSGFromL2Muon::FastTSGFromL2Muon(const edm::ParameterSet& cfg)
0017     : thePtCut(cfg.getParameter<double>("PtCut")),
0018       theL2CollectionLabel(cfg.getParameter<edm::InputTag>("MuonCollectionLabel")),
0019       theSeedCollectionLabels(cfg.getParameter<std::vector<edm::InputTag> >("SeedCollectionLabels")),
0020       theSimTrackCollectionLabel(cfg.getParameter<edm::InputTag>("SimTrackCollectionLabel")),
0021       simTrackToken_(consumes<edm::SimTrackContainer>(theSimTrackCollectionLabel)),
0022       l2TrackToken_(consumes<reco::TrackCollection>(theL2CollectionLabel)) {
0023   produces<L3MuonTrajectorySeedCollection>();
0024
0025   for (auto& seed : theSeedCollectionLabels)
0026     seedToken_.emplace_back(consumes<edm::View<TrajectorySeed> >(seed));
0027
0028   edm::ParameterSet regionBuilderPSet = cfg.getParameter<edm::ParameterSet>("MuonTrackingRegionBuilder");
0029   theRegionBuilder = std::make_unique<MuonTrackingRegionBuilder>(regionBuilderPSet, consumesCollector());
0030 }
0031
0032 void FastTSGFromL2Muon::produce(edm::Event& ev, const edm::EventSetup& es) {
0033   // Initialize the output product
0034   std::unique_ptr<L3MuonTrajectorySeedCollection> result(new L3MuonTrajectorySeedCollection());
0035
0036   // Region builder
0037   theRegionBuilder->setEvent(ev, es);
0038
0039   // Retrieve the Monte Carlo truth (SimTracks)
0040   const edm::Handle<edm::SimTrackContainer>& theSimTracks = ev.getHandle(simTrackToken_);
0041
0042   // Retrieve L2 muon collection
0043   const edm::Handle<reco::TrackCollection>& l2muonH = ev.getHandle(l2TrackToken_);
0044
0045   // Retrieve Seed collection
0046   unsigned seedCollections = theSeedCollectionLabels.size();
0047   std::vector<edm::Handle<edm::View<TrajectorySeed> > > theSeeds;
0048   theSeeds.resize(seedCollections);
0049   unsigned seed_size = 0;
0050   for (unsigned iseed = 0; iseed < seedCollections; ++iseed) {
0051     ev.getByToken(seedToken_[iseed], theSeeds[iseed]);
0052     seed_size += theSeeds[iseed]->size();
0053   }
0054
0055   // Loop on L2 muons
0056   unsigned int imu = 0;
0057   unsigned int imuMax = l2muonH->size();
0058   edm::LogVerbatim("FastTSGFromL2Muon") << "Found " << imuMax << " L2 muons";
0059   for (; imu != imuMax; ++imu) {
0060     // Make a ref to l2 muon
0061     reco::TrackRef muRef(l2muonH, imu);
0062
0063     // Cut on muons with low momenta
0064     if (muRef->pt() < thePtCut || muRef->innerMomentum().Rho() < thePtCut || muRef->innerMomentum().R() < 2.5)
0065       continue;
0066
0067     // Define the region of interest
0068     std::unique_ptr<RectangularEtaPhiTrackingRegion> region = theRegionBuilder->region(muRef);
0069
0070     // Copy the collection of seeds (ahem, this is time consuming!)
0071     std::vector<TrajectorySeed> tkSeeds;
0072     std::set<unsigned> tkIds;
0073     tkSeeds.reserve(seed_size);
0074     for (unsigned iseed = 0; iseed < seedCollections; ++iseed) {
0075       edm::Handle<edm::View<TrajectorySeed> > aSeedCollection = theSeeds[iseed];
0076       unsigned nSeeds = aSeedCollection->size();
0077       for (unsigned seednr = 0; seednr < nSeeds; ++seednr) {
0078         // The seed
0079         const BasicTrajectorySeed* aSeed = &((*aSeedCollection)[seednr]);
0080
0081         // The SimTrack id associated to the first hit of the Seed
0082         int simTrackId = static_cast<FastTrackerRecHit const&>(*aSeed->recHits().begin()).simTrackId(0);
0083
0084         // Track already associated to a seed
0085         std::set<unsigned>::iterator tkId = tkIds.find(simTrackId);
0086         if (tkId != tkIds.end())
0087           continue;
0088
0089         const SimTrack& theSimTrack = (*theSimTracks)[simTrackId];
0090
0091         if (clean(muRef, region.get(), aSeed, theSimTrack))
0092           tkSeeds.push_back(*aSeed);
0093         tkIds.insert(simTrackId);
0094
0095       }  // End loop on seeds
0096
0097     }  // End loop on seed collections
0098
0099     // Now create the Muon Trajectory Seed
0100     unsigned int is = 0;
0101     unsigned int isMax = tkSeeds.size();
0102     for (; is != isMax; ++is) {
0103       result->push_back(L3MuonTrajectorySeed(tkSeeds[is], muRef));
0104     }  // End of tk seed loop
0105
0106   }  // End of l2 muon loop
0107
0108   edm::LogVerbatim("FastTSGFromL2Muon") << "Found " << result->size() << " seeds for muons";
0109
0110   //put in the event
0111   ev.put(std::move(result));
0112 }
0113
0114 bool FastTSGFromL2Muon::clean(reco::TrackRef muRef,
0115                               RectangularEtaPhiTrackingRegion* region,
0116                               const BasicTrajectorySeed* aSeed,
0117                               const SimTrack& theSimTrack) {
0118   // Eta cleaner
0119   const PixelRecoRange<float>& etaRange = region->etaRange();
0120   double etaSeed = theSimTrack.momentum().Eta();
0121   double etaLimit = (fabs(fabs(etaRange.max()) - fabs(etaRange.mean())) < 0.05)
0122                         ? 0.05
0123                         : fabs(fabs(etaRange.max()) - fabs(etaRange.mean()));
0124   bool inEtaRange = etaSeed >= (etaRange.mean() - etaLimit) && etaSeed <= (etaRange.mean() + etaLimit);
0125   if (!inEtaRange)
0126     return false;
0127
0128   // Phi cleaner
0129   const TkTrackingRegionsMargin<float>& phiMargin = region->phiMargin();
0130   double phiSeed = theSimTrack.momentum().Phi();
0131   double phiLimit = (phiMargin.right() < 0.05) ? 0.05 : phiMargin.right();
0132   bool inPhiRange = (fabs(deltaPhi(phiSeed, double(region->direction().phi()))) < phiLimit);
0133   if (!inPhiRange)
0134     return false;
0135
0136   // pt cleaner
0137   double ptSeed = std::sqrt(theSimTrack.momentum().Perp2());
0138   double ptMin = (region->ptMin() > 3.5) ? 3.5 : region->ptMin();
0139   bool inPtRange = ptSeed >= ptMin && ptSeed <= 2 * (muRef->pt());
0140   return inPtRange;
0141 }