L3MuonIsolationProducer/plugins/L3MuonCombinedRelativeIsolationProducer.cc

0001 #include "RecoMuon/L3MuonIsolationProducer/plugins/L3MuonCombinedRelativeIsolationProducer.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 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0017 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0018
0019 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0020 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
0021
0022 #include "RecoMuon/MuonIsolation/interface/Range.h"
0023 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
0024
0025 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractor.h"
0026 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractorFactory.h"
0027
0028 #include "RecoMuon/L3MuonIsolationProducer/interface/L3NominalEfficiencyConfigurator.h"
0029
0030 #include <string>
0031
0032 using namespace edm;
0033 using namespace std;
0034 using namespace reco;
0035 using namespace muonisolation;
0036
0037 /// constructor with config
0038 L3MuonCombinedRelativeIsolationProducer::L3MuonCombinedRelativeIsolationProducer(const ParameterSet& par)
0039     : theConfig(par),
0040       theMuonCollectionLabel(par.getParameter<InputTag>("inputMuonCollection")),
0041       optOutputIsoDeposits(par.getParameter<bool>("OutputMuIsoDeposits")),
0042       useCaloIso(par.existsAs<bool>("UseCaloIso") ? par.getParameter<bool>("UseCaloIso") : true),
0043       useRhoCorrectedCaloDeps(par.existsAs<bool>("UseRhoCorrectedCaloDeposits")
0044                                   ? par.getParameter<bool>("UseRhoCorrectedCaloDeposits")
0045                                   : false),
0046       theCaloDepsLabel(par.existsAs<InputTag>("CaloDepositsLabel") ? par.getParameter<InputTag>("CaloDepositsLabel")
0047                                                                    : InputTag("hltL3CaloMuonCorrectedIsolations")),
0048       theTrackPt_Min(-1),
0049       printDebug(par.getParameter<bool>("printDebug")) {
0050   LogDebug("RecoMuon|L3MuonCombinedRelativeIsolationProducer") << " L3MuonCombinedRelativeIsolationProducer CTOR";
0051
0052   theMuonCollectionToken = consumes<RecoChargedCandidateCollection>(theMuonCollectionLabel);
0053   theCaloDepsToken = consumes<edm::ValueMap<float>>(theCaloDepsLabel);
0054
0055   if (optOutputIsoDeposits) {
0056     produces<reco::IsoDepositMap>("trkIsoDeposits");
0057     if (useRhoCorrectedCaloDeps == false)  // otherwise, calo deposits have been previously computed
0058       produces<reco::IsoDepositMap>("caloIsoDeposits");
0059     //produces<std::vector<double> >("combinedRelativeIsoDeposits");
0060     produces<edm::ValueMap<double>>("combinedRelativeIsoDeposits");
0061   }
0062   produces<edm::ValueMap<bool>>();
0063
0064   //
0065   // Extractor
0066   //
0067   // Calorimeters (ONLY if not previously computed)
0068   //
0069   if (useCaloIso && (useRhoCorrectedCaloDeps == false)) {
0070     edm::ParameterSet caloExtractorPSet = theConfig.getParameter<edm::ParameterSet>("CaloExtractorPSet");
0071
0072     std::string caloExtractorName = caloExtractorPSet.getParameter<std::string>("ComponentName");
0073     caloExtractor = std::unique_ptr<reco::isodeposit::IsoDepositExtractor>{
0074         IsoDepositExtractorFactory::get()->create(caloExtractorName, caloExtractorPSet, consumesCollector())};
0075     //std::string caloDepositType = caloExtractorPSet.getUntrackedParameter<std::string>("DepositLabel"); // N.B. Not used in the following!
0076   }
0077
0078   // Tracker
0079   //
0080   edm::ParameterSet trkExtractorPSet = theConfig.getParameter<edm::ParameterSet>("TrkExtractorPSet");
0081
0082   theTrackPt_Min = theConfig.getParameter<double>("TrackPt_Min");
0083   std::string trkExtractorName = trkExtractorPSet.getParameter<std::string>("ComponentName");
0084   trkExtractor = std::unique_ptr<reco::isodeposit::IsoDepositExtractor>{
0085       IsoDepositExtractorFactory::get()->create(trkExtractorName, trkExtractorPSet, consumesCollector())};
0086   //std::string trkDepositType = trkExtractorPSet.getUntrackedParameter<std::string>("DepositLabel"); // N.B. Not used in the following!
0087
0088   //
0089   // Cuts for track isolation
0090   //
0091   edm::ParameterSet cutsPSet = theConfig.getParameter<edm::ParameterSet>("CutsPSet");
0092   std::string cutsName = cutsPSet.getParameter<std::string>("ComponentName");
0093   if (cutsName == "SimpleCuts") {
0094     theCuts = Cuts(cutsPSet);
0095   } else if (
0096       //        (cutsName== "L3NominalEfficiencyCuts_PXLS" && depositType=="PXLS")
0097       //     || (cutsName== "L3NominalEfficiencyCuts_TRKS" && depositType=="TRKS")
0098       //! test cutsName only. The depositType is informational only (has not been used so far) [VK]
0099       (cutsName == "L3NominalEfficiencyCuts_PXLS") || (cutsName == "L3NominalEfficiencyCuts_TRKS")) {
0100     theCuts = L3NominalEfficiencyConfigurator(cutsPSet).cuts();
0101   } else {
0102     LogError("L3MuonCombinedRelativeIsolationProducer::beginJob") << "cutsName: " << cutsPSet << " is not recognized:"
0103                                                                   << " theCuts not set!";
0104   }
0105   LogTrace("") << theCuts.print();
0106
0107   // (kludge) additional cut on the number of tracks
0108   theMaxNTracks = cutsPSet.getParameter<int>("maxNTracks");
0109   theApplyCutsORmaxNTracks = cutsPSet.getParameter<bool>("applyCutsORmaxNTracks");
0110 }
0111
0112 /// destructor
0113 L3MuonCombinedRelativeIsolationProducer::~L3MuonCombinedRelativeIsolationProducer() {
0114   LogDebug("RecoMuon|L3MuonCombinedRelativeIsolationProducer") << " L3MuonCombinedRelativeIsolationProducer DTOR";
0115 }
0116
0117 /// ParameterSet descriptions
0118 void L3MuonCombinedRelativeIsolationProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0119   edm::ParameterSetDescription desc;
0120   desc.add<bool>("UseRhoCorrectedCaloDeposits", false);
0121   desc.add<bool>("UseCaloIso", true);
0122   desc.add<edm::InputTag>("CaloDepositsLabel", edm::InputTag("hltL3CaloMuonCorrectedIsolations"));
0123   desc.add<edm::InputTag>("inputMuonCollection", edm::InputTag("hltL3MuonCandidates"));
0124   desc.add<bool>("OutputMuIsoDeposits", true);
0125   desc.add<double>("TrackPt_Min", -1.0);
0126   desc.add<bool>("printDebug", false);
0127   edm::ParameterSetDescription cutsPSet;
0128   {
0129     cutsPSet.add<std::vector<double>>("ConeSizes", std::vector<double>(1, 0.24));
0130     cutsPSet.add<std::string>("ComponentName", "SimpleCuts");
0131     cutsPSet.add<std::vector<double>>("Thresholds", std::vector<double>(1, 0.1));
0132     cutsPSet.add<int>("maxNTracks", -1);
0133     cutsPSet.add<std::vector<double>>("EtaBounds", std::vector<double>(1, 2.411));
0134     cutsPSet.add<bool>("applyCutsORmaxNTracks", false);
0135   }
0136   desc.add<edm::ParameterSetDescription>("CutsPSet", cutsPSet);
0137   edm::ParameterSetDescription trkExtractorPSet;
0138   {
0139     trkExtractorPSet.add<double>("Chi2Prob_Min", -1.0);
0140     trkExtractorPSet.add<double>("Chi2Ndof_Max", 1.0E64);
0141     trkExtractorPSet.add<double>("Diff_z", 0.2);
0142     trkExtractorPSet.add<edm::InputTag>("inputTrackCollection", edm::InputTag("hltPixelTracks"));
0143     trkExtractorPSet.add<double>("ReferenceRadius", 6.0);
0144     trkExtractorPSet.add<edm::InputTag>("BeamSpotLabel", edm::InputTag("hltOnlineBeamSpot"));
0145     trkExtractorPSet.add<std::string>("ComponentName", "PixelTrackExtractor");
0146     trkExtractorPSet.add<double>("DR_Max", 0.24);
0147     trkExtractorPSet.add<double>("Diff_r", 0.1);
0148     trkExtractorPSet.add<bool>("VetoLeadingTrack", true);
0149     trkExtractorPSet.add<double>("DR_VetoPt", 0.025);
0150     trkExtractorPSet.add<double>("DR_Veto", 0.01);
0151     trkExtractorPSet.add<unsigned int>("NHits_Min", 0);
0152     trkExtractorPSet.add<double>("Pt_Min", -1.0);
0153     trkExtractorPSet.addUntracked<std::string>("DepositLabel", "PXLS");
0154     trkExtractorPSet.add<std::string>("BeamlineOption", "BeamSpotFromEvent");
0155     trkExtractorPSet.add<bool>("PropagateTracksToRadius", true);
0156     trkExtractorPSet.add<double>("PtVeto_Min", 2.0);
0157   }
0158   desc.add<edm::ParameterSetDescription>("TrkExtractorPSet", trkExtractorPSet);
0159   edm::ParameterSetDescription caloExtractorPSet;
0160   {
0161     caloExtractorPSet.add<double>("DR_Veto_H", 0.1);
0162     caloExtractorPSet.add<bool>("Vertex_Constraint_Z", false);
0163     caloExtractorPSet.add<double>("Threshold_H", 0.5);
0164     caloExtractorPSet.add<std::string>("ComponentName", "CaloExtractor");
0165     caloExtractorPSet.add<double>("Threshold_E", 0.2);
0166     caloExtractorPSet.add<double>("DR_Max", 0.24);
0167     caloExtractorPSet.add<double>("DR_Veto_E", 0.07);
0168     caloExtractorPSet.add<double>("Weight_E", 1.5);
0169     caloExtractorPSet.add<bool>("Vertex_Constraint_XY", false);
0170     caloExtractorPSet.addUntracked<std::string>("DepositLabel", "EcalPlusHcal");
0171     caloExtractorPSet.add<edm::InputTag>("CaloTowerCollectionLabel", edm::InputTag("hltTowerMakerForMuons"));
0172     caloExtractorPSet.add<double>("Weight_H", 1.0);
0173   }
0174   desc.add<edm::ParameterSetDescription>("CaloExtractorPSet", caloExtractorPSet);
0175   descriptions.add("hltL3MuonIsolations", desc);
0176 }
0177
0178 void L3MuonCombinedRelativeIsolationProducer::produce(Event& event, const EventSetup& eventSetup) {
0179   std::string metname = "RecoMuon|L3MuonCombinedRelativeIsolationProducer";
0180
0181   if (printDebug)
0182     std::cout << " L3 Muon Isolation producing..."
0183               << " BEGINING OF EVENT "
0184               << "================================" << std::endl;
0185
0186   // Take the SA container
0187   if (printDebug)
0188     std::cout << " Taking the muons: " << theMuonCollectionLabel << std::endl;
0189   Handle<RecoChargedCandidateCollection> muons;
0190   event.getByToken(theMuonCollectionToken, muons);
0191
0192   // Take calo deposits with rho corrections (ONLY if previously computed)
0193   Handle<edm::ValueMap<float>> caloDepWithCorrMap;
0194   if (useRhoCorrectedCaloDeps)
0195     event.getByToken(theCaloDepsToken, caloDepWithCorrMap);
0196
0197   auto caloDepMap = std::make_unique<reco::IsoDepositMap>();
0198   auto trkDepMap = std::make_unique<reco::IsoDepositMap>();
0199
0200   auto comboIsoDepMap = std::make_unique<edm::ValueMap<bool>>();
0201
0202   //auto combinedRelativeDeps = std::make_unique<std::vector<double>>();
0203   auto combinedRelativeDepMap = std::make_unique<edm::ValueMap<double>>();
0204
0205   //
0206   // get Vetos and deposits
0207   //
0208   unsigned int nMuons = muons->size();
0209
0210   IsoDeposit::Vetos trkVetos(nMuons);
0211   std::vector<IsoDeposit> trkDeps(nMuons);
0212
0213   // IsoDeposit::Vetos caloVetos(nMuons);
0214   // std::vector<IsoDeposit> caloDeps(nMuons);
0215   // std::vector<float> caloCorrDeps(nMuons, 0.);  // if calo deposits with corrections available
0216
0217   IsoDeposit::Vetos caloVetos;
0218   std::vector<IsoDeposit> caloDeps;
0219   std::vector<float> caloCorrDeps;  // if calo deposits with corrections available
0220
0221   if (useCaloIso && useRhoCorrectedCaloDeps) {
0222     caloCorrDeps.resize(nMuons, 0.);
0223   } else if (useCaloIso) {
0224     caloVetos.resize(nMuons);
0225     caloDeps.resize(nMuons);
0226   }
0227
0228   std::vector<double> combinedRelativeDeps(nMuons, 0.);
0229   std::vector<bool> combinedRelativeIsos(nMuons, false);
0230
0231   for (unsigned int i = 0; i < nMuons; i++) {
0232     RecoChargedCandidateRef candref(muons, i);
0233     TrackRef mu = candref->track();
0234
0235     trkDeps[i] = trkExtractor->deposit(event, eventSetup, *mu);
0236     trkVetos[i] = trkDeps[i].veto();
0237
0238     if (useCaloIso && useRhoCorrectedCaloDeps) {
0239       caloCorrDeps[i] = (*caloDepWithCorrMap)[candref];
0240     } else if (useCaloIso) {
0241       caloDeps[i] = caloExtractor->deposit(event, eventSetup, *mu);
0242       caloVetos[i] = caloDeps[i].veto();
0243     }
0244   }
0245
0246   //
0247   // add here additional vetos
0248   //
0249   //.....
0250
0251   //
0252   // actual cut step
0253   //
0254
0255   if (printDebug)
0256     std::cout << "Looping over deposits...." << std::endl;
0257   for (unsigned int iMu = 0; iMu < nMuons; ++iMu) {
0258     if (printDebug)
0259       std::cout << "Muon number = " << iMu << std::endl;
0260     TrackRef mu = (*muons)[iMu].track();
0261
0262     // cuts
0263     const Cuts::CutSpec& cut = theCuts(mu->eta());
0264
0265     if (printDebug)
0266       std::cout << "CUTDEBUG: Muon eta = " << mu->eta() << std::endl
0267                 << "CUTDEBUG: Muon pt  = " << mu->pt() << std::endl
0268                 << "CUTDEBUG: minEta   = " << cut.etaRange.min() << std::endl
0269                 << "CUTDEBUG: maxEta   = " << cut.etaRange.max() << std::endl
0270                 << "CUTDEBUG: consize  = " << cut.conesize << std::endl
0271                 << "CUTDEBUG: thresho  = " << cut.threshold << std::endl;
0272
0273     const IsoDeposit& trkDeposit = trkDeps[iMu];
0274     if (printDebug)
0275       std::cout << trkDeposit.print();
0276     std::pair<double, int> trkIsoSumAndCount = trkDeposit.depositAndCountWithin(cut.conesize, trkVetos, theTrackPt_Min);
0277
0278     double caloIsoSum = 0.;
0279     if (useCaloIso && useRhoCorrectedCaloDeps) {
0280       caloIsoSum = caloCorrDeps[iMu];
0281       if (caloIsoSum < 0.)
0282         caloIsoSum = 0.;
0283       if (printDebug)
0284         std::cout << "Rho-corrected calo deposit (min. 0) = " << caloIsoSum << std::endl;
0285     } else if (useCaloIso) {
0286       const IsoDeposit& caloDeposit = caloDeps[iMu];
0287       if (printDebug)
0288         std::cout << caloDeposit.print();
0289       caloIsoSum = caloDeposit.depositWithin(cut.conesize, caloVetos);
0290     }
0291
0292     double trkIsoSum = trkIsoSumAndCount.first;
0293     int count = trkIsoSumAndCount.second;
0294
0295     double muPt = mu->pt();
0296     if (muPt < 1.)
0297       muPt = 1.;
0298     double combinedRelativeDeposit = ((trkIsoSum + caloIsoSum) / muPt);
0299     bool result = (combinedRelativeDeposit < cut.threshold);
0300     if (theApplyCutsORmaxNTracks)
0301       result |= count <= theMaxNTracks;
0302     if (printDebug)
0303       std::cout << "  trk dep in cone:  " << trkIsoSum << "  with count " << count << std::endl
0304                 << "  calo dep in cone: " << caloIsoSum << std::endl
0305                 << "  muPt: " << muPt << std::endl
0306                 << "  relIso:  " << combinedRelativeDeposit << std::endl
0307                 << "  is isolated: " << result << std::endl;
0308
0309     combinedRelativeIsos[iMu] = result;
0310     //combinedRelativeDeps->push_back(combinedRelativeDeposit);
0311     combinedRelativeDeps[iMu] = combinedRelativeDeposit;
0312   }
0313
0314   //
0315   // store
0316   //
0317   if (optOutputIsoDeposits) {
0318     reco::IsoDepositMap::Filler depFillerTrk(*trkDepMap);
0319     depFillerTrk.insert(muons, trkDeps.begin(), trkDeps.end());
0320     depFillerTrk.fill();
0321     event.put(std::move(trkDepMap), "trkIsoDeposits");
0322
0323     if (useCaloIso && (useRhoCorrectedCaloDeps == false)) {
0324       reco::IsoDepositMap::Filler depFillerCalo(*caloDepMap);
0325       depFillerCalo.insert(muons, caloDeps.begin(), caloDeps.end());
0326       depFillerCalo.fill();
0327       event.put(std::move(caloDepMap), "caloIsoDeposits");
0328     }
0329
0330     //event.put(std::move(combinedRelativeDeps, "combinedRelativeIsoDeposits");
0331     edm::ValueMap<double>::Filler depFillerCombRel(*combinedRelativeDepMap);
0332     depFillerCombRel.insert(muons, combinedRelativeDeps.begin(), combinedRelativeDeps.end());
0333     depFillerCombRel.fill();
0334     event.put(std::move(combinedRelativeDepMap), "combinedRelativeIsoDeposits");
0335   }
0336   edm::ValueMap<bool>::Filler isoFiller(*comboIsoDepMap);
0337   isoFiller.insert(muons, combinedRelativeIsos.begin(), combinedRelativeIsos.end());
0338   isoFiller.fill();
0339   event.put(std::move(comboIsoDepMap));
0340
0341   if (printDebug)
0342     std::cout << " END OF EVENT "
0343               << "================================";
0344 }
0345
0346 #include "FWCore/Framework/interface/MakerMacros.h"
0347 DEFINE_FWK_MODULE(L3MuonCombinedRelativeIsolationProducer);