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::beginRun(edm::Run const& run, edm::EventSetup const& es) {
0033   //region builder
0034 }
0035
0036 void FastTSGFromL2Muon::produce(edm::Event& ev, const edm::EventSetup& es) {
0037   // Initialize the output product
0038   std::unique_ptr<L3MuonTrajectorySeedCollection> result(new L3MuonTrajectorySeedCollection());
0039
0040   // Region builder
0041   theRegionBuilder->setEvent(ev, es);
0042
0043   // Retrieve the Monte Carlo truth (SimTracks)
0044   const edm::Handle<edm::SimTrackContainer>& theSimTracks = ev.getHandle(simTrackToken_);
0045
0046   // Retrieve L2 muon collection
0047   const edm::Handle<reco::TrackCollection>& l2muonH = ev.getHandle(l2TrackToken_);
0048
0049   // Retrieve Seed collection
0050   unsigned seedCollections = theSeedCollectionLabels.size();
0051   std::vector<edm::Handle<edm::View<TrajectorySeed> > > theSeeds;
0052   theSeeds.resize(seedCollections);
0053   unsigned seed_size = 0;
0054   for (unsigned iseed = 0; iseed < seedCollections; ++iseed) {
0055     ev.getByToken(seedToken_[iseed], theSeeds[iseed]);
0056     seed_size += theSeeds[iseed]->size();
0057   }
0058
0059   // Loop on L2 muons
0060   unsigned int imu = 0;
0061   unsigned int imuMax = l2muonH->size();
0062   edm::LogVerbatim("FastTSGFromL2Muon") << "Found " << imuMax << " L2 muons";
0063   for (; imu != imuMax; ++imu) {
0064     // Make a ref to l2 muon
0065     reco::TrackRef muRef(l2muonH, imu);
0066
0067     // Cut on muons with low momenta
0068     if (muRef->pt() < thePtCut || muRef->innerMomentum().Rho() < thePtCut || muRef->innerMomentum().R() < 2.5)
0069       continue;
0070
0071     // Define the region of interest
0072     std::unique_ptr<RectangularEtaPhiTrackingRegion> region = theRegionBuilder->region(muRef);
0073
0074     // Copy the collection of seeds (ahem, this is time consuming!)
0075     std::vector<TrajectorySeed> tkSeeds;
0076     std::set<unsigned> tkIds;
0077     tkSeeds.reserve(seed_size);
0078     for (unsigned iseed = 0; iseed < seedCollections; ++iseed) {
0079       edm::Handle<edm::View<TrajectorySeed> > aSeedCollection = theSeeds[iseed];
0080       unsigned nSeeds = aSeedCollection->size();
0081       for (unsigned seednr = 0; seednr < nSeeds; ++seednr) {
0082         // The seed
0083         const BasicTrajectorySeed* aSeed = &((*aSeedCollection)[seednr]);
0084
0085         // The SimTrack id associated to the first hit of the Seed
0086         int simTrackId = static_cast<FastTrackerRecHit const&>(*aSeed->recHits().begin()).simTrackId(0);
0087
0088         // Track already associated to a seed
0089         std::set<unsigned>::iterator tkId = tkIds.find(simTrackId);
0090         if (tkId != tkIds.end())
0091           continue;
0092
0093         const SimTrack& theSimTrack = (*theSimTracks)[simTrackId];
0094
0095         if (clean(muRef, region.get(), aSeed, theSimTrack))
0096           tkSeeds.push_back(*aSeed);
0097         tkIds.insert(simTrackId);
0098
0099       }  // End loop on seeds
0100
0101     }  // End loop on seed collections
0102
0103     // Now create the Muon Trajectory Seed
0104     unsigned int is = 0;
0105     unsigned int isMax = tkSeeds.size();
0106     for (; is != isMax; ++is) {
0107       result->push_back(L3MuonTrajectorySeed(tkSeeds[is], muRef));
0108     }  // End of tk seed loop
0109
0110   }  // End of l2 muon loop
0111
0112   edm::LogVerbatim("FastTSGFromL2Muon") << "Found " << result->size() << " seeds for muons";
0113
0114   //put in the event
0115   ev.put(std::move(result));
0116 }
0117
0118 bool FastTSGFromL2Muon::clean(reco::TrackRef muRef,
0119                               RectangularEtaPhiTrackingRegion* region,
0120                               const BasicTrajectorySeed* aSeed,
0121                               const SimTrack& theSimTrack) {
0122   // Eta cleaner
0123   const PixelRecoRange<float>& etaRange = region->etaRange();
0124   double etaSeed = theSimTrack.momentum().Eta();
0125   double etaLimit = (fabs(fabs(etaRange.max()) - fabs(etaRange.mean())) < 0.05)
0126                         ? 0.05
0127                         : fabs(fabs(etaRange.max()) - fabs(etaRange.mean()));
0128   bool inEtaRange = etaSeed >= (etaRange.mean() - etaLimit) && etaSeed <= (etaRange.mean() + etaLimit);
0129   if (!inEtaRange)
0130     return false;
0131
0132   // Phi cleaner
0133   const TkTrackingRegionsMargin<float>& phiMargin = region->phiMargin();
0134   double phiSeed = theSimTrack.momentum().Phi();
0135   double phiLimit = (phiMargin.right() < 0.05) ? 0.05 : phiMargin.right();
0136   bool inPhiRange = (fabs(deltaPhi(phiSeed, double(region->direction().phi()))) < phiLimit);
0137   if (!inPhiRange)
0138     return false;
0139
0140   // pt cleaner
0141   double ptSeed = std::sqrt(theSimTrack.momentum().Perp2());
0142   double ptMin = (region->ptMin() > 3.5) ? 3.5 : region->ptMin();
0143   bool inPtRange = ptSeed >= ptMin && ptSeed <= 2 * (muRef->pt());
0144   return inPtRange;
0145 }