Muons/plugins/FastTSGFromIOHit.cc

0001 #include "FastSimulation/Muons/plugins/FastTSGFromIOHit.h"
0002
0003 /** \class FastTSGFromIOHit
0004  *
0005  *  Emulate TSGFromIOHit in RecoMuon
0006  *
0007  *  \author Adam Everett - Purdue University
0008  */
0009
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0013 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
0014 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeed.h"
0015 #include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
0016 #include "DataFormats/Math/interface/deltaPhi.h"
0017
0018 FastTSGFromIOHit::FastTSGFromIOHit(const edm::ParameterSet& iConfig, edm::ConsumesCollector& iC) {
0019   theCategory = "FastSimulation|Muons||FastTSGFromIOHit";
0020
0021   thePtCut = iConfig.getParameter<double>("PtCut");
0022
0023   simTracksTk = iC.consumes<edm::SimTrackContainer>(iConfig.getParameter<edm::InputTag>("SimTrackCollectionLabel"));
0024   const auto& seedLabels = iConfig.getParameter<std::vector<edm::InputTag> >("SeedCollectionLabels");
0025   for (const auto& seedLabel : seedLabels) {
0026     seedsTks.push_back(iC.consumes<TrajectorySeedCollection>(seedLabel));
0027   }
0028 }
0029
0030 FastTSGFromIOHit::~FastTSGFromIOHit() { LogTrace(theCategory) << " FastTSGFromIOHit dtor called "; }
0031
0032 void FastTSGFromIOHit::trackerSeeds(const TrackCand& staMuon,
0033                                     const TrackingRegion& region,
0034                                     const TrackerTopology* tTopo,
0035                                     std::vector<TrajectorySeed>& result) {
0036   // Make a ref to l2 muon
0037   reco::TrackRef muRef(staMuon.second);
0038
0039   // Cut on muons with low momenta
0040   if (muRef->pt() < thePtCut || muRef->innerMomentum().Rho() < thePtCut || muRef->innerMomentum().R() < 2.5) {
0041     return;
0042   }
0043
0044   // Retrieve the Monte Carlo truth (SimTracks)
0045   edm::Handle<edm::SimTrackContainer> simTracks;
0046   getEvent()->getByToken(simTracksTk, simTracks);
0047
0048   // Retrieve Seed collection
0049   std::vector<edm::Handle<TrajectorySeedCollection> > seedCollections;
0050   seedCollections.resize(seedsTks.size());
0051   for (unsigned iSeed = 0; iSeed < seedsTks.size(); iSeed++) {
0052     getEvent()->getByToken(seedsTks[iSeed], seedCollections[iSeed]);
0053   }
0054
0055   // cast the tracking region
0056   const RectangularEtaPhiTrackingRegion& regionRef = dynamic_cast<const RectangularEtaPhiTrackingRegion&>(region);
0057
0058   // select and store seeds
0059   std::set<unsigned> simTrackIds;
0060   for (const auto& seeds : seedCollections) {
0061     for (const auto& seed : *seeds) {
0062       // Find the simtrack corresponding to the seed
0063       int simTrackId = static_cast<FastTrackerRecHit const&>(*seed.recHits().begin()).simTrackId(0);
0064       const SimTrack& simTrack = (*simTracks)[simTrackId];
0065
0066       // skip if simTrack already associated to a seed
0067       if (simTrackIds.find(simTrackId) != simTrackIds.end()) {
0068         continue;
0069       }
0070       simTrackIds.insert(simTrackId);
0071
0072       // filter seed
0073       if (!clean(muRef, regionRef, &seed, simTrack)) {
0074         continue;
0075       }
0076
0077       // store results
0078       result.push_back(L3MuonTrajectorySeed(seed, muRef));
0079     }
0080   }
0081 }
0082
0083 bool FastTSGFromIOHit::clean(reco::TrackRef muRef,
0084                              const RectangularEtaPhiTrackingRegion& region,
0085                              const TrajectorySeed* aSeed,
0086                              const SimTrack& theSimTrack) {
0087   // Eta cleaner
0088   const PixelRecoRange<float>& etaRange = region.etaRange();
0089
0090   double etaSeed = theSimTrack.momentum().Eta();
0091   double etaLimit = (fabs(fabs(etaRange.max()) - fabs(etaRange.mean())) < 0.05)
0092                         ? 0.05
0093                         : fabs(fabs(etaRange.max()) - fabs(etaRange.mean()));
0094   bool inEtaRange = etaSeed >= (etaRange.mean() - etaLimit) && etaSeed <= (etaRange.mean() + etaLimit);
0095   if (!inEtaRange) {
0096     return false;
0097   }
0098
0099   // Phi cleaner
0100   const TkTrackingRegionsMargin<float>& phiMargin = region.phiMargin();
0101   double phiSeed = theSimTrack.momentum().Phi();
0102   double phiLimit = (phiMargin.right() < 0.05) ? 0.05 : phiMargin.right();
0103   bool inPhiRange = (fabs(deltaPhi(phiSeed, double(region.direction().phi()))) < phiLimit);
0104   if (!inPhiRange) {
0105     return false;
0106   }
0107
0108   // pt cleaner
0109   double ptSeed = std::sqrt(theSimTrack.momentum().Perp2());
0110   double ptMin = (region.ptMin() > 3.5) ? 3.5 : region.ptMin();
0111   bool inPtRange = ptSeed >= ptMin && ptSeed <= 2 * (muRef->pt());
0112   if (!inPtRange) {
0113     return false;
0114   }
0115   return true;
0116 }