Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-05-10 03:54:14

0001 #include "DQMOffline/Trigger/interface/HLTTauRefProducer.h"
0002 #include "TLorentzVector.h"
0003 // TAU includes
0004 #include "DataFormats/TauReco/interface/PFTau.h"
0005 #include "DataFormats/TauReco/interface/PFTauDiscriminator.h"
0006 #include "DataFormats/TauReco/interface/TauDiscriminatorContainer.h"
0007 // ELECTRON includes
0008 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0009 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0010 #include "DataFormats/EgammaCandidates/interface/Electron.h"
0011 #include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
0012 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0013 #include "DataFormats/TrackReco/interface/Track.h"
0014 // MUON includes
0015 #include "DataFormats/MuonReco/interface/Muon.h"
0016 #include "DataFormats/JetReco/interface/CaloJet.h"
0017 #include "TLorentzVector.h"
0018 #include "DataFormats/Math/interface/deltaR.h"
0019 //CaloTower includes
0020 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0021 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0022 #include "DataFormats/CaloTowers/interface/CaloTowerDefs.h"
0023 #include "Math/GenVector/VectorUtil.h"
0024 
0025 using namespace edm;
0026 using namespace reco;
0027 using namespace std;
0028 
0029 HLTTauRefProducer::HLTTauRefProducer(const edm::ParameterSet& iConfig) {
0030   //One Parameter Set per Collection
0031   {
0032     auto const& pfTau = iConfig.getUntrackedParameter<edm::ParameterSet>("PFTaus");
0033     PFTaus_ = consumes<reco::PFTauCollection>(pfTau.getUntrackedParameter<InputTag>("PFTauProducer"));
0034     auto discs = pfTau.getUntrackedParameter<vector<InputTag>>("PFTauDiscriminators");
0035     auto discConts = pfTau.getUntrackedParameter<vector<InputTag>>("PFTauDiscriminatorContainers");
0036     PFTauDisContWPs_ = pfTau.getUntrackedParameter<vector<std::string>>("PFTauDiscriminatorContainerWPs");
0037     if (discConts.size() != PFTauDisContWPs_.size())
0038       throw cms::Exception("Configuration") << "HLTTauRefProducer: Input parameters PFTauDiscriminatorContainers and "
0039                                                "PFTauDiscriminatorContainerWPs must have the same number of entries!\n";
0040     for (auto const& tag : discs) {
0041       PFTauDis_.push_back(consumes<reco::PFTauDiscriminator>(tag));
0042     }
0043     for (auto const& tag : discConts) {
0044       PFTauDisCont_.push_back(consumes<reco::TauDiscriminatorContainer>(tag));
0045     }
0046     doPFTaus_ = pfTau.getUntrackedParameter<bool>("doPFTaus", false);
0047     ptMinPFTau_ = pfTau.getUntrackedParameter<double>("ptMin", 15.);
0048     etaMinPFTau_ = pfTau.getUntrackedParameter<double>("etaMin", -2.5);
0049     etaMaxPFTau_ = pfTau.getUntrackedParameter<double>("etaMax", 2.5);
0050     phiMinPFTau_ = pfTau.getUntrackedParameter<double>("phiMin", -3.15);
0051     phiMaxPFTau_ = pfTau.getUntrackedParameter<double>("phiMax", 3.15);
0052   }
0053 
0054   {
0055     auto const& electrons = iConfig.getUntrackedParameter<edm::ParameterSet>("Electrons");
0056     Electrons_ = consumes<reco::GsfElectronCollection>(electrons.getUntrackedParameter<InputTag>("ElectronCollection"));
0057     doElectrons_ = electrons.getUntrackedParameter<bool>("doElectrons", false);
0058     e_ctfTrackCollectionSrc_ = electrons.getUntrackedParameter<InputTag>("TrackCollection");
0059     e_ctfTrackCollection_ = consumes<reco::TrackCollection>(e_ctfTrackCollectionSrc_);
0060     ptMinElectron_ = electrons.getUntrackedParameter<double>("ptMin", 15.);
0061     e_doTrackIso_ = electrons.getUntrackedParameter<bool>("doTrackIso", false);
0062     e_trackMinPt_ = electrons.getUntrackedParameter<double>("ptMinTrack", 1.5);
0063     e_lipCut_ = electrons.getUntrackedParameter<double>("lipMinTrack", 1.5);
0064     e_minIsoDR_ = electrons.getUntrackedParameter<double>("InnerConeDR", 0.02);
0065     e_maxIsoDR_ = electrons.getUntrackedParameter<double>("OuterConeDR", 0.6);
0066     e_isoMaxSumPt_ = electrons.getUntrackedParameter<double>("MaxIsoVar", 0.02);
0067   }
0068 
0069   {
0070     auto const& muons = iConfig.getUntrackedParameter<edm::ParameterSet>("Muons");
0071     Muons_ = consumes<reco::MuonCollection>(muons.getUntrackedParameter<InputTag>("MuonCollection"));
0072     doMuons_ = muons.getUntrackedParameter<bool>("doMuons", false);
0073     ptMinMuon_ = muons.getUntrackedParameter<double>("ptMin", 15.);
0074   }
0075 
0076   {
0077     auto const& jets = iConfig.getUntrackedParameter<edm::ParameterSet>("Jets");
0078     Jets_ = consumes<reco::CaloJetCollection>(jets.getUntrackedParameter<InputTag>("JetCollection"));
0079     doJets_ = jets.getUntrackedParameter<bool>("doJets");
0080     ptMinJet_ = jets.getUntrackedParameter<double>("etMin");
0081   }
0082 
0083   {
0084     auto const& towers = iConfig.getUntrackedParameter<edm::ParameterSet>("Towers");
0085     Towers_ = consumes<CaloTowerCollection>(towers.getUntrackedParameter<InputTag>("TowerCollection"));
0086     doTowers_ = towers.getUntrackedParameter<bool>("doTowers");
0087     ptMinTower_ = towers.getUntrackedParameter<double>("etMin");
0088     towerIsol_ = towers.getUntrackedParameter<double>("towerIsolation");
0089   }
0090 
0091   {
0092     auto const& photons = iConfig.getUntrackedParameter<edm::ParameterSet>("Photons");
0093     Photons_ = consumes<reco::PhotonCollection>(photons.getUntrackedParameter<InputTag>("PhotonCollection"));
0094     doPhotons_ = photons.getUntrackedParameter<bool>("doPhotons");
0095     ptMinPhoton_ = photons.getUntrackedParameter<double>("etMin");
0096     photonEcalIso_ = photons.getUntrackedParameter<double>("ECALIso");
0097   }
0098 
0099   {
0100     auto const& met = iConfig.getUntrackedParameter<edm::ParameterSet>("MET");
0101     MET_ = consumes<reco::CaloMETCollection>(met.getUntrackedParameter<InputTag>("METCollection"));
0102     doMET_ = met.getUntrackedParameter<bool>("doMET", false);
0103     ptMinMET_ = met.getUntrackedParameter<double>("ptMin", 15.);
0104   }
0105 
0106   etaMin_ = iConfig.getUntrackedParameter<double>("EtaMin", -2.5);
0107   etaMax_ = iConfig.getUntrackedParameter<double>("EtaMax", 2.5);
0108   phiMin_ = iConfig.getUntrackedParameter<double>("PhiMin", -3.15);
0109   phiMax_ = iConfig.getUntrackedParameter<double>("PhiMax", 3.15);
0110 
0111   //recoCollections
0112   produces<LorentzVectorCollection>("PFTaus");
0113   produces<LorentzVectorCollection>("Electrons");
0114   produces<LorentzVectorCollection>("Muons");
0115   produces<LorentzVectorCollection>("Jets");
0116   produces<LorentzVectorCollection>("Photons");
0117   produces<LorentzVectorCollection>("Towers");
0118   produces<LorentzVectorCollection>("MET");
0119 }
0120 
0121 void HLTTauRefProducer::produce(edm::StreamID iID, edm::Event& iEvent, edm::EventSetup const&) const {
0122   if (doPFTaus_)
0123     doPFTaus(iID, iEvent);
0124   if (doElectrons_)
0125     doElectrons(iEvent);
0126   if (doMuons_)
0127     doMuons(iEvent);
0128   if (doJets_)
0129     doJets(iEvent);
0130   if (doPhotons_)
0131     doPhotons(iEvent);
0132   if (doTowers_)
0133     doTowers(iEvent);
0134   if (doMET_)
0135     doMET(iEvent);
0136 }
0137 
0138 void HLTTauRefProducer::doPFTaus(edm::StreamID iID, edm::Event& iEvent) const {
0139   auto product_PFTaus = make_unique<LorentzVectorCollection>();
0140 
0141   edm::Handle<PFTauCollection> pftaus;
0142   if (iEvent.getByToken(PFTaus_, pftaus)) {
0143     // Retrieve ID container indices if config history changes, in particular for the first event.
0144     if (streamCache(iID)->first != iEvent.processHistoryID()) {
0145       streamCache(iID)->first = iEvent.processHistoryID();
0146       streamCache(iID)->second.resize(PFTauDisContWPs_.size());
0147       for (size_t i = 0; i < PFTauDisCont_.size(); ++i) {
0148         auto const aHandle = iEvent.getHandle(PFTauDisCont_[i]);
0149         auto const aProv = aHandle.provenance();
0150         if (aProv == nullptr)
0151           aHandle.whyFailed()->raise();
0152         const auto& psetsFromProvenance = edm::parameterSet(aProv->stable(), iEvent.processHistory());
0153         if (psetsFromProvenance.exists("workingPoints")) {
0154           auto const idlist = psetsFromProvenance.getParameter<std::vector<std::string>>("workingPoints");
0155           bool found = false;
0156           for (size_t j = 0; j < idlist.size(); ++j) {
0157             if (PFTauDisContWPs_[i] == idlist[j]) {
0158               found = true;
0159               streamCache(iID)->second[i] = j;
0160             }
0161           }
0162           if (!found)
0163             throw cms::Exception("Configuration")
0164                 << "HLTTauRefProducer: Requested working point '" << PFTauDisContWPs_[i] << "' not found!\n";
0165         } else if (psetsFromProvenance.exists("IDWPdefinitions")) {
0166           auto const idlist = psetsFromProvenance.getParameter<std::vector<edm::ParameterSet>>("IDWPdefinitions");
0167           bool found = false;
0168           for (size_t j = 0; j < idlist.size(); ++j) {
0169             if (PFTauDisContWPs_[i] == idlist[j].getParameter<std::string>("IDname")) {
0170               found = true;
0171               streamCache(iID)->second[i] = j;
0172             }
0173           }
0174           if (!found)
0175             throw cms::Exception("Configuration")
0176                 << "HLTTauRefProducer: Requested working point '" << PFTauDisContWPs_[i] << "' not found!\n";
0177         } else
0178           throw cms::Exception("Configuration")
0179               << "HLTTauRefProducer: No suitable ID list found in provenace config!\n";
0180       }
0181     }
0182     for (unsigned int i = 0; i < pftaus->size(); ++i) {
0183       auto const& pftau = (*pftaus)[i];
0184       if (pftau.pt() > ptMinPFTau_ && pftau.eta() > etaMinPFTau_ && pftau.eta() < etaMaxPFTau_ &&
0185           pftau.phi() > phiMinPFTau_ && pftau.phi() < phiMaxPFTau_) {
0186         reco::PFTauRef thePFTau{pftaus, i};
0187         bool passAll{true};
0188 
0189         for (auto const& token : PFTauDis_) {
0190           edm::Handle<reco::PFTauDiscriminator> pftaudis;
0191           if (iEvent.getByToken(token, pftaudis)) {
0192             if ((*pftaudis)[thePFTau] < 0.5) {
0193               passAll = false;
0194               break;
0195             }
0196           } else {
0197             passAll = false;
0198             break;
0199           }
0200         }
0201 
0202         int idx = 0;
0203         for (auto const& token : PFTauDisCont_) {
0204           edm::Handle<reco::TauDiscriminatorContainer> pftaudis;
0205           if (iEvent.getByToken(token, pftaudis)) {
0206             //WP vector not filled if prediscriminator in RecoTauDiscriminator failed.
0207             if ((*pftaudis)[thePFTau].workingPoints.empty() ||
0208                 !(*pftaudis)[thePFTau].workingPoints.at(streamCache(iID)->second[idx])) {
0209               passAll = false;
0210               break;
0211             }
0212           } else {
0213             passAll = false;
0214             break;
0215           }
0216           idx++;
0217         }
0218         if (passAll) {
0219           product_PFTaus->emplace_back(pftau.px(), pftau.py(), pftau.pz(), pftau.energy());
0220         }
0221       }
0222     }
0223   }
0224   iEvent.put(std::move(product_PFTaus), "PFTaus");
0225 }
0226 
0227 void HLTTauRefProducer::doElectrons(edm::Event& iEvent) const {
0228   auto product_Electrons = make_unique<LorentzVectorCollection>();
0229 
0230   edm::Handle<reco::TrackCollection> pCtfTracks;
0231   if (!iEvent.getByToken(e_ctfTrackCollection_, pCtfTracks)) {
0232     edm::LogInfo("") << "Error! Can't get " << e_ctfTrackCollectionSrc_.label() << " by label. ";
0233     iEvent.put(std::move(product_Electrons), "Electrons");
0234     return;
0235   }
0236 
0237   edm::Handle<GsfElectronCollection> electrons;
0238   if (iEvent.getByToken(Electrons_, electrons)) {
0239     for (size_t i = 0; i < electrons->size(); ++i) {
0240       edm::Ref<reco::GsfElectronCollection> electronRef(electrons, i);
0241       auto const& electron = (*electrons)[i];
0242       if (electron.pt() > ptMinElectron_ && fabs(electron.eta()) < etaMax_) {
0243         if (e_doTrackIso_) {
0244           double sum_of_pt_ele{};
0245           for (auto const& tr : *pCtfTracks) {
0246             double const lip{electron.gsfTrack()->dz() - tr.dz()};
0247             if (tr.pt() > e_trackMinPt_ && fabs(lip) < e_lipCut_) {
0248               double dphi{fabs(tr.phi() - electron.trackMomentumAtVtx().phi())};
0249               if (dphi > acos(-1.)) {
0250                 dphi = 2 * acos(-1.) - dphi;
0251               }
0252               double const deta{fabs(tr.eta() - electron.trackMomentumAtVtx().eta())};
0253               double const dr_ctf_ele{sqrt(deta * deta + dphi * dphi)};
0254               if ((dr_ctf_ele > e_minIsoDR_) && (dr_ctf_ele < e_maxIsoDR_)) {
0255                 double const cft_pt_2{tr.pt() * tr.pt()};
0256                 sum_of_pt_ele += cft_pt_2;
0257               }
0258             }
0259           }
0260           double const isolation_value_ele{sum_of_pt_ele /
0261                                            (electron.trackMomentumAtVtx().Rho() * electron.trackMomentumAtVtx().Rho())};
0262           if (isolation_value_ele < e_isoMaxSumPt_) {
0263             product_Electrons->emplace_back(electron.px(), electron.py(), electron.pz(), electron.energy());
0264           }
0265         } else {
0266           product_Electrons->emplace_back(electron.px(), electron.py(), electron.pz(), electron.energy());
0267         }
0268       }
0269     }
0270   }
0271   iEvent.put(std::move(product_Electrons), "Electrons");
0272 }
0273 
0274 void HLTTauRefProducer::doMuons(edm::Event& iEvent) const {
0275   auto product_Muons = make_unique<LorentzVectorCollection>();
0276 
0277   edm::Handle<MuonCollection> muons;
0278   if (iEvent.getByToken(Muons_, muons)) {
0279     for (auto const& muon : *muons) {
0280       if (muon.pt() > ptMinMuon_ && muon.eta() > etaMin_ && muon.eta() < etaMax_ && muon.phi() > phiMin_ &&
0281           muon.phi() < phiMax_) {
0282         product_Muons->emplace_back(muon.px(), muon.py(), muon.pz(), muon.energy());
0283       }
0284     }
0285   }
0286   iEvent.put(std::move(product_Muons), "Muons");
0287 }
0288 
0289 void HLTTauRefProducer::doJets(edm::Event& iEvent) const {
0290   auto product_Jets = make_unique<LorentzVectorCollection>();
0291 
0292   edm::Handle<CaloJetCollection> jets;
0293   if (iEvent.getByToken(Jets_, jets)) {
0294     for (auto const& jet : *jets) {
0295       if (jet.et() > ptMinJet_ && jet.eta() > etaMin_ && jet.eta() < etaMax_ && jet.phi() > phiMin_ &&
0296           jet.phi() < phiMax_) {
0297         product_Jets->emplace_back(jet.px(), jet.py(), jet.pz(), jet.energy());
0298       }
0299     }
0300   }
0301   iEvent.put(std::move(product_Jets), "Jets");
0302 }
0303 
0304 void HLTTauRefProducer::doTowers(edm::Event& iEvent) const {
0305   auto product_Towers = make_unique<LorentzVectorCollection>();
0306 
0307   edm::Handle<CaloTowerCollection> towers;
0308   if (iEvent.getByToken(Towers_, towers)) {
0309     for (auto const& tower1 : *towers) {
0310       if (tower1.pt() > ptMinTower_ && tower1.eta() > etaMin_ && tower1.eta() < etaMax_ && tower1.phi() > phiMin_ &&
0311           tower1.phi() < phiMax_) {
0312         //calculate isolation
0313         double isolET{};
0314         for (auto const& tower2 : *towers) {
0315           if (ROOT::Math::VectorUtil::DeltaR(tower1.p4(), tower2.p4()) < 0.5) {
0316             isolET += tower2.pt();
0317           }
0318           isolET -= tower1.pt();
0319         }
0320         if (isolET < towerIsol_) {
0321           product_Towers->emplace_back(tower1.px(), tower1.py(), tower1.pz(), tower1.energy());
0322         }
0323       }
0324     }
0325   }
0326   iEvent.put(std::move(product_Towers), "Towers");
0327 }
0328 
0329 void HLTTauRefProducer::doPhotons(edm::Event& iEvent) const {
0330   auto product_Gammas = make_unique<LorentzVectorCollection>();
0331 
0332   edm::Handle<PhotonCollection> photons;
0333   if (iEvent.getByToken(Photons_, photons)) {
0334     for (auto const& photon : *photons) {
0335       if (photon.ecalRecHitSumEtConeDR04() < photonEcalIso_ && photon.et() > ptMinPhoton_ && photon.eta() > etaMin_ &&
0336           photon.eta() < etaMax_ && photon.phi() > phiMin_ && photon.phi() < phiMax_) {
0337         product_Gammas->emplace_back(photon.px(), photon.py(), photon.pz(), photon.energy());
0338       }
0339     }
0340   }
0341   iEvent.put(std::move(product_Gammas), "Photons");
0342 }
0343 
0344 void HLTTauRefProducer::doMET(edm::Event& iEvent) const {
0345   auto product_MET = make_unique<LorentzVectorCollection>();
0346 
0347   edm::Handle<reco::CaloMETCollection> met;
0348   if (iEvent.getByToken(MET_, met) && !met->empty()) {
0349     auto const& metMom = met->front().p4();
0350     product_MET->emplace_back(metMom.Px(), metMom.Py(), 0, metMom.Pt());
0351   }
0352   iEvent.put(std::move(product_MET), "MET");
0353 }
0354 
0355 #include "FWCore/Framework/interface/MakerMacros.h"
0356 DEFINE_FWK_MODULE(HLTTauRefProducer);