MuonIsolation/plugins/JetExtractor.cc

0001 #include "JetExtractor.h"
0002
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006
0007 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0008 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0009 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0010 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0011
0012 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0013
0014 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0015 #include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h"
0016 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0017
0018 #include "DataFormats/Math/interface/deltaR.h"
0019
0020 using namespace edm;
0021 using namespace std;
0022 using namespace reco;
0023 using namespace muonisolation;
0024 using reco::isodeposit::Direction;
0025
0026 JetExtractor::JetExtractor() {}
0027
0028 JetExtractor::JetExtractor(const ParameterSet& par, edm::ConsumesCollector&& iC)
0029     : theJetCollectionToken(iC.consumes<CaloJetCollection>(par.getParameter<edm::InputTag>("JetCollectionLabel"))),
0030       thePropagatorName(par.getParameter<std::string>("PropagatorName")),
0031       theFieldToken(iC.esConsumes()),
0032       theThreshold(par.getParameter<double>("Threshold")),
0033       theDR_Veto(par.getParameter<double>("DR_Veto")),
0034       theDR_Max(par.getParameter<double>("DR_Max")),
0035       theExcludeMuonVeto(par.getParameter<bool>("ExcludeMuonVeto")),
0036       theService(nullptr),
0037       theAssociator(nullptr),
0038       thePrintTimeReport(par.getUntrackedParameter<bool>("PrintTimeReport")) {
0039   ParameterSet serviceParameters = par.getParameter<ParameterSet>("ServiceParameters");
0040   theService = std::make_unique<MuonServiceProxy>(serviceParameters, edm::ConsumesCollector(iC));
0041
0042   //  theAssociatorParameters = new TrackAssociatorParameters(par.getParameter<edm::ParameterSet>("TrackAssociatorParameters"), iC_);
0043   theAssociatorParameters = std::make_unique<TrackAssociatorParameters>();
0044   theAssociatorParameters->loadParameters(par.getParameter<edm::ParameterSet>("TrackAssociatorParameters"), iC);
0045   theAssociator = std::make_unique<TrackDetectorAssociator>();
0046 }
0047
0048 JetExtractor::~JetExtractor() {}
0049
0050 void JetExtractor::fillVetos(const edm::Event& event, const edm::EventSetup& eventSetup, const TrackCollection& muons) {
0051   //   LogWarning("JetExtractor")
0052   //     <<"fillVetos does nothing now: IsoDeposit provides enough functionality\n"
0053   //     <<"to remove a deposit at/around given (eta, phi)";
0054 }
0055
0056 IsoDeposit JetExtractor::deposit(const Event& event, const EventSetup& eventSetup, const Track& muon) const {
0057   theService->update(eventSetup);
0058   theAssociator->setPropagator(&*(theService->propagator(thePropagatorName)));
0059
0060   typedef IsoDeposit::Veto Veto;
0061   IsoDeposit::Direction muonDir(muon.eta(), muon.phi());
0062
0063   IsoDeposit depJet(muonDir);
0064
0065   auto const& bField = eventSetup.getData(theFieldToken);
0066
0067   reco::TransientTrack tMuon(muon, &bField);
0068   FreeTrajectoryState iFTS = tMuon.initialFreeState();
0069   TrackDetMatchInfo mInfo = theAssociator->associate(event, eventSetup, iFTS, *theAssociatorParameters);
0070
0071   reco::isodeposit::Direction vetoDirection(mInfo.trkGlobPosAtHcal.eta(), mInfo.trkGlobPosAtHcal.phi());
0072   depJet.setVeto(Veto(vetoDirection, theDR_Veto));
0073
0074   edm::Handle<CaloJetCollection> caloJetsH;
0075   event.getByToken(theJetCollectionToken, caloJetsH);
0076
0077   //use calo towers
0078   CaloJetCollection::const_iterator jetCI = caloJetsH->begin();
0079   for (; jetCI != caloJetsH->end(); ++jetCI) {
0080     double deltar0 = reco::deltaR(muon, *jetCI);
0081     if (deltar0 > theDR_Max)
0082       continue;
0083     if (jetCI->et() < theThreshold)
0084       continue;
0085
0086     //should I make a separate config option for this?
0087     std::vector<CaloTowerPtr> jetConstituents = jetCI->getCaloConstituents();
0088
0089     std::vector<DetId>::const_iterator crossedCI = mInfo.crossedTowerIds.begin();
0090     std::vector<CaloTowerPtr>::const_iterator jetTowCI = jetConstituents.begin();
0091
0092     double sumEtExcluded = 0;
0093     for (; jetTowCI != jetConstituents.end(); ++jetTowCI) {
0094       bool isExcluded = false;
0095       double deltaRLoc = reco::deltaR(vetoDirection, *jetCI);
0096       if (deltaRLoc < theDR_Veto) {
0097         isExcluded = true;
0098       }
0099       for (; !isExcluded && crossedCI != mInfo.crossedTowerIds.end(); ++crossedCI) {
0100         if (crossedCI->rawId() == (*jetTowCI)->id().rawId()) {
0101           isExcluded = true;
0102         }
0103       }
0104       if (isExcluded)
0105         sumEtExcluded += (*jetTowCI)->et();
0106     }
0107     if (theExcludeMuonVeto) {
0108       if (jetCI->et() - sumEtExcluded < theThreshold)
0109         continue;
0110     }
0111
0112     double depositEt = jetCI->et();
0113     if (theExcludeMuonVeto)
0114       depositEt = depositEt - sumEtExcluded;
0115
0116     reco::isodeposit::Direction jetDir(jetCI->eta(), jetCI->phi());
0117     depJet.addDeposit(jetDir, depositEt);
0118   }
0119
0120   std::vector<const CaloTower*>::const_iterator crossedCI = mInfo.crossedTowers.begin();
0121   double muSumEt = 0;
0122   for (; crossedCI != mInfo.crossedTowers.end(); ++crossedCI) {
0123     muSumEt += (*crossedCI)->et();
0124   }
0125   depJet.addCandEnergy(muSumEt);
0126
0127   return depJet;
0128 }