File indexing completed on 2024-04-06 12:21:35
0001 #include "L1Trigger/Phase2L1Taus/interface/L1HPSPFTauBuilder.h"
0002 #include "FWCore/Utilities/interface/Exception.h" // cms::Exception
0003 #include "DataFormats/Math/interface/deltaR.h" // reco::deltaR
0004 #include <regex> // sd::regex_replace
0005 #include <TMath.h> // TMath::Pi()
0006 #include <string> // std::string
0007 #include <algorithm> // std::max(), std::sort()
0008 #include <cmath> // std::fabs
0009
0010 namespace {
0011 std::string getSignalConeSizeFormula(const edm::ParameterSet& cfg) {
0012 return std::regex_replace(cfg.getParameter<std::string>("signalConeSize"), std::regex("pt"), "x");
0013 }
0014 }
0015
0016 L1HPSPFTauBuilder::L1HPSPFTauBuilder(const edm::ParameterSet& cfg)
0017 : signalConeSizeFormula_(getSignalConeSizeFormula(cfg)),
0018 minSignalConeSize_(cfg.getParameter<double>("minSignalConeSize")),
0019 maxSignalConeSize_(cfg.getParameter<double>("maxSignalConeSize")),
0020 useStrips_(cfg.getParameter<bool>("useStrips")),
0021 stripSizeEta_(cfg.getParameter<double>("stripSizeEta")),
0022 stripSizePhi_(cfg.getParameter<double>("stripSizePhi")),
0023 isolationConeSize_(cfg.getParameter<double>("isolationConeSize")),
0024 debug_(cfg.getUntrackedParameter<bool>("debug", false)) {
0025 assert(maxSignalConeSize_ >= minSignalConeSize_);
0026
0027 isolationConeSize2_ = isolationConeSize_ * isolationConeSize_;
0028
0029 if (debug_) {
0030 std::cout << "setting Quality cuts for signal PFCands:" << std::endl;
0031 }
0032 edm::ParameterSet cfg_signalQualityCuts = cfg.getParameter<edm::ParameterSet>("signalQualityCuts");
0033 signalQualityCutsDzCutDisabled_ = readL1PFTauQualityCuts(cfg_signalQualityCuts, "disabled", debug_);
0034 signalQualityCutsDzCutEnabledPrimary_ = readL1PFTauQualityCuts(cfg_signalQualityCuts, "enabled_primary", debug_);
0035 if (debug_) {
0036 std::cout << "setting Quality cuts for isolation PFCands:" << std::endl;
0037 }
0038 edm::ParameterSet cfg_isolationQualityCuts = cfg.getParameter<edm::ParameterSet>("isolationQualityCuts");
0039 isolationQualityCutsDzCutDisabled_ = readL1PFTauQualityCuts(cfg_isolationQualityCuts, "disabled", debug_);
0040 isolationQualityCutsDzCutEnabledPrimary_ =
0041 readL1PFTauQualityCuts(cfg_isolationQualityCuts, "enabled_primary", debug_);
0042 isolationQualityCutsDzCutEnabledPileup_ = readL1PFTauQualityCuts(cfg_isolationQualityCuts, "enabled_pileup", debug_);
0043 }
0044
0045 void L1HPSPFTauBuilder::reset() {
0046 signalConeSize_ = 0.;
0047 signalConeSize2_ = 0.;
0048
0049 l1PFCandProductID_ = edm::ProductID();
0050 isPFCandSeeded_ = false;
0051 l1PFCandSeed_ = l1t::PFCandidateRef();
0052 isJetSeeded_ = false;
0053 l1JetSeed_ = reco::CaloJetRef();
0054 l1PFTauSeedEta_ = 0.;
0055 l1PFTauSeedPhi_ = 0.;
0056 l1PFTauSeedZVtx_ = 0.;
0057 sumAllL1PFCandidatesPt_ = 0.;
0058 primaryVertex_ = l1t::VertexWordRef();
0059 l1PFTau_ = l1t::HPSPFTau();
0060
0061 stripP4_ = reco::Particle::LorentzVector(0., 0., 0., 0.);
0062
0063 signalAllL1PFCandidates_.clear();
0064 signalChargedHadrons_.clear();
0065 signalElectrons_.clear();
0066 signalNeutralHadrons_.clear();
0067 signalPhotons_.clear();
0068 signalMuons_.clear();
0069
0070 stripAllL1PFCandidates_.clear();
0071 stripElectrons_.clear();
0072 stripPhotons_.clear();
0073
0074 isoAllL1PFCandidates_.clear();
0075 isoChargedHadrons_.clear();
0076 isoElectrons_.clear();
0077 isoNeutralHadrons_.clear();
0078 isoPhotons_.clear();
0079 isoMuons_.clear();
0080
0081 sumAllL1PFCandidates_.clear();
0082 sumChargedHadrons_.clear();
0083 sumElectrons_.clear();
0084 sumNeutralHadrons_.clear();
0085 sumPhotons_.clear();
0086 sumMuons_.clear();
0087
0088 sumChargedIsoPileup_ = 0.;
0089 }
0090
0091 void L1HPSPFTauBuilder::setL1PFCandProductID(const edm::ProductID& l1PFCandProductID) {
0092 l1PFCandProductID_ = l1PFCandProductID;
0093 }
0094
0095 void L1HPSPFTauBuilder::setVertex(const l1t::VertexWordRef& primaryVertex) { primaryVertex_ = primaryVertex; }
0096
0097 void L1HPSPFTauBuilder::setL1PFTauSeed(const l1t::PFCandidateRef& l1PFCandSeed) {
0098 if (debug_) {
0099 std::cout << "<L1HPSPFTauBuilder::setL1PFTauSeed>:" << std::endl;
0100 std::cout << "seeding HPSPFTau with ChargedPFCand:";
0101 printPFCand(std::cout, *l1PFCandSeed, primaryVertex_);
0102 }
0103
0104 l1PFCandSeed_ = l1PFCandSeed;
0105 l1PFTauSeedEta_ = l1PFCandSeed->eta();
0106 l1PFTauSeedPhi_ = l1PFCandSeed->phi();
0107 if (l1PFCandSeed->charge() != 0 && l1PFCandSeed->pfTrack().isNonnull()) {
0108 l1PFTauSeedZVtx_ = l1PFCandSeed->pfTrack()->vertex().z();
0109 isPFCandSeeded_ = true;
0110 }
0111 }
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152 void L1HPSPFTauBuilder::setL1PFTauSeed(const reco::CaloJetRef& l1JetSeed,
0153 const std::vector<l1t::PFCandidateRef>& l1PFCands) {
0154 if (debug_) {
0155 std::cout << "<L1HPSPFTauBuilder::setL1PFTauSeed>:" << std::endl;
0156 std::cout << "seeding HPSPFTau with Jet:";
0157 std::cout << " pT = " << l1JetSeed->pt() << ", eta = " << l1JetSeed->eta() << ", phi = " << l1JetSeed->phi()
0158 << std::endl;
0159 }
0160
0161 l1JetSeed_ = l1JetSeed;
0162 reco::Candidate::LorentzVector l1PFTauSeed_p4;
0163 float l1PFTauSeedZVtx = 0.;
0164 bool l1PFTauSeed_hasVtx = false;
0165 float max_chargedPFCand_pt = -1.;
0166 for (const auto& l1PFCand : l1PFCands) {
0167 double dR = reco::deltaR(l1PFCand->eta(), l1PFCand->phi(), l1JetSeed->eta(), l1JetSeed->phi());
0168 if (dR > 0.4)
0169 continue;
0170 if (l1PFCand->id() == l1t::PFCandidate::ChargedHadron || l1PFCand->id() == l1t::PFCandidate::Electron ||
0171 l1PFCand->id() == l1t::PFCandidate::Photon || l1PFCand->id() == l1t::PFCandidate::Muon) {
0172 l1PFTauSeed_p4 += l1PFCand->p4();
0173 if (l1PFCand->charge() != 0 && l1PFCand->pfTrack().isNonnull() && l1PFCand->pt() > max_chargedPFCand_pt) {
0174 l1PFTauSeedZVtx = l1PFCand->pfTrack()->vertex().z();
0175 l1PFTauSeed_hasVtx = true;
0176 max_chargedPFCand_pt = l1PFCand->pt();
0177 }
0178 }
0179 }
0180 if (l1PFTauSeed_p4.pt() > 1. && l1PFTauSeed_hasVtx) {
0181 l1PFTauSeedEta_ = l1PFTauSeed_p4.eta();
0182 l1PFTauSeedPhi_ = l1PFTauSeed_p4.phi();
0183 l1PFTauSeedZVtx_ = l1PFTauSeedZVtx;
0184 isJetSeeded_ = true;
0185 }
0186 }
0187
0188 void L1HPSPFTauBuilder::addL1PFCandidates(const std::vector<l1t::PFCandidateRef>& l1PFCands) {
0189 if (debug_) {
0190 std::cout << "<L1HPSPFTauBuilder::addL1PFCandidates>:" << std::endl;
0191 }
0192
0193
0194
0195
0196 if (!(isPFCandSeeded_ || isJetSeeded_))
0197 return;
0198
0199 for (const auto& l1PFCand : l1PFCands) {
0200 if (!isWithinIsolationCone(*l1PFCand))
0201 continue;
0202 sumAllL1PFCandidates_.push_back(l1PFCand);
0203 if (l1PFCand->id() == l1t::PFCandidate::ChargedHadron) {
0204 sumChargedHadrons_.push_back(l1PFCand);
0205 } else if (l1PFCand->id() == l1t::PFCandidate::Electron) {
0206 sumElectrons_.push_back(l1PFCand);
0207 } else if (l1PFCand->id() == l1t::PFCandidate::NeutralHadron) {
0208 sumNeutralHadrons_.push_back(l1PFCand);
0209 } else if (l1PFCand->id() == l1t::PFCandidate::Photon) {
0210 sumPhotons_.push_back(l1PFCand);
0211 } else if (l1PFCand->id() == l1t::PFCandidate::Muon) {
0212 sumMuons_.push_back(l1PFCand);
0213 }
0214 }
0215
0216 for (const auto& l1PFCand : sumAllL1PFCandidates_) {
0217 sumAllL1PFCandidatesPt_ += l1PFCand->pt();
0218 }
0219 std::vector<double> emptyV;
0220 std::vector<double> sumAllL1PFCandidatesPt(1);
0221 sumAllL1PFCandidatesPt[0] = sumAllL1PFCandidatesPt_;
0222
0223 signalConeSize_ = signalConeSizeFormula_.evaluate(sumAllL1PFCandidatesPt, emptyV);
0224
0225 if (signalConeSize_ < minSignalConeSize_)
0226 signalConeSize_ = minSignalConeSize_;
0227 if (signalConeSize_ > maxSignalConeSize_)
0228 signalConeSize_ = maxSignalConeSize_;
0229 signalConeSize2_ = signalConeSize_ * signalConeSize_;
0230
0231 for (const auto& l1PFCand : sumAllL1PFCandidates_) {
0232 if (debug_) {
0233 printPFCand(std::cout, *l1PFCand, primaryVertex_);
0234 }
0235
0236 bool isSignalPFCand = false;
0237 bool isStripPFCand = false;
0238 bool isElectron_or_Photon =
0239 l1PFCand->id() == l1t::PFCandidate::Electron || l1PFCand->id() == l1t::PFCandidate::Photon;
0240 bool isChargedHadron = l1PFCand->id() == l1t::PFCandidate::ChargedHadron;
0241 if (isWithinSignalCone(*l1PFCand) && !(isChargedHadron && signalChargedHadrons_.size() > 3)) {
0242 isSignalPFCand = true;
0243 }
0244 if (isElectron_or_Photon && isWithinStrip(*l1PFCand)) {
0245 if (useStrips_) {
0246 isSignalPFCand = true;
0247 }
0248 isStripPFCand = true;
0249 }
0250 bool passesSignalQualityCuts = isSelected(signalQualityCutsDzCutEnabledPrimary_, *l1PFCand, l1PFTauSeedZVtx_);
0251 if (isSignalPFCand && passesSignalQualityCuts) {
0252 signalAllL1PFCandidates_.push_back(l1PFCand);
0253 if (l1PFCand->id() == l1t::PFCandidate::ChargedHadron) {
0254 signalChargedHadrons_.push_back(l1PFCand);
0255 } else if (l1PFCand->id() == l1t::PFCandidate::Electron) {
0256 signalElectrons_.push_back(l1PFCand);
0257 } else if (l1PFCand->id() == l1t::PFCandidate::NeutralHadron) {
0258 signalNeutralHadrons_.push_back(l1PFCand);
0259 } else if (l1PFCand->id() == l1t::PFCandidate::Photon) {
0260 signalPhotons_.push_back(l1PFCand);
0261 } else if (l1PFCand->id() == l1t::PFCandidate::Muon) {
0262 signalMuons_.push_back(l1PFCand);
0263 }
0264 }
0265 if (isStripPFCand && passesSignalQualityCuts) {
0266 stripAllL1PFCandidates_.push_back(l1PFCand);
0267 if (l1PFCand->id() == l1t::PFCandidate::Electron) {
0268 stripElectrons_.push_back(l1PFCand);
0269 stripP4_ += l1PFCand->p4();
0270 } else if (l1PFCand->id() == l1t::PFCandidate::Photon) {
0271 stripPhotons_.push_back(l1PFCand);
0272 stripP4_ += l1PFCand->p4();
0273 } else
0274 assert(0);
0275 }
0276
0277 bool isIsolationPFCand = isWithinIsolationCone(*l1PFCand) && !isSignalPFCand;
0278 bool passesIsolationQualityCuts = isSelected(isolationQualityCutsDzCutEnabledPrimary_, *l1PFCand, l1PFTauSeedZVtx_);
0279 if (isIsolationPFCand && passesIsolationQualityCuts) {
0280 isoAllL1PFCandidates_.push_back(l1PFCand);
0281 if (l1PFCand->id() == l1t::PFCandidate::ChargedHadron) {
0282 isoChargedHadrons_.push_back(l1PFCand);
0283 } else if (l1PFCand->id() == l1t::PFCandidate::Electron) {
0284 isoElectrons_.push_back(l1PFCand);
0285 } else if (l1PFCand->id() == l1t::PFCandidate::NeutralHadron) {
0286 isoNeutralHadrons_.push_back(l1PFCand);
0287 } else if (l1PFCand->id() == l1t::PFCandidate::Photon) {
0288 isoPhotons_.push_back(l1PFCand);
0289 } else if (l1PFCand->id() == l1t::PFCandidate::Muon) {
0290 isoMuons_.push_back(l1PFCand);
0291 }
0292 }
0293
0294 if (debug_) {
0295 std::cout << "dR = " << reco::deltaR(l1PFCand->eta(), l1PFCand->phi(), l1PFTauSeedEta_, l1PFTauSeedPhi_) << ":"
0296 << " isSignalPFCand = " << isSignalPFCand << ", isStripPFCand = " << isStripPFCand
0297 << " (passesSignalQualityCuts = " << passesSignalQualityCuts << "),"
0298 << " isIsolationPFCand = " << isIsolationPFCand
0299 << " (passesIsolationQualityCuts = " << passesIsolationQualityCuts << ")" << std::endl;
0300 }
0301 }
0302
0303 for (const auto& l1PFCand : l1PFCands) {
0304 if (!isWithinIsolationCone(*l1PFCand))
0305 continue;
0306
0307 if (l1PFCand->charge() != 0 && isSelected(isolationQualityCutsDzCutEnabledPileup_, *l1PFCand, l1PFTauSeedZVtx_)) {
0308 sumChargedIsoPileup_ += l1PFCand->pt();
0309 }
0310 }
0311 }
0312
0313
0314
0315 bool L1HPSPFTauBuilder::isWithinSignalCone(const l1t::PFCandidate& l1PFCand) {
0316 if (isPFCandSeeded_ || isJetSeeded_) {
0317 double deltaEta = l1PFCand.eta() - l1PFTauSeedEta_;
0318 double deltaPhi = l1PFCand.phi() - l1PFTauSeedPhi_;
0319 if ((deltaEta * deltaEta + deltaPhi * deltaPhi) < signalConeSize2_)
0320 return true;
0321 }
0322 return false;
0323 }
0324
0325 bool L1HPSPFTauBuilder::isWithinStrip(const l1t::PFCandidate& l1PFCand) {
0326 if (isPFCandSeeded_ || isJetSeeded_) {
0327 double deltaEta = l1PFCand.eta() - l1PFTauSeedEta_;
0328 double deltaPhi = l1PFCand.phi() - l1PFTauSeedPhi_;
0329 if (std::fabs(deltaEta) < stripSizeEta_ && std::fabs(deltaPhi) < stripSizePhi_)
0330 return true;
0331 }
0332 return false;
0333 }
0334
0335 bool L1HPSPFTauBuilder::isWithinIsolationCone(const l1t::PFCandidate& l1PFCand) {
0336 double deltaEta = l1PFCand.eta() - l1PFTauSeedEta_;
0337 double deltaPhi = l1PFCand.phi() - l1PFTauSeedPhi_;
0338 if ((deltaEta * deltaEta + deltaPhi * deltaPhi) < isolationConeSize2_)
0339 return true;
0340 else
0341 return false;
0342 }
0343
0344 void L1HPSPFTauBuilder::buildL1PFTau() {
0345 reco::Particle::LorentzVector l1PFTau_p4;
0346 for (const auto& l1PFCand : signalAllL1PFCandidates_) {
0347 if (l1PFCand->id() == l1t::PFCandidate::ChargedHadron || l1PFCand->id() == l1t::PFCandidate::Electron ||
0348 l1PFCand->id() == l1t::PFCandidate::Photon) {
0349 l1PFTau_p4 += l1PFCand->p4();
0350 if (l1PFCand->charge() != 0 &&
0351 (l1PFTau_.leadChargedPFCand().isNull() || l1PFCand->pt() > l1PFTau_.leadChargedPFCand()->pt())) {
0352 l1PFTau_.setLeadChargedPFCand(l1PFCand);
0353 }
0354 }
0355 }
0356 if (l1PFTau_.leadChargedPFCand().isNonnull() && l1PFTau_.leadChargedPFCand()->pfTrack().isNonnull()) {
0357 l1PFTau_.setZ(l1PFTau_.leadChargedPFCand()->pfTrack()->vertex().z());
0358
0359 l1PFTau_.setP4(l1PFTau_p4);
0360
0361 l1PFTau_.setSeedChargedPFCand(l1PFCandSeed_);
0362 l1PFTau_.setSeedJet(l1JetSeed_);
0363
0364 l1PFTau_.setSignalAllL1PFCandidates(convertToRefVector(signalAllL1PFCandidates_));
0365 l1PFTau_.setSignalChargedHadrons(convertToRefVector(signalChargedHadrons_));
0366 l1PFTau_.setSignalElectrons(convertToRefVector(signalElectrons_));
0367 l1PFTau_.setSignalNeutralHadrons(convertToRefVector(signalNeutralHadrons_));
0368 l1PFTau_.setSignalPhotons(convertToRefVector(signalPhotons_));
0369 l1PFTau_.setSignalMuons(convertToRefVector(signalMuons_));
0370
0371 l1PFTau_.setStripAllL1PFCandidates(convertToRefVector(stripAllL1PFCandidates_));
0372 l1PFTau_.setStripElectrons(convertToRefVector(stripElectrons_));
0373 l1PFTau_.setStripPhotons(convertToRefVector(stripPhotons_));
0374
0375 l1PFTau_.setIsoAllL1PFCandidates(convertToRefVector(isoAllL1PFCandidates_));
0376 l1PFTau_.setIsoChargedHadrons(convertToRefVector(isoChargedHadrons_));
0377 l1PFTau_.setIsoElectrons(convertToRefVector(isoElectrons_));
0378 l1PFTau_.setIsoNeutralHadrons(convertToRefVector(isoNeutralHadrons_));
0379 l1PFTau_.setIsoPhotons(convertToRefVector(isoPhotons_));
0380 l1PFTau_.setIsoMuons(convertToRefVector(isoMuons_));
0381
0382 l1PFTau_.setSumAllL1PFCandidates(convertToRefVector(sumAllL1PFCandidates_));
0383 l1PFTau_.setSumChargedHadrons(convertToRefVector(sumChargedHadrons_));
0384 l1PFTau_.setSumElectrons(convertToRefVector(sumElectrons_));
0385 l1PFTau_.setSumNeutralHadrons(convertToRefVector(sumNeutralHadrons_));
0386 l1PFTau_.setSumPhotons(convertToRefVector(sumPhotons_));
0387 l1PFTau_.setSumMuons(convertToRefVector(sumMuons_));
0388
0389 l1PFTau_.setPrimaryVertex(primaryVertex_);
0390
0391 if (l1PFTau_.signalChargedHadrons().size() > 1) {
0392 if (stripP4_.pt() < 5.)
0393 l1PFTau_.setTauType(l1t::HPSPFTau::kThreeProng0Pi0);
0394 else
0395 l1PFTau_.setTauType(l1t::HPSPFTau::kThreeProng1Pi0);
0396 } else {
0397 if (stripP4_.pt() < 5.)
0398 l1PFTau_.setTauType(l1t::HPSPFTau::kOneProng0Pi0);
0399 else
0400 l1PFTau_.setTauType(l1t::HPSPFTau::kOneProng1Pi0);
0401 }
0402
0403 l1PFTau_.setStripP4(stripP4_);
0404
0405 l1PFTau_.setSumAllL1PFCandidatesPt(sumAllL1PFCandidatesPt_);
0406 l1PFTau_.setSignalConeSize(signalConeSize_);
0407 l1PFTau_.setisolationConeSize(isolationConeSize_);
0408
0409 double sumChargedIso = 0.;
0410 double sumNeutralIso = 0.;
0411 for (const auto& l1PFCand : isoAllL1PFCandidates_) {
0412 if (l1PFCand->charge() != 0) {
0413 sumChargedIso += l1PFCand->pt();
0414 } else if (l1PFCand->id() == l1t::PFCandidate::Photon) {
0415 sumNeutralIso += l1PFCand->pt();
0416 }
0417 }
0418 l1PFTau_.setSumChargedIso(sumChargedIso);
0419 l1PFTau_.setSumNeutralIso(sumNeutralIso);
0420 const double weightNeutralIso = 1.;
0421 const double offsetNeutralIso = 0.;
0422 l1PFTau_.setSumCombinedIso(sumChargedIso + weightNeutralIso * (sumNeutralIso - offsetNeutralIso));
0423 l1PFTau_.setSumChargedIsoPileup(sumChargedIsoPileup_);
0424
0425 if (l1PFTau_.sumChargedIso() < 20.0) {
0426 l1PFTau_.setPassVLooseIso(true);
0427 }
0428 if (l1PFTau_.sumChargedIso() < 10.0) {
0429 l1PFTau_.setPassLooseIso(true);
0430 }
0431 if (l1PFTau_.sumChargedIso() < 5.0) {
0432 l1PFTau_.setPassMediumIso(true);
0433 }
0434 if (l1PFTau_.sumChargedIso() < 2.5) {
0435 l1PFTau_.setPassTightIso(true);
0436 }
0437
0438 if (l1PFTau_p4.pt() != 0) {
0439 if (l1PFTau_.sumChargedIso() / l1PFTau_p4.pt() < 0.40) {
0440 l1PFTau_.setPassVLooseRelIso(true);
0441 }
0442 if (l1PFTau_.sumChargedIso() / l1PFTau_p4.pt() < 0.20) {
0443 l1PFTau_.setPassLooseRelIso(true);
0444 }
0445 if (l1PFTau_.sumChargedIso() / l1PFTau_p4.pt() < 0.10) {
0446 l1PFTau_.setPassMediumRelIso(true);
0447 }
0448 if (l1PFTau_.sumChargedIso() / l1PFTau_p4.pt() < 0.05) {
0449 l1PFTau_.setPassTightRelIso(true);
0450 }
0451 }
0452 }
0453 }
0454
0455 l1t::PFCandidateRefVector L1HPSPFTauBuilder::convertToRefVector(const std::vector<l1t::PFCandidateRef>& l1PFCands) {
0456 l1t::PFCandidateRefVector l1PFCandsRefVector(l1PFCandProductID_);
0457 for (const auto& l1PFCand : l1PFCands) {
0458 l1PFCandsRefVector.push_back(l1PFCand);
0459 }
0460 return l1PFCandsRefVector;
0461 }