L3MuonIsolationProducer/src/L3MuonIsolationProducer.cc

0001 #include "L3MuonIsolationProducer.h"
0002
0003 // Framework
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/ConsumesCollector.h"
0007 #include "DataFormats/Common/interface/Handle.h"
0008
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 #include "DataFormats/Common/interface/AssociationMap.h"
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0016
0017 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0018 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
0019
0020 #include "RecoMuon/MuonIsolation/interface/Range.h"
0021 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
0022
0023 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractor.h"
0024 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractorFactory.h"
0025
0026 #include "L3NominalEfficiencyConfigurator.h"
0027
0028 #include <string>
0029
0030 using namespace edm;
0031 using namespace std;
0032 using namespace reco;
0033 using namespace muonisolation;
0034
0035 /// constructor with config
0036 L3MuonIsolationProducer::L3MuonIsolationProducer(const ParameterSet& par)
0037     : theConfig(par),
0038       theMuonCollectionLabel(par.getParameter<InputTag>("inputMuonCollection")),
0039       optOutputIsoDeposits(par.getParameter<bool>("OutputMuIsoDeposits")),
0040       theTrackPt_Min(-1) {
0041   LogDebug("RecoMuon|L3MuonIsolationProducer") << " L3MuonIsolationProducer CTOR";
0042
0043   theMuonCollectionToken = consumes<RecoChargedCandidateCollection>(theMuonCollectionLabel);
0044
0045   if (optOutputIsoDeposits)
0046     produces<reco::IsoDepositMap>();
0047   produces<edm::ValueMap<bool>>();
0048
0049   //
0050   // Extractor
0051   //
0052   edm::ParameterSet extractorPSet = theConfig.getParameter<edm::ParameterSet>("ExtractorPSet");
0053   //! get min pt for the track to go into sumPt
0054   theTrackPt_Min = theConfig.getParameter<double>("TrackPt_Min");
0055   std::string extractorName = extractorPSet.getParameter<std::string>("ComponentName");
0056   theExtractor = std::unique_ptr<reco::isodeposit::IsoDepositExtractor>{
0057       IsoDepositExtractorFactory::get()->create(extractorName, extractorPSet, consumesCollector())};
0058   std::string depositType = extractorPSet.getUntrackedParameter<std::string>("DepositLabel");
0059
0060   //
0061   // Cuts
0062   //
0063   edm::ParameterSet cutsPSet = theConfig.getParameter<edm::ParameterSet>("CutsPSet");
0064   std::string cutsName = cutsPSet.getParameter<std::string>("ComponentName");
0065   if (cutsName == "SimpleCuts") {
0066     theCuts = Cuts(cutsPSet);
0067   } else if (
0068       //        (cutsName== "L3NominalEfficiencyCuts_PXLS" && depositType=="PXLS")
0069       //     || (cutsName== "L3NominalEfficiencyCuts_TRKS" && depositType=="TRKS")
0070       //! test cutsName only. The depositType is informational only (has not been used so far) [VK]
0071       (cutsName == "L3NominalEfficiencyCuts_PXLS") || (cutsName == "L3NominalEfficiencyCuts_TRKS")) {
0072     theCuts = L3NominalEfficiencyConfigurator(cutsPSet).cuts();
0073   } else {
0074     LogError("L3MuonIsolationProducer::beginJob") << "cutsName: " << cutsPSet << " is not recognized:"
0075                                                   << " theCuts not set!";
0076   }
0077   LogTrace("") << theCuts.print();
0078
0079   // (kludge) additional cut on the number of tracks
0080   theMaxNTracks = cutsPSet.getParameter<int>("maxNTracks");
0081   theApplyCutsORmaxNTracks = cutsPSet.getParameter<bool>("applyCutsORmaxNTracks");
0082 }
0083
0084 /// destructor
0085 L3MuonIsolationProducer::~L3MuonIsolationProducer() {
0086   LogDebug("RecoMuon|L3MuonIsolationProducer") << " L3MuonIsolationProducer DTOR";
0087 }
0088
0089 void L3MuonIsolationProducer::produce(Event& event, const EventSetup& eventSetup) {
0090   std::string metname = "RecoMuon|L3MuonIsolationProducer";
0091
0092   LogDebug(metname) << " L3 Muon Isolation producing..."
0093                     << " BEGINING OF EVENT "
0094                     << "================================";
0095
0096   // Take the SA container
0097   LogTrace(metname) << " Taking the muons: " << theMuonCollectionLabel;
0098   Handle<TrackCollection> muons;
0099   event.getByToken(theMuonCollectionToken, muons);
0100
0101   auto depMap = std::make_unique<reco::IsoDepositMap>();
0102   auto isoMap = std::make_unique<edm::ValueMap<bool>>();
0103
0104   //
0105   // get Vetos and deposits
0106   //
0107   unsigned int nMuons = muons->size();
0108
0109   IsoDeposit::Vetos vetos(nMuons);
0110
0111   std::vector<IsoDeposit> deps(nMuons);
0112   std::vector<bool> isos(nMuons, false);
0113
0114   for (unsigned int i = 0; i < nMuons; i++) {
0115     TrackRef mu(muons, i);
0116     deps[i] = theExtractor->deposit(event, eventSetup, *mu);
0117     vetos[i] = deps[i].veto();
0118   }
0119
0120   //
0121   // add here additional vetos
0122   //
0123   //.....
0124
0125   //
0126   // actual cut step
0127   //
0128   for (unsigned int iMu = 0; iMu < nMuons; ++iMu) {
0129     const reco::Track* mu = &(*muons)[iMu];
0130
0131     const IsoDeposit& deposit = deps[iMu];
0132     LogTrace(metname) << deposit.print();
0133
0134     const Cuts::CutSpec& cut = theCuts(mu->eta());
0135     std::pair<double, int> sumAndCount = deposit.depositAndCountWithin(cut.conesize, vetos, theTrackPt_Min);
0136
0137     double value = sumAndCount.first;
0138     int count = sumAndCount.second;
0139
0140     bool result = (value < cut.threshold);
0141     if (theApplyCutsORmaxNTracks)
0142       result |= count <= theMaxNTracks;
0143     LogTrace(metname) << "deposit in cone: " << value << "with count " << count << " is isolated: " << result;
0144
0145     isos[iMu] = result;
0146   }
0147
0148   //
0149   // store
0150   //
0151   if (optOutputIsoDeposits) {
0152     reco::IsoDepositMap::Filler depFiller(*depMap);
0153     depFiller.insert(muons, deps.begin(), deps.end());
0154     depFiller.fill();
0155     event.put(std::move(depMap));
0156   }
0157   edm::ValueMap<bool>::Filler isoFiller(*isoMap);
0158   isoFiller.insert(muons, isos.begin(), isos.end());
0159   isoFiller.fill();
0160   event.put(std::move(isoMap));
0161
0162   LogTrace(metname) << " END OF EVENT "
0163                     << "================================";
0164 }