Muons/plugins/FastTSGFromL2Muon.cc

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